Characteristics of ionospheric anomalies before the earthquake in Indonesia on August 5，2018
-
摘要: 基于2018年8月5日印度尼西亚地震震中（116.45°E，8.33°S）附近的澳大利亚达尔文站和位于震中磁力线共轭区域的武汉站的电离层测高仪观测数据，以及相同区域全球卫星导航系统接收机的观测数据，对该地震震前的电离层异常扰动特征进行了分析。结果显示，震前同时观测到了电离层F2层临界频率（foF2）和总电子含量（TEC）时间序列的异常。基于岩石圈-大气层-电离层直流电场耦合模式和美国国家大气研究中心的热层-电离层环流耦合模式（NCAR/HAO TIEGCM）对震前震中及其对应磁力线共轭点区域的电子密度异常进行了模拟，模拟结果表明，在地震异常电场的作用下，地震区域及其对应半球的磁共轭区域的TEC和电离层F2层峰值电子密度NmF2发生了明显扰动。
-
关键词:
- 岩石圈-大气层-电离层耦合 /
- foF2 /
- TEC /
- TIEGCM模拟
Abstract: Based on the observations of the ionosondes at the Darwin station in Australia near the epicenter （116.45°E, 8.33°S） of the earthquake in Indonesia on August 5, 2018 and at the Wuhan station in the magnetically conjugated area of the epicenter, and the observations of the GNSS receiver in the corresponding regions, this paper analyzed the characteristics and the mechanism of the ionospheric anomaly disturbance. The results show that both the critical frequency （ foF2） and time series of the ionospheric total electron content （TEC） were observed to be abnormal before the earthquake. Based on the lithosphere-atmosphere-ionosphere DC electric field coupling model and thermosphere-ionosphere-electrodynamics general circulation model （TIEGCM） from the National Center for Atmospheric Research （NCAR/HAO）, the electron density anomalies of the seismic area and the magnetically conjugated area before the earthquake are simulated. The simulation results show that TEC and maximum electron density NmF2 of the seismic area and the magnetically conjugated area are significantly disturbed by the abnormal electric field before the earthquake.-
Keywords:
- lithosphere-atmosphere-ionosphere coupling /
- foF2 /
- TEC /
- TIEGCM simulation
-
引言
龙门山断裂带位于四川盆地西部，青藏高原东缘，南起四川天全，总体以北东40°—50°方向延伸至广元，全长绵延约500 km，宽约30—40 km (李勇等，2009a)(图 1).整个龙门山断裂带主要由3条北东走向的断裂组成：汶川—茂县断裂 (后山断裂)、北川—映秀断裂 (中央断裂) 和江油—灌县断裂 (前山断裂); 3条断裂呈大致平行延伸，经过长期且多阶段的构造演化作用，该区地质构造也愈加复杂 (徐杰等，2010；李志伟等，2011).作为稳定的扬子地块与活动的松潘—甘孜地块之间的过渡带，龙门山断裂带不仅在地形上显示出剧烈的变化，其两侧的地壳厚度 (含沉积盖层，下同) 也存在着相似的变化模式，即地壳厚度从松潘—甘孜地块之下的57—58 km减少至四川盆地之下的43—44 km，地形和地壳厚度的巨大差异表明该区发生了强烈的构造变形 (李志伟等，2011).
![]()
图 1 龙门山区域主要活动构造及深地震测深剖面位置图F1:江油—灌县断裂；F2:北川—映秀断裂；F3:汶川—茂县断裂Figure 1. Main active tectonic settings and location of the deep seismic sounding profile in Longmenshan regionF1: Jiangyou-Guanxian fault; F2: Beichuan-Yingxiu fault; F3: Wenchuan-Maoxian fault据历史地震资料记载，中等强度及以上的地震主要集中在龙门山断裂带的南段及中南段，且龙门山断裂带的3条断裂上均已发生过M6以上的中等强度地震. 李勇等 (2009b)根据历史地震记录分析推测，龙门山断裂带近几十年来已进入地震高发期，3条主断裂均具备发生更大强度地震的潜在危险.
龙门山断裂带长期以来一直是地球科学家们关注的焦点. 2008年5月12日，位于龙门山中央的映秀—北川断裂上的汶川县映秀镇发生MS8.0地震 (杨智娴等，2012).此次地震引起了世界范围内地球科学界的极大关注.在此之前，关于龙门山地区的地下地质构造，地学界已开展了大量的研究工作. 蒋福珍和方剑 (2001)使用重力场分离和密度反演方法对康滇地区的地壳构造进行了研究; 李勇等 (2005)给出了龙门山断裂带的均衡重力异常及其对青藏高原东缘山脉地壳隆升的约束; 王谦身等 (2009)研究了四川中西部地区的地壳结构与重力均衡; 李立和金国元 (1987)以及孙洁等 (2003)使用大地电磁测深方法得到了攀西裂谷带及龙门山断裂带地壳上地幔的电性结构；张先等 (1996)获得了四川盆地及其西部边缘震区的居里等温面；张培震等 (2008)和徐锡伟等 (2008)对龙门山断裂带的构造演化过程以及地质成因进行了深入的探讨；刘建华等 (1989)，Huang等 (2002)，Wang等 (2003)和雷建设等 (2009)利用大量天然地震走时资料确定了该区的三维速度结构；刘启元等 (2009)和李志伟等 (2011)利用地震层析成像方法和远震接收函数反演了地下速度结构、地壳厚度分布、平均泊松比和波速比，并且通过接收函数的偏振分析研究了龙门山断裂带两侧地壳的各向异性差异 (楼海等，2008；Zhang et al，2009；Robert et al，2010；郑勇等，2013)；刘启元等 (2009)和邓文泽等 (2014)通过对汶川地震余震序列的研究确定了龙门山断裂带中上地壳速度结构和断裂构造，从而探讨了该区域深部构造和汶川地震的形成原因 (李志伟等，2011).
自20世纪80年代中期起，龙门山地区实施了一系列壳幔构造的探测计划，布设了多条人工地震剖面，仅横穿龙门山断裂带中段的测线就有两条.但是，由于当时观测系统不完善，地震仪器落后，一直未能得到翔实可靠的地下深部构造数据资料，难以对龙门山地区的地壳结构进行细致的研究.此外，李志伟等 (2011)关于该区天然地震的观测研究显示，低速层分布和形态分布等方面存在着不一致.而汶川地震的发生，推翻了之前对于龙门山断裂带的一些认知，提出了一系列新的科学问题.为了更加深入细致地研究该区的壳幔精细速度结构，中国地震局地球物理勘探中心布设了一条横穿龙门山断裂带的人工深地震测深剖面.
本文拟对横穿龙门山断裂带的这条人工深地震测深剖面数据进行反演计算和解释.首先用每一炮的多个震相反演一维速度模型，利用多炮一维模型插值出一个较为粗略的二维初始速度模型；然后使用射线追踪法计算反射、折射波的理论到时，并且根据理论走时与观测走时的拟合程度 (试错法)，调整不同层位的速度和深度，得到一个较为接近真实地壳上地幔顶部的二维速度模型.在试错法建立的速度模型基础上，本文提出一种对深地震测深剖面数据进行联合迭代反演 (simultaneous iterative reconstruction tomography，简写为SIRT) 的方法, 并使用该方法对横跨龙门山断裂带中段的这条深地震宽角反射/折射剖面进行反演，以期获取该剖面的二维地壳速度结构.
1. 人工地震测深剖面位置与观测系统
2008年汶川MS8.0地震发生后，中国地震局地球物理勘探中心布设了一条横穿龙门山断裂带的人工深地震测深剖面.该剖面东起四川遂宁 (105°30′02″E，30°36′51″N)，西至阿坝 (101°48′49″E，32°55′08″N)，全长约500 km (图 1).由于整条测线横穿四川盆地、龙门山断裂带和川西北高原，因此对该测线的分析研究能够加深对该区域速度结构特征和深部构造的认识，推动对龙门山断裂带发展演化过程的研究.
该测线横穿龙门山断裂带中段，途经绵竹、茂县等极震区，向西北经过黑水，最终到达阿坝州.整条测线布设12个炮点，炸药量为600—2800 kg (嘉世旭等，2014)，采用井下爆破的方式，井深范围为15—61 m.龙门山断裂带及其附近的炮点较密集，最小炮检距约为10 km，而在远离龙门山断裂带的四川盆地和高原地区炮检距最大约为110 km.
图 2给出了测线所采用的宽角反射/折射观测系统，通过该系统可以实现对地下介质进行多次覆盖和交叉采样，保证壳内波组可以连续对比追踪 (张恩会等，2013)，从而压制部分干扰波，提高地震记录的信噪比.本文选取12炮中震相清晰、容易识别且较为典型的6炮数据进行处理分析，具体列于表 1.
![]()
图 2 龙门山断裂带中段及两侧宽角反射/折射观测系统及检波器高程Figure 2. Layout of wide-angle reflection/ refraction observatory system of the middle Longmenshan fault zone and elevation of receivers along the profile表 1 炮点参数一览表Table 1. Shot parameters list炮点
编号桩号
/km炮点坐标 药量
/kg平均井深
/m高程
/m地理位置 东经/° 北纬/° SP1 113.478 105°30′02″ 30°36′51″ 1800 61 280 四川省遂宁市新桥镇黄连沱村 SP5 255.298 104°18′26″ 31°23′14″ 800 20 486 四川省绵竹县绵远乡双狮村 SP8 290.191 104°03′13″ 31°37′07″ 1000 10 1310 四川省绵竹县清平乡长河坝村 SP9 320.638 103°48′10″ 31°47′26″ 600 26.5 1 712 四川省茂县沟口乡五家村 SP12 467.784 102°24′14″ 32°24′44″ 2688 20.8 3596 四川省红原县龙日坝乡 SP13 545.360 101°48′49″ 32°55′08″ 2814 42.5 3468 四川省阿坝县麦昆乡 2. 主要震相及其特征
整条测线的原始数据分为12个组段，首先，从其中挑选出同一炮点的数据，重新分组形成共炮点数据文件；然后，对每一炮的地震数据进行分析，识别出清晰可靠的震相；最后，进行走时拾取，得到观测走时数据.
识别震相时，不仅要对多炮的对应波组进行对比，而且要分析震相的运动学和动力学特征.震相识别是所有工作的基础，它不仅是构建初始速度模型的依据，同时还对反演的可靠性有决定作用.因此，对于不太清晰的或者无法确定的震相应采取宁缺毋滥的原则予以剔除.最终分析识别得到的主要震相为：结晶基底折射波Pg，莫霍面反射波PM，折射波Pn和地壳内反射波Pi(图 3).测线以龙门山断裂带为界，东侧为扬子地块 (四川盆地)，西侧为松潘—甘孜地块 (川西北高原)，两大地块之间存在着不同的地下地质构造和速度结构，因而使得两侧的地震波震相、走时和振幅也存在较大差异.
![]()
3 深地震测深地震记录截面及走时拾取图 (折合速度为6.0 km/s)黑色短线表示拾取的走时. (a) SP1; (b) SP5; (c) SP83. Seismic record sections and picked travel times for the shots SP1 (a), SP5 (b) and SP8 (c)Reduced travel time is 6.0 km/s, and black short lines represent picked travel times![]()
3 深地震测深地震记录截面及走时拾取图 (折合速度为6.0 km/s)黑色短线表示拾取的走时. (d) SP9; (e) SP12; (f) SP133. Seismic record sections and picked travel times for the shots SP9 (d), SP12 (e) and SP13 (f)Reduced travel time is 6.0 km/s, and black short lines represent picked travel times2.1 结晶基底折射波Pg震相
Pg震相是地壳内结晶基底的折射波.测线东侧 (100—270 km桩号段) 追踪区间为10—100 km，以清晰可靠的初至波出现 (图 3).在炮点附近视速度迅速增加，远离炮点时视速度逐渐趋于稳定，大约为5.8—6.2 km/s，折合到时约为0—1.8 s，且其能量随着炮检距的增大而逐渐减弱.这些特征体现了扬子地块稳定的结构和较厚的沉积层.测线中部 (270—460 km桩号段) 位于龙门山断裂带，复杂的地下地质构造使得Pg波组的追踪范围大大缩小，仅在炮点附近可以追踪到，折合到时降至约0.1—0.4 s.测线西侧 (460—600 km桩号段) 的Pg波折合到时约为0.3—1.0 s，显示了较浅的沉积层厚度和高速的地下速度结构. Pg震相波组的上述特征是由地表沉积盖层厚度和地壳结晶层顶部介质速度结构所决定的，而局部折合到时的超前、滞后一般是由地表局部的隆起和凹陷构造所引起.
2.2 莫霍界面反射波PM震相
作为壳幔分界面莫霍面的反射波，由于上下岩性差异较大，速度差异也较大，PM波在整个可追踪段震相清晰且振幅能量强，自炮检距80 km处出现，最远可达300 km，可连续对比追踪.与其它壳内反射波震相的比较表明，PM波的时距曲线更为陡峭，视速度更大.
2.3 上地幔顶部折射波Pn震相
Pn波是来自上地幔的折射波，能够反映上地幔顶部的速度结构.由于该波组的能量较弱，连续对比追踪比较困难，除了SP8可以在炮点两侧追踪到外，其余炮仅能在单侧追踪. Pn波出现的最小距离为140 km，最远可达300 km，追踪长度最小仅有十几千米，显示了复杂的上地幔地质构造.以龙门山断裂带为界，测线西段川西北高原的折合到时在约190 km处为0；而在东段的四川盆地，其折合到时在约260 km处为0，显示了高速稳定的上地幔速度结构，视速度明显大于川西北高原.
2.4 地壳内反射波Pi震相
Pi波是地壳内各层的反射波组，其中P1和P3分别为上地壳底部C1界面和中地壳底部C3界面的反射波震相.在所选6炮地震记录截面 (图 3) 上可以看出，地壳内反射波Pi震相较弱；川西北高原可识别的波组震相比四川盆地明显要多，说明川西北高原的地下地质结构更为复杂.
位于四川盆地的炮点SP1西侧和SP5东侧的震相显示，有3组反射波震相较明显，为P1，P3和P4波组.这3组反射波出现在炮检距为50—180 km范围内，随着炮检距的增大，折合到时逐渐减小，其时距曲线相对龙门山断裂带西侧较为陡峭，与四川盆地较为简单的层位分布和较高的地下介质速度相对应.位于龙门山断裂带的SP8和SP9炮点，其附近的反射波震相较凌乱，震相追踪距离较短，折合到时不是单纯地随炮检距增大，这反映了龙门山断裂带复杂的地下速度结构和层内界面展布.位于川西北高原的反射波组震相增多，主要的震相有P1，P2，P3，P4和P5，追踪范围为50—270 km.相较于四川盆地，其时距曲线较平缓，说明壳内介质速度比四川盆地要低.
3. 速度模型的构建
构建地下速度模型需要自上而下确定地壳各层的埋深和速度值.地表的速度值利用炮点附近的Pg波走时得到，界面深度及其上覆层介质的平均速度由x2-t2方法反演得到.在此基础上，构造每一层的速度梯度模型，然后通过走时拟合获取各炮点下方的地壳上地幔一维速度模型.
龙门山断裂带构造变形剧烈.横跨该断裂带的地壳和上地幔速度结构存在强烈的横向不均匀性，因此使用简单的一维地壳模型不可能很好地描述真实的地下速度结构.从图 4可以看出，位于该区的SP5的右支 (图 4b)，SP8的两支 (图 4c) 和SP9的左支 (图 4d) 均较难拟合，比较突出的是SP8左支和SP9左支的P1震相，其走时曲线出现凹凸弯曲，表明其下方存在比较大的界面起伏.此外，SP5和SP8炮点的PM震相理论走时与实测走时的残差较大，说明该区下方莫霍面起伏较大. SP1位于龙门山东侧的四川盆地中部，SP12和SP13位于川西北高原，一维拟合即可取得较好的结果，说明这两个地区的地下结构呈水平层状分布.在所有的实测走时中，会出现局部走时点的突跳，这与地下存在的不规则高速体、低速体有关，一维拟合无法处理，需要建立二维地壳速度结构模型进行研究.
![]()
4 单炮基于一维速度模型的走时拟合图. (a) SP1; (b) SP5；(c) SP8小矩形框内曲线表示一维速度模型；实线代表理论走时，红点代表拾取走时4. Computed travel-time data of single shot based on 1-D velocity model(a) SP1; (b) SP5; (c) SP8. The curves in the small frame represent the 1-D velocity model. Solid lines represent theory travel time, red dots represent observed travel time![]()
4 单炮基于一维速度模型的走时拟合图. (d) SP9; (e) SP12；(f) SP13小矩形框内曲线表示一维速度模型；实线代表理论走时，红点代表拾取走时4. Computed travel-time data of single shot based on 1-D velocity model(d) SP9; (e) SP12; (f) SP13. The curves in the small frame represent the 1-D velocity model. Solid lines represent theory travel time, red dots represent observed travel time首先利用每一炮的一维速度模型，使用插值法构造出一个较为粗略的二维初始速度模型；然后应用“剥皮法”，由浅入深地对沉积层、上地壳、中地壳、下地壳和上地幔顶部，采用射线追踪法计算反射、折射波的理论到时 (Červený V, Pšenčík，1984)；最后根据理论走时与观测走时的拟合程度，调整不同层位的速度和深度，即可得到一个较为接近真实地壳上地幔的二维速度模型.
4. 反演计算
由于试错法存在耗时和不能给出定量误差的缺点，仅仅借助试错法，很难得到高精度的地壳速度结构.为此，本文将使用SIRT方法对观测走时进行反演.首先对二维初始模型进行网格化并建立射线路径上任意节点处的速度与周边网格节点速度之间的关系，然后推导出走时相对于网格节点速度的偏导数，最后建立走时反演方程组.
4.1 理论方法
首先使用矩形网格对二维速度模型进行网格化.每一层介质为一个独立的网格单元，单元的左、右边界与研究区域的左、右边界重合，网格的上、下边界将所在层位的上、下两个界面包含在内.每个层位的矩形网格尺度可以不一致，网格的大小取决于二维初始模型速度的纵向和横向变化，速度梯度较大时采用较密的网格，反之，则采用较稀疏的网格.
对之前得到的二维初始速度模型进行横向和纵向线性插值，可以得到用于反演的二维网格化速度模型.每个网格有4个节点，节点速度用vij表示，其中i为沿x方向 (水平方向) 的节点序号，j为沿深度方向的节点序号.取其中一个矩形网格进行分析，网格的4个节点的速度分别为v11, v12, v21, v22(图 5).不同射线可能穿过网格的不同部位，射线在该网格内的长度也不同，如图 5中射线1和射线2.由于正演射线追踪得到的射线可能在某网格内的射线节点过少，常常会出现射线2的情况 (V1，V2，…, Vn代表沿射线的节点)，这可能导致反演时所求的方程组欠定.在这种情况下，需要对网格内的射线点进行线性插值，以保证最终反演的方程组为适定或超定.
![]()
图 5 速度网格单元空心圈表示射线点，红色圆点表示插值点Figure 5. Velocity grid cellOpen circles represent ray points, and red dots represent interpolation points对网格内的射线节点采用反距离平方法进行线性插值.根据射线追踪结果可知，射线2在图 5所示的矩形网格内没有节点，它在网格外的两个节点分别为V1和V2.为了得到射线2在矩形网格内任意一段射线两个节点 (图中任意两个相邻红色圆点, 例如A和B) 的速度参数，首先要在V1和V2的连线上确定射线段两节端点A和B的坐标，然后分别求出每一个节点到矩形4个顶点的距离.设节点A到矩形4个顶点的距离为d1, d2, d3和d4，节点B到矩形4个顶点的距离为e1, e2, e3和e4，4个节点的速度依次为V11, V12, V21, V22(图 5)，那么，如果用4个网格节点速度表示射线节点速度，A和B两节点的速度为
(1) 
(2) 射线路径上任意相邻两点A和B的走时差tB-tA对矩形节点处速度的偏导数为
(3) 式中
其中s1为射线节点A与B之间的路径长度.同样可得∂tB/∂vA和∂tB/∂vB. ∂vA/∂vij和∂vB/∂vij可由式 (1) 和 (2) 对vij求偏导获取.
对于每一条射线，沿射线路径周围节点速度的改变会导致射线终点走时的改变.用Δt1, Δt2, …, Δtn分别表示第1条，第2条，…, 第n条射线的走时增量，可得
(4) 上式中每一条射线的走时是该条射线上所有射线段走时的总和.将方程组 (4) 等号左右两边互换，可用矩阵形式表示为
(5) 此即为用于反演的方程组，其右边为观测走时与理论走时的差，左边系数矩阵的每个元素是理论走时对模型节点速度的偏导数.在完成射线追踪后, 系数矩阵可由式 (1)，(2) 和 (3) 求出；等号左边的矢量 (ΔV11, ΔV12, …, Δvmn)T是模型节点速度的改变量，可以通过反演计算求得.方程组 (5) 的系数矩阵是一个稀疏矩阵，只有那些有射线穿过的网格的相关节点的偏导数不为零，其余节点的偏导数均为零.
4.2 反演计算
大型稀疏矩阵方程组 (5) 可通过许多方法求解，例如，直接解法，雅克比共轭梯度法，不完全柯列斯基 (Cholesky) 共轭梯度法, 预条件共轭梯度法等.这些方法有的要求高内存，有的要求计算速度快，有的可能存在收敛和稳定性问题.本文将采用层析成像理论中的代数迭代法和联合迭代法 (SIRT).利用代数迭代法计算一条射线之后就可以将走时残差所对应的速度变化分配到射线路径上的相关节点；利用SIRT法则是计算完所有射线后将走时残差所对应的速度变化量分配到相关的节点上.由于SIRT法要求模型更新的次数较少，所以比较稳健.本文主要使用SIRT进行迭代反演.
将走时拟合后得到的二维速度模型作为初始反演模型，使用上述反演方法对2668个走时数据 (见表 2) 进行反演，可以进一步缩小理论走时与观测走时之间的误差. Pg波组的到时记录用来约束地壳顶部速度结构，P1-P5反射波组用来约束对应的壳内C1-C5界面，PM和Pn分别用来约束莫霍面和上地幔顶部速度结构.迭代反演进行13次后，所有检波器实测走时与理论走时残差的平方和呈减小趋势，如图 6所示.可见，在反演的初期，残差平方和急剧减小，随后出现震荡，并最终趋于收敛，其中在第11次反演过程中残差平方和最小，为25.7257 s2.
表 2 用于反演的全部震相及其走时数据统计Table 2. Statistics of the phases and travel times used in inversion炮点 Pg P1 P2 P3 P4 P5 PM Pn SP1 31 27 - 26 - - 48 12 SP5 136 83 63 35 56 - 72 6 SP8 105 106 94 58 69 40 63 13 SP9 155 101 89 73 84 37 - 13 SP12 81 73 88 94 91 77 85 11 SP13 65 36 26 55 33 55 61 10 合计 573 426 360 347 359 209 329 65 注：“-”表示无此震相. ![]()
图 6 迭代反演的走时残差平方和与迭代次数变化关系曲线Figure 6. Relationship between the sum of travel time residual squares and iteration number of inversion图 7给出了反演前与反演迭代第11次的对比结果，可以看出，速度模型有了较大的改善. SP1(图 7a) 和SP13(图 7f) 炮点模型反演前实测走时与理论走时的残差较小，拟合效果比较好；反演后速度模型变化不大. SP5(图 7b) 反演前左支走时残差较小，右支P2-PM震相走时残差较大; 反演后右支走时残差得到了压制. SP8(图 7c) 反演前右支走时残差较小，左支残差主要出现在P1和P2震相，理论走时偏小；反演后走时残差变小. SP9(图 7d) 反演前左支P1和P2震相及右支的P3-PM震相残差较大，左支理论走时偏小，右支理论走时偏大；反演后拟合效果变好. SP12(图 7e) 反演前左支P2震相理论走时偏小，P3-PM震相理论走时偏大, 反演后走时残差变小.
![]()
7 基于二维速度模型的走时拟合图. (a) SP1; (b) SP5左图为反演迭代前结果，右图为反演迭代第11次时结果.黑色符号表示理论走时，红色符号表示实测走时，箭头表示炮点位置7. Travel-times fitting based on 2-D velocity model. (a) SP1; (b) SP5Left panels represent the results befor inversion, right panels represent the results after 11 iteration inversions. Black symbols represent theory travel times, red symbols represent observed travel times, arrows represent shot location![]()
7 基于二维速度模型的走时拟合图. (c) SP8; (d) SP9; (e) SP12; (f) SP13左图为反演迭代前结果，右图为反演迭代第11次时结果.黑色符号表示理论走时，红色符号表示实测走时，箭头表示炮点位置7. Travel-times fitting based on 2-D velocity model. (c) SP8; (d) SP9; (e) SP12; (f) SP13Left panels represent the results befor inversion, right panels represent the results after 11 iteration inversions. Black symbols represent theory travel times, red symbols represent observed travel times, arrows represent shot location5. 壳幔速度结构与构造特征
通过反演得到了二维P波速度结构剖面图, 如图 8所示.可以看出，整条测线下方的地壳结构以龙门山为界，两侧的厚度和速度结构均存在较大的差异.测线东段 (100—270 km桩号段) 位于四川盆地，地壳厚度约为43—44 km；测线西段 (460—600 km桩号段) 位于川西北高原，地壳厚度约为57—58 km；测线中段 (270—460 km桩号段) 位于松潘—甘孜构造转换带，强烈的构造活动使得该区界面起伏明显并伴有错断，地壳厚度介于东、西两段之间.
![]()
图 8 宽角反射/折射剖面二维P波速度结构图红色星形表示汶川地震震源在剖面上的投影，白色水平箭头表示块体运动方向F1：江油—灌县断裂；F2：北川—映秀断裂；F3：汶川—茂县断裂Figure 8. 2-D P wave velocity model along the wide-angle reflection/refraction profile The red star represents the projection of hypocenter of the Wenchuan earthquakeThe red star represents the projection of hypocenter of the Wenchuan earthquake on the profile, and the white horizontal arrows represent the direction of block motion F1: Jiangyou-Guanxian fault; F2: Beichuan-Yingxiu fault; F3: Wenchuan-Maoxian fault5.1 沉积盖层
沿着测线自西向东，东、西两段沉积盖层的特征存在巨大差异.川西北高原沉积盖层的厚度为2.5—3.0 km，介质速度为3.60—5.25 km/s；沿测线向东经过龙日坝后，沉积盖层逐渐变薄；而从桩号280 km左右开始，盖层厚度急剧增大；当测线到达德阳市时，盖层厚度达到5.5—6.5 km，并延伸至遂宁，介质速度为3.55—5.34 km/s.龙门山断裂带两侧沉积盖层厚度的差异与扬子板块和松潘—甘孜板块的地质年龄相吻合.
5.2 上地壳结构
上地壳位于沉积盖层G与C1界面之间，由基底回折波Pg和P1反射波组约束控制.川西北高原上地壳厚度为15—19 km，速度为5.94—6.16 km/s，在10 km深处存在约3 km厚的低速层；四川盆地上地壳厚度为12—15 km，速度为5.87—6.18 km/s.在松潘—甘孜构造转换带的东部，上地壳界面自西向东逐渐抬升，层厚变小，层内出现不规则低速块体，速度在横向上存在差异；其西部层厚逐渐加大，界面趋于水平，层上部存在高速块体.
5.3 中地壳结构
中地壳位于C1界面与C3界面之间，由反射波组P2和P3控制约束.西部高原地区中地壳厚度为12—14 km，速度为6.24—6.31 km/s；其间存在一个明显的速度间断面C2，速度由6.16 km/s突变为6.31 km/s，该间断面将中地壳分为两层，且均为高低速互层；龙门山褶皱带下方存在一长150 km厚40 km的低速体.四川盆地中地壳厚度为6—7 km，相较于川西北高原大幅度变薄，速度为6.29—6.40 km/s，其埋深起伏较大，在龙门山前陆盆地其埋深为18—25 km，在成都平原埋深为21—28 km.
5.4 下地壳结构
下地壳位于C3界面与莫霍面之间，由反射波组P4，P5和莫霍面反射波PM共同控制约束.川西北高原下地壳厚22—25 km，速度值偏低，约为6.46—6.91 km/s；四川盆地下地壳厚16—18 km，速度值较高，为6.66—7.20 km/s.下地壳中存在两个速度间断面：一个为C4间断面，该间断面横贯东西，并且其西部的速度跃变值为0.15 km/s，其东部的速度跃变值为0.25 km/s；另一个速度间断面C5仅存在于龙门山构造带西部之下，其速度从6.38 km/s跃变为6.61 km/s.
5.5 莫霍面及上地幔
莫霍面作为壳幔边界，其深度沿测线变化很大，由川西北高原的53—55 km，自西向东经龙门山褶皱带，逐渐抬升至四川盆地的43—44 km.上地幔顶部的速度Pn在横向上存在比较剧烈的变化，川西北高原速度约为7.95 km/s，而四川盆地高达8.10 km/s.
将汶川地震的震源垂直投影至测量剖面上，如图 8红色星形所示，其在地面的投影坐标为 (103.7830°E，31.6879°N)，震源深度为 (15.5±0.3) km，桩号为332.9 km，可见其位于上地壳底部的地层抬升区.
6. 讨论与结论
本文研究的测线横跨西部松潘—甘孜地块、中部龙门山断裂带和东部扬子地块，在地壳厚度、壳内速度结构方面存在着巨大的差异，低速层、高速体的分布特征也有其独特的地质成因.测线剖面显示：沉积盖层在东部四川盆地的分布明显厚于中部褶皱带和西部高原，而中部褶皱带部分地区出现基岩裸露，表明扬子板块的地质年龄要老于松潘—甘孜板块，中部褶皱带处于活跃的地质构造期；地壳内构造转换带两侧的地层分界面近于水平层状分布，转换带西侧的中、下地壳内各存在一个层间速度间断面，而东侧则不存在类似的间断面.对应上、中、下地壳的分界面，西侧比东侧的速度要小，而且自上而下逐渐增大，下地壳底部莫霍面的速度西侧比东侧要低0.3 km/s.在厚度方面，西侧比东侧的上、中、下地壳厚度分别厚3—7 km，6—8 km，6—9 km，可见地壳的增厚主要由被改造的中下地壳介质形变增厚所致 (李志伟等，2011；嘉世旭等，2014)，这有别于蔡学林等 (1999)一文中地壳厚度基本不变的结果.在构造转换带内，存在着薄厚不等的低速层，自西向东有增厚的趋势.川西北高原的松潘—甘孜地块在向东挤压刚性的扬子地块时受到阻挡，使得转换带内的低速层堆积、增厚，中下地壳物质向上抬升，形成了上陡下缓、向北西倾斜的构造，从而使测线东、西两侧的莫霍面产生了较大的高差，地壳厚度由四川盆地的43—44 km增厚至川西北高原的57—58 km.
龙门山断裂带的3条主断裂即前山断裂F1、中央断裂F2和后山断裂F3，分别向下深切结晶地壳，而断裂F1和F3正是自上而下并向西侧倾斜的巨型铲式低速通道的包络 (图 8).这是由于西部松潘—甘孜地块自西向东运动，受到刚性扬子地块的阻挡，沿铲式断裂向上爬升，便在断层上盘距地表 15 km左右深处出现最大剪切应力的极值区，最终应变能以地震的形式得以释放，与杨智娴等 (2012)对汶川MW7.9地震震源位置的测定结果相吻合. 陶玮等 (2011)利用有限元方法对该区的铲式断裂进行模拟，其结果也显示最大剪应力变化的极值区位于深约14—18 km地壳内.
龙门山断裂带由南西向北东延伸约500 km，不同地段的深部构造存在着差异.以往的深地震测深测线主要布测于龙门山南段，且观测系统布设不够完善.横穿龙门山断裂带的这条测线完善了龙门山深地震测深资料，改善了观测系统的布设，从而得到了精细的地下速度结构，同时为其它地质和地球物理研究提供了参考和依据.
本文提出的反演方法是在试错法建立的速度模型的基础上，对试错法的一种补充.试错法在试错的后期残差收敛较慢，新的反演方法大大减少了工作量和时间；试错法无法对走时残差进行量化，本文使用的方法则弥补了试错法的缺点，使误差处于可控的范围内.然而，该方法尚有很多需要改进的地方，例如网格内射线点速度的插值方法，可以考虑使用更多的网格节点进行插值；对射线上的点进行样条插值，以期更精确地表示射线追踪得到的射线.
中国地震局地球物理勘探中心为本文提供了人工地震数据，参与该项目野外施工的工作人员克服困难，取得了高质量的野外地震资料，中国地震局地球物理勘探中心嘉世旭研究员和刘志研究员为初期的资料处理给予了有益的指导帮助，审稿专家提出了诸多改进意见，作者在此一并表示感谢.
-
![]()
图 1 印度尼西亚地震震中、磁共轭点、电离层测高仪站、IGS站的地理位置分布
Figure 1. Geographic locations of Indonesia earthquake epicenter，magnetically conjugated point，ionosondes and IGS stations
![]()
图 2 2018年7月1日到8月12日Dst指数（a），Kp指数（b）和F10.7 （c）的变化
Figure 2. Variation of the indices Dst （a），Kp （b） and F10.7 （c） from July 1 to August 12，2018
![]()
图 3 2018年8月3日至9日武汉站和达尔文站的foF2观测变化
（a） foF2原始数据；（b） foF2标准差
Figure 3. The variations of foF2 observed by Wuhan station and Darwin station from August 3 to 9，2018
（a） Original data of foF2；（b） Standrad deviation of foF2
![]()
图 4 2018年8月3日至9日武汉站和达尔文站的foF2异常扰动观测
红色、绿色曲线分别表示原始foF2数据减去月中值、月均值后再减去两倍标准差的结果，蓝色方框为武汉站和达尔文站的foF2同时呈现异常的时间点
Figure 4. The abnormal disturbances of foF2 observed by the stations Wuhan and Darwin from August 3 to 9，2018
The red curve indicates the result that the original foF2 data minuses the 31-day median and then subtracts twice the standard deviation. The green curve indicates the result that the original foF2 data minuses the 31-day mean and then subtracts twice the standard deviation. The blue box indicates the time point at which the foF2 in Wuhan and Darwin stations are simultaneously anomalous
![]()
图 5 30号卫星观测到各站点的δTEC序列
Figure 5. The δTEC sequences from different stations observed by the satellite No.30
![]()
图 6 30号卫星观测到的13个站点的IPP分布
Figure 6. IPP distribution of 13 sites observed by the satellite No.30
![]()
图 7 地震区域2018年8月4日（年积日216） UTC 03:40电离层底部z＝90 km处的异常水平电场分布
（a） 总电场强度E；（b） 磁南北向电场强度ESN；（c） 磁东西向电场强度EEW
Figure 7. Distribution of the abnormal horizontal electric field at UTC 03:40 on August 4，2018（216 day of the year） at the bottom of the ionosphere （z＝90 km）
（a） Total electric field intensity E；（b） The electric field intensity ESN in magnetic north-south direction；（c） The electric field intensity EEW in magnetic east-west direction
![]()
图 8 基于TIEGCM模型对插入异常电场四小时后地震及磁共轭区域TEC情况的模拟结果
Figure 8. Simulation results of TEC evolution in the seismically active area and magnetically conjugated area four hours after insertion of an abnormal electric field based on the model TIEGCM
-
丁宗华,吴健,孙树计,陈金松,班盼盼. 2010. 汶川大地震前电离层参量的变化特征与分析[J]. 地球物理学报,<bold>53</bold>(1):30–38. doi: 10.3969/j.issn.0001-5733.2010.01.004 Ding Z H,Wu J,Sun S J,Chen J S,Ban P P. 2010. The variation of ionosphere on some days before the Wenchuan earthquake[J]. <italic>Chinese Journal of Geophysics</italic>,<bold>53</bold>(1):30–38 (in Chinese).
刘静,万卫星,黄建平,张学民,赵庶凡,欧阳新艳,泽仁志玛. 2011. 智利8.8级地震的震前电子浓度扰动[J]. 地球物理学报,<bold>54</bold>(11):2717–2725. Liu J,Wan W X,Huang J P,Zhang X M,Zhao S F,Ouyang X Y,Zeren Z M. 2011. Electron density perturbation before Chile <italic>M</italic>8.8 earthquake[J]. <italic>Chinese Journal of Geophysics</italic>,<bold>54</bold>(11):2717–2725 (in Chinese).
马新欣,林湛,陈化然,金红林,焦立果,刘晓灿. 2014. 基于GPS和COSMIC数据分析汶川地震TEC和<italic>N</italic><sub>mF2</sub>扰动[J]. 地球物理学报,<bold>57</bold>(8):2415–2422. doi: 10.6038/cjg20140803 Ma X X,Lin Z,Chen H R,Jin H L,Jiao L G,Liu X C. 2014. Analysis on ionospheric perturbation of TEC and <italic>N</italic><sub>mF2</sub> based on GPS and COSMIC data before and after the Wenchuan earthquake[J]. <italic>Chinese Journal of Geophysics</italic>,<bold>57</bold>(8):2415–2422 (in Chinese).
滕荣荣,赵谊,刘耀炜,马玉川. 2010. 强震地壳逸出氡和电离层异常耦合关系的讨论[J]. 中国地震,<bold>26</bold>(1):60–72. doi: 10.3969/j.issn.1001-4683.2010.01.006 Teng R R,Zhao Y,Liu Y W,Ma Y C. 2010. The study on the strong earthquakes coupling relationship between abnormal radon escaped from crustal and ionosphere[J]. <italic>Earthquake Research in China</italic>,<bold>26</bold>(1):60–72 (in Chinese).
余涛,毛田,王云冈,王劲松. 2009. 汶川特大地震前电离层主要参量变化[J]. 科学通报,<bold>54</bold>(4):493–499. Yu T,Mao T,Wang Y G,Wang J S. 2009. Study of the ionospheric anomaly before the Wenchuan earthquake[J]. <italic>Chinese Science Bulletin</italic>,<bold>54</bold>(6):1080–1086.
张明敏,刘智敏,刘盼,从建锋. 2018. 九寨沟7.0级地震前电离层TEC异常分析[J]. 测绘工程,<bold>27</bold>(12):24–30. Zhang M M,Liu Z M,Liu P,Cong J F. 2018. Analysis of ionospheric TEC anomalies before the Jiuzhaigou <italic>M</italic><sub>S</sub>7.0 earth-quake[J]. <italic>Engineering of Surveying and Mapping</italic>,<bold>27</bold>(12):24–30 (in Chinese).
Carter B A,Kellerman A C,Kane T A,Dyson P L,Norman R,Zhang K. 2013. Ionospheric precursors to large earthquakes:A case study of the 2011 Japanese Tohoku earthquake[J]. <italic>J Atmos Solar-Terr Phys</italic>,<bold>102</bold>:290–297. doi: 10.1016/j.jastp.2013.06.006
Chuo Y J,Liu J Y,Pulinets S A,Chen Y I. 2002. The ionospheric perturbations prior to the Chi-Chi and Chia-Yi earthquakes[J]. <italic>J Geodyn</italic>,<bold>33</bold>(4/5):509–517.
Dabas R S,Das R M,Sharma K,Pillal K G M. 2007. Ionospheric precursors observed over low latitudes during some of the recent major earthquakes[J]. <italic>J Atmos Solar-Terr Phys</italic>,<bold>69</bold>(15):1813–1824. doi: 10.1016/j.jastp.2007.09.005
Freund F T,Kulahci I G,Cyr G,Ling J,Winnick M,Tregloan-Reed J,Freund M M. 2009. Air ionization at rock surfaces and pre-earthquake signals[J]. <italic>J Atmos Solar-Terr Phys</italic>,<bold>71</bold>(17/18):1824–1834.
Hayakawa M,Molchanov O A,NASDA/UEC Team. 2004. Achievements of NASDA’s earthquake remote sensing frontier project[J]. <italic>Terr Atmos Ocean Sci</italic>,<bold>15</bold>(3):311–327. doi: 10.3319/TAO.2004.15.3.311(EP)
Heinicke J,Koch U,Martinelli G. 1995. CO<sub>2</sub> and radon measurements in the Vogtland area (Germany):A contribution to earthquake prediction research[J]. <italic>Geophys Res Lett</italic>,<bold>22</bold>(7):771–774. doi: 10.1029/94GL03074
Kong J,Yao Y B,Zhou C,Liu Y,Zhai C Z,Wang Z M,Liu L. 2018. Tridimensional reconstruction of the co-seismic ionospheric disturbance around the time of 2015 Nepal earthquake[J]. <italic>J Geodyn</italic>,<bold>92</bold>(11):1255–1266. doi: 10.1007/s00190-018-1117-3
Leonard R S,Barnes R A. 1965. Observation of ionospheric disturbances following the Alaska earthquake[J]. <italic>J Geophys Res</italic>,<bold>70</bold>(5):1250–1253.
Liperovsky V A,Pokhotelov O A,Meister C V,Liperovskaya E V. 2008. Physical models of coupling in the lithosphere-atmosphere-ionosphere system before earthquakes[J]. <italic>Geomagn Aeron</italic>,<bold>48</bold>(6):795–806. doi: 10.1134/S0016793208060133
Liu J Y,Chen Y I,Pulinets S A,Tsai Y B,Chuo Y J. 2000. Seismo-ionospheric signatures prior to <italic>M</italic>≥6.0 Taiwan earthquakes[J]. <italic>Geophys Res Lett</italic>,<bold>27</bold>(19):3113–3116. doi: 10.1029/2000GL011395
Liu J Y,Chuo Y J,Shan S J,Tsai Y B,Chen Y I,Pulinets S A,Yu S B. 2004. Pre-earthquake ionospheric anomalies registered by continuous GPS TEC measurements[J]. <italic>Annal Geophys</italic>,<bold>22</bold>(5):1585–1593. doi: 10.5194/angeo-22-1585-2004
Oyama K I,Kakinami Y,Liu J Y,Kamogawa M,Kodama T. 2008. Reduction of electron temperature in low-latitude ionosphere at 600 km before and after large earthquakes[J]. <italic>J Geophys Res</italic>:<italic>Space Phys</italic>,<bold>113</bold>:A11317.
Popov K V,Liperovsky V A,Meister C V,Biagi P F,Liperovskaya E V,Silina A S. 2004. On ionospheric precursors of earthquakes in scales of 2−3 h[J]. <italic>Phys Chem Earth</italic>,<bold>29</bold>(4/5/6/7/8/9):529–535.
Pulinets S A,Boyarchuk K A,Hegai V V,Kim V P,Lomonosov A M. 2000. Quasielectrostatic model of atmosphere-thermosphere-ionosphere coupling[J]. <italic>Adv Space Res</italic>,<bold>26</bold>(8):1209–1218. doi: 10.1016/S0273-1177(99)01223-5
Pulinets S A,Legen’ka A D. 2003. Spatial-temporal characteristics of large scale disturbances of electron density observed in the ionospheric F-region before strong earthquakes[J]. <italic>Cosmic Res</italic>,<bold>41</bold>(3):221–229. doi: 10.1023/A:1024046814173
Pulinets S A,Legen’ka A D,Gaivoronskaya T V,Depuev V K. 2003. Main phenomenological features of ionospheric precursors of strong earthquakes[J]. <italic>J Atmos Solar-Terr Phys</italic>,<bold>65</bold>(16/18):1337–1347.
Pulinets S. 2004. Ionospheric precursors of earthquakes:Recent advances in theory and practical applications[J]. <italic>Terr Atmos Ocean Sci</italic>,<bold>15</bold>(3):413–435. doi: 10.3319/TAO.2004.15.3.413(EP)
Pulinets S,Ouzounov D. 2011. Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) model:An unified concept for earthquake precursors validation[J]. <italic>J Asian Earth Sci</italic>,<bold>41</bold>(4/5):371–382.
Pulinets S,Davidenko D. 2014. Ionospheric precursors of earthquakes and Global Electric Circuit[J]. <italic>Adv Space Res</italic>,<bold>53</bold>(5):709–723. doi: 10.1016/j.asr.2013.12.035
Richmond A D,Ridley E C,Roble R G. 1992. A thermosphere/ionosphere general circulation model with coupled electrodynamics[J]. <italic>Geophys Res Lett</italic>,<bold>19</bold>(6):601–604. doi: 10.1029/92GL00401
Richmond A D. 1995. Ionospheric electrodynamics using magnetic apex coordinates[J]. <italic>J Geomagn Geoelectr</italic>,<bold>47</bold>(2):191–212. doi: 10.5636/jgg.47.191
Sharma K,Dabas R S,Sarkar S K,Das R M,Ravindran S,Gwal A K. 2010. Anomalous enhancement of ionospheric F<sub>2</sub> layer critical frequency and total electron content over low latitudes before three recent major earthquakes in China[J]. <italic>J Geophys Res</italic>:<italic>Space Phys</italic>,<bold>115</bold>(A11):A11313.
Sorokin V M,Chmyrev V M,Yaschenko A K. 2005. Theoretical model of DC electric field formation in the ionosphere stimulated by seismic activity[J]. <italic>J Atmos Solar-Terr Phys</italic>,<bold>67</bold>(14):1259–1268. doi: 10.1016/j.jastp.2005.07.013
Sorokin V M,Yaschenko A K,Hayakawa M. 2007. A perturbation of DC electric field caused by light ion adhesion to aerosols during the growth in seismic-related atmospheric radioactivity[J]. <italic>Nat Hazards Earth Syst Sci</italic>,<bold>7</bold>(1):155–163. doi: 10.5194/nhess-7-155-2007
Tariq M A,Shah M,Hernández-Pajares M,Iqbal T. 2019. Pre-earthquake ionospheric anomalies before three major earthquakes by GPS-TEC and GIM-TEC data during 2015−2017[J]. <italic>Adv Space Res</italic>,<bold>63</bold>(7):2088–2099. doi: 10.1016/j.asr.2018.12.028
Virk H S,Singh B. 1994. Radon recording of Uttarkashi earthquake[J]. <italic>Geophys Res Lett</italic>,<bold>21</bold>(8):737–740. doi: 10.1029/94GL00310
Zhao B Q,Wang M,Yu T,Wan W X,Lei J H,Liu L B,Ning B Q. 2008. Is an unusual large enhancement of ionospheric electron density linked with the 2008 great Wenchuan earthquake?[J]. <italic>J Geophys Res</italic>:<italic>Space Phys</italic>,<bold>113</bold>:A11304.
Zhou C,Liu Y,Zhao S F,Liu J,Zhang X M,Huang J P,Shen X H,Ni B B,Zhao Z Y. 2017. An electric field penetration model for seismo-ionospheric research[J]. <italic>Adv Space Res</italic>,<bold>60</bold>(10):2217–2232. doi: 10.1016/j.asr.2017.08.007
-
期刊类型引用(2)
1. Yingkang LI，Jianwei GAO，Jian HAN，Yu'e YANG. Geophysical evidence for thrusting of crustal materials from orogenic belts over both sides of the Yangtze Block and its geological significance. Science China(Earth Sciences). 2019(05): 812-831 . 
必应学术
2. 李英康，高建伟，韩健，杨予鄂. 扬子块体两侧造山带地壳推覆的地球物理证据及其地质意义. 中国科学:地球科学. 2019(04): 687-705 . 百度学术
其他类型引用(1)
下载:
计量
- 文章访问数: 938
- HTML全文浏览量: 404
- PDF下载量: 49
- 被引次数: 3
