利用中国GPS站网对地震波引发的大尺度电离层扰动的观测

郝永强， 李泉翰， 郭建广， 张效信， 杨光林， 张东和， 肖佐

1 北京大学地球与空间科学学院, 北京 100871 2 国家空间天气监测预警中心， 北京 100081

0 引言

在地面和低层大气中发生的一些现象，如地震、海啸、火山等，其释放的巨大能量以声波或重力波的形式向高层大气传播，导致电离层的扰动.基于近年来的大量研究，人们认识到地震引发电离层扰动主要有两种途径(Astafyeva， 2019)：一是由于震中断层导致的地面或海面起伏，驱动产生大气中的声波和重力波，二是随着地震波(主要是瑞利波)的传播，在地震波所到之处激发大气中的次声波.前一种大气波动导致的电离层扰动由震中向四周传播，典型速度为几百到1000 m·s-1，这是由声波或重力波在大气和电离层中的传播速度决定的；后一种方式激发的电离层扰动跟随瑞利波的传播，表现出的速度与瑞利波相同，约为3～4 km·s-1.

对地震引起的电离层扰动的研究开始于20世纪60年代(Bolt, 1964; Davies and Baker, 1965)，近年来不断发展的观测技术和观测网络帮助人们形成了更全面的认识.2008年汶川地震后，利用中国区域的观测台网开展了一些研究，并发展了电磁监测试验卫星(Zhao and Hao, 2015; Huang et al., 2018; Yan et al., 2018; Wang et al., 2019).2011年3月11日日本大地震后，地震发生后的电离层扰动被多种观测设备记录下来，如电离层高频多普勒(Chum et al., 2012; Hao et al., 2012)、地磁场(Hao et al., 2013; Yen et al., 2015)、测高仪(Maruyama and Shinagawa, 2014)、GPS接收机测量的电离层总电子含量(total electron content, TEC)等(Astafyeva et al., 2011; Liu et al., 2011; Rolland et al., 2011a; Komjathy et al., 2012; Occhipinti et al., 2013)等.综合这些观测发现，地震引起的大气扰动在从地面向电离层传播的过程中，依次影响了电离层电流/磁场和电子密度(Liu et al., 2016).

此次地震在震中附近(近场)引发的电离层扰动，被日本区域密集的GPS接收机网络清晰地记录了下来.分析GPS TEC扰动发生的时间与其距震中的距离之间的关系，可发现扰动的传播速度主要有两个分支，分别对应前述的震中断层和瑞利波两种扰动激发的来源 (如Liu et al., 2011; Galvan et al., 2012; Kakinami et al., 2013; Rolland et al., 2011a).这次地震还引发了强烈的海啸，海面的起伏导致了以震中为中心的同心环状向外传播的电离层扰动，以及在震中上空出现低电子密度区等有趣的特征(Tsugawa et al., 2011; Kakinami et al., 2012; Astafyeva, et al., 2013; Zettergren et al., 2017; 唐龙等, 2017a, b).密集的GPS TEC数据清楚地展示了这些现象的观测特征.

这次海底地震引发的海啸横越了太平洋，使得地震的能量波及到远离震中(远场)的区域.当海啸波到达美国西海岸时，其产生的大气重力波扰动了上空的电离层，被美国的GPS观测网记录下来(Azeem et al., 2017).然而在陆地上，地震的能量主要以瑞利波的形式传播到远场区域.在地震波传播的同时，地表震动通过大气耦合至电离层高度，所以可观测到电离层扰动与瑞利波同步地以3～4 km·s-1的水平速度传播到距震中数千公里之外(Chum et al., 2012; Hao et al., 2012, 2013; Liu et al., 2016).这些研究主要基于为数不多的观测台站，能够确定扰动传播到远场的大致速度，但却无法获知电离层扰动在远场的大尺度空间结构.

要对远场扰动的大尺度空间结构进行研究，需要观测站网与震中的相对位置合适，并具有较高的台站密度.对于2011年3月11日的日本大地震，中国区域内的GPS接收机网络提供了宝贵的观测机会，本文基于数百台GPS接收机观测数据，分析了此次地震后中国上空电离层出现的扰动，首次揭示了瑞利波在远场导致的大尺度电离层扰动的空间结构.

1 数据预处理和扰动提取

2011年3月11日05∶46∶24 UT，在日本本州东海岸发生了MW9级的大地震，震中位于东经142.4°，北纬38.3°.本研究用到的数据来自中国境内约470台GPS接收机，地理位置的分布如图1所示，可见接收机距震中大圆距离大多为20°～40°，即约2200～4400 km之间.在地震发生前后，位于接收机网络上空的GPS卫星主要有G09、G15、G18、G21、G27等.接收机记录GPS卫星L1、L2两个频率上导航信号，利用信号的伪距和载波相位信息，解算了电离层垂直TEC.基于电离层薄壳假设，将TEC测量值的经纬度位置映射到400 km高度处.具体解算方法和过程参看Zhang等(2009, 2010).卫星的截止仰角设置为 30°，以减少多径效应等干扰因素.

图1 2011年3月11日日本地震震中和本文所用的GPS卫星、接收机的地理位置红星表示地震震中位置，虚线是距震中大圆距离分别为20°、30°、40°处，蓝点是接收机台站位置，品红点是地震发生时上空5颗GPS卫星投影到地面的位置.Fig.1 Map of the epicenter of 11 March 2011 Tohoku earthquake and locations of GPS satellites/receiversRed star is the location of epicenter, and dotted arcs denote the great circle distance of 20°, 30°, and 40° from the epicenter. Blue dots are locations of GPS receivers, and magenta dots are locations of 5 GPS satellites projected to the surface.

为了提取地震导致的TEC变化，对于每台接收机测量到的TEC时间序列进行滤波拟合以取得其背景值，此处采用滑动窗内多项式最小二乘拟合的Savitzky-Golay方法.用原观测值减去拟合得到的背景值序列，即得到扰动量的时间序列.由于不同接收机数据的时间分辨率不同，有的是30 s，有的是1 min，Savitzky-Golay滤波拟合的时间窗口都取为6 min，这样提取出的主要是周期短于6 min的短时扰动.图2a显示了对mnzi站数据进行处理的过程，叠加在原始TEC观测值上的实线是拟合得到的背景趋势，相减后可见在06∶20 UT左右发生短时的TEC扰动.

图2 提取TEC扰动的方法示例(a)和对部分台站数据进行处理的结果(b)(a)中上为原始TEC观测值和拟合得到的背景趋势(实线)，下为去掉背景后的扰动时间序列.图中数据均来自于G15卫星信号，并标示了台站代号(纬度，经度).Fig.2 Demonstration of the detection of TEC disturbances (a) and the results for several GPS stations (b)The upper part of panel (a) displays the original TEC data and fitted trend (solid line), and the lower part shows the time series of disturbance after removal of the trend. The TEC data are all from signals of G15 satellite, and the codes (latitude, longitude) of GPS stations are denoted.

图2b显示了部分接收机从G15卫星信号中得到的TEC扰动序列，可见扰动最大幅度可超过0.5 TECU(TEC单位，1 TECU=1016el/m2).图中由下而上，台站距震中距离越来越远，在距震中越远的台站，TEC扰动发生的时间越晚.

2 结果

在G09、G15、G18、G21、G27这几颗卫星中，G21卫星信号中测得的TEC扰动最少，幅度也小，这是由于对扰动的探测依赖于接收机-卫星连线的方向，具体详见后文讨论部分.此处的分析主要基于G09、G15、G18、G27四颗卫星的数据.

2.1 电离层扰动的起始时间

以G27卫星的数据为例，图3在时间与距震中的大圆距离坐标下组织TEC数据.从TEC扰动发生的时间—距离关系，可看到在距震中越远的地方扰动开始得越晚，说明这些扰动应与地震的发生和能量的向外传播有关.

图3 从G27卫星数据中提取的TEC扰动发生的时间-距离关系三条虚线分别表示地震波P波(品红)、S波(红色)、瑞利波(蓝色)随时间到达不同的震中距离(延时11 min).Fig.3 Time-distance diagram of the TEC disturbances from G27 satelliteThe three dashed lines denote the arrival of seismic P (magenta), S (red) and Rayleigh (blue) waves (delayed 11 min).

图3中自左向右，前两条斜线分别标出地震P波、S波到达不同距离处的时间作为参考，此处采用美国地质调查局给出的典型P波、S波传播时间和大圆距离的关系(参看https:∥www.usgs.gov/natural-hazards/earthquake-hazards/science/earthquake-travel-times)；第三条斜线给出瑞利波的传播，其速度依据该区域内多个地震仪记录拟合的结果，取为3.77 km·s-1(Hao et al., 2013).地震波能量从地表向上传输到电离层需要一定的时间，基于不同的观测手段和观测高度，以往的研究发现该时间为4～14 min不等(Tsai et al., 2011; Liu and Sun, 2011; Hao et al., 2012, 2013; Liu et al., 2016).本文所用GPS TEC反映的是电离层较高高度处的等离子体，通过比较瑞利波到达时间和TEC扰动开始的时间，发现该延时取为11 min时与观测符合较好.所以在图3中，将P波、S波和瑞利波的到达时间推迟了11 min显示(三条虚线)，以便与TEC扰动相对照.

图3显示，TEC扰动的开始时间与瑞利波的传播斜线基本符合，说明电离层扰动是随着瑞利波的传播同步发生的：瑞利波以3.77 km·s-1的速度从震中向外传播，在到达任一位置处都会激发大气波动，之后波动上传到电离层高度并导致TEC扰动.由图3可以看到，该现象的发生覆盖中国大部分区域，可一直影响到距震中大圆距离40°(超过4000 km)的中国西部地区.在更远距离的台站(参见图1，主要位于新疆)，则未观测到明显的扰动.

2.2 密集GPS网络观测到的电离层扰动空间结构

上一节给出了TEC扰动在距震中距离这一个维度上的时延、速度等特征，如果将数据显示在经纬度坐标下，则将揭示TEC扰动的二维空间结构信息.按图1所示，大部分GPS接收机位于中国中东部地区，但分布很不均匀，比如在江苏、山西两省接收机密集且均匀覆盖，但在其他区域也存在较大空白.图4首先给出江苏、山西两区域接收机的观测数据.

在经纬度坐标下，画出每一时刻G09、G15、G18、G27四颗卫星的数据，每30 s可得到一幅TEC扰动图，图4中对山西、江苏两区域的数据分别选择显示了四个代表性时刻.由于这四颗卫星相对于接收机分别处于不同的方位角(参看图1中的卫星位置)，TEC数据映射到电离层的位置会有所错开，从而覆盖了更大的空间区域.以图4e为例，山西地区接收机的TEC数据点分为了两群，左边一群为接收G18卫星信号所得到的，而右边一群则是G15、G27、G09卫星的数据重叠在一起.

图4 江苏和山西两区域接收机观测到的TEC扰动经纬度分布(a—d) 江苏； (e—h) 山西.观测数据来自G09、G15、G18、G27四颗卫星，虚线为延时11 min显示的地震波S波(红色)、瑞利波(蓝色)的位置.Fig.4 Geographic distribution of TEC disturbances from receivers in JiangSu and ShanXi(a—d) JiangSu; (e—h) ShanXi. The data are from signals of G09, G15, G18, and G27 satellites. The dashed curves denote arrival locations of seismic S (red), Rayleigh (blue) waves (11 min delayed).

按2.1节所述，简化地认为地震波P波、S波、瑞利波从震中向各个方向均匀地传播，并假设地震波能量从地面传播到电离层需要11 min.与图3类似的，将地震波的到达的位置延时11 min后用虚线标示在图4中.在图4中可见，地震波到达中国大陆东海岸后，该地上空最早在06∶07UT出现TEC扰动(图4b);在06∶10UT，TEC扰动呈现条纹结构，走向大致为西北-东南(图4c)；由06∶10UT和06∶12UT图可见扰动继续向南向西运动，条纹走向渐趋于东西方向，并逐渐减弱(图4c—d).山西地区的接收机也观测到类似扰动，如图4(e—h)所示，西北-东南的扰动条纹在瑞利波到达后出现，并同样向南向西运动.

2.3 电离层扰动的大尺度空间结构

图5给出所有台站观测的TEC扰动在经纬度坐标下的二维分布，与图4类似的，汇集G09、G15、G18、G27四颗卫星的数据后，提高了空间覆盖，中国中东部地区的台站空白区域大多得到了填补，从而可给出完整的电离层扰动的大尺度空间结构.

图5 大尺度TEC扰动的经纬度分布三条虚线分别表示地震波P波(品红)、S波(红色)、瑞利波(蓝色)的位置(延时11 min).震中和卫星位置的标识与图1相同.Fig.5 Geographic distribution of large-scale TEC disturbancesThe dashed curves denote arrival locations of seismic P (magenta), S (red), Rayleigh (blue) waves (11 min delayed). Locations of the epicenter and GPS satellites are marked in the same way as in Fig.1.

在图5中把地震波到达与电离层扰动的位置相对照，可以看到：在06∶05UT，地震波尚未到达观测区域，没有扰动发生；在06∶09 UT，在瑞利波(蓝色虚线)到达中国大陆后，电离层中的扰动最早从东部沿海地区开始出现；随着瑞利波的传播，扰动区域逐渐扩展到更西和更南的区域，并形成多个条纹状的扰动结构.条纹的走向在中纬度主要为西北-东南走向，而到了低于北纬30°的区域，条纹大致为东西走向.

3 讨论和结论

与一个固定的点源激发大气波动不同，地震波是在不断传播的，而成为运动的激发源.在该过程中激发的一系列大气波动会叠加起来，形成沿瑞利波传播方向倾斜的波前.由于瑞利波的波速(3～4 km·s-1)比大气声速(340 m·s-1)高很多，该波前与地面的夹角约为5～10°，在电离层高度声速约800 m·s-1，则夹角增大至10～20°(Rolland et al., 2011b; Kakinami et al., 2013; Inchin et al., 2020).接收机-卫星连线方向与波前方向越平行，则越能敏感地捕捉到扰动带来的等离子体变化.反之，如果两个方向接近垂直，接收机-卫星连线将同时穿过扰动的波峰波谷，因为TEC测量的是信号路径上的积分电子含量，所以正负扰动抵消，此时会无法探测到明显的扰动现象.G21卫星与接收机连线的方向就属于这种情况，所以如前文所述G21卫星数据中只探测到少而弱的TEC扰动.

瑞利波从震中向四周传播，理应在各个方向上都形成上述的波前结构，相应的电离层扰动是以震中为圆心的环状结构，模拟研究也再现了这种环形TEC扰动的主要特征(Rolland et al., 2011b; Inchin et al., 2020).对于此次地震，在震中附近的GPS网络确实观测到圆环状的、从震中向四周传播的TEC扰动，速度与瑞利波速相近(Tsugawa et al., 2011; Kakinami et al., 2013).但是，在远场类似的观测证据还比较少，如在台湾地区的GPS网观测到大致东西走向的TEC扰动(Chen et al., 2011).本文图4则显示，江苏地区扰动先开始，条带起初是偏西北—东南方向，后来趋于大致东西方向；山西地区的扰动主要为西北—东南走向.结合这两个区域的观测和相对于震中的方位可以看出，扰动开始时的形态似乎与瑞利波形成的大气波动的叠加波前有关，即以震中为圆心的环状结构.但图5显示，在更大的空间尺度上，TEC扰动结构多以东西向为主，并非以震中为圆心的圆环结构，这与前述在震中附近观测到的现象完全不同.

图4和图5中扰动条带的运动方向主要向南，这是由于地磁场的倾角特点会导致扰动的南北不对称，即大气扰动会在朝向赤道一侧(在北半球即为向南)产生较强的电离层扰动(Rolland et al., 2011b; Kakinami et al., 2013).但需要强调的是，在南方(低纬地区)出现的扰动并不是从北方传播过去的，而是由瑞利波在当地激发产生为主，原因如图3所示，TEC扰动的出现是与瑞利波的传播同步的，它们具有相同的约3.77 km·s-1的水平速度；并且从图5显示的时间切片可知，条带结构传播速度较慢，不足以迅速地从北方传播到低纬地区.这不同于Chen等(2011)的研究，他们认为台湾上空观测到的TEC扰动直接来源于震中产生的次声波，并估算其水平传播速度为1000～1700 m·s-1.

图5显示条带在中纬度产生时更偏向西北—东南走向，之后逐渐转向，特别是在低纬度区域，大致为东西走向.这说明扰动波前在传播时发生了逆时针的旋转，该过程可能与地磁场作用下的中性-离子耦合有关.在对大气波动导致的电离层行扰的研究中发现：因为在电离层F区离子磁回旋频率远高于离子-分子碰撞频率，中性分子推动等离子体沿磁力线运动的效率是最高的，同时在磁场方向上离子曳力对大气波动的阻尼作用最弱，这最终导致大气波动和电离层行扰的传播都以沿磁场方向为主(Liu and Yeh, 1969; Kelley and Miller, 1997).本文观测到的TEC扰动结构的传播和转向，可能也经历了类似的过程.

综合本文的观测事实可以推断：在距震中不同距离处，瑞利波到达后其能量从地表经大气耦合到达电离层高度，并导致TEC扰动.但观测到的电离层大尺度扰动并不是以震中为圆心的圆环结构，说明它不是第一时间由瑞利波产生的，而是扰动产生后在传播过程中被改变形成的形态：瑞利波激发的电离层扰动在电离层中自由传播，其方向在地磁场作用下转为南向为主，形成东西走向的TEC条带.本文首次给出了地震远场瑞利波产生的电离层大尺度扰动结构，但这些扰动激发和传播的具体过程还有待更深入的研究.

致谢感谢子午工程提供部分GPS接收机观测数据.