2015, Vol. 58
2. 武汉大学电子信息学院, 武汉 430072
2. Department of Space Physics, School of Electronic Information, Wuhan University, Wuhan 430072, China
VLF electric and magnetic field spectral data from DEMETER satellite is used to analyze the diurnal, seasonal, and annual variations of ionospheric electromagnetic responses from 2005 to 2009, which were excited by the NWC communication station located in Australia. The Maxwell equation in the ionosphere is solved using the transfer matrix method to simulate the propagation characteristics of VLF electromagnetic waves excited by NWC transmitter. Compared with the observation result and simulation result, the electromagnetic response is obtained with different ionospheric backgrounds.
The result shows the electric field at daytime is almost 40 dB larger than at nighttime, and the magnetic field is 15 dB. According to the data analysis result, it is found that seasonal and annual variations of the electromagnetic field is not very significant and solar activity has some effect on the response, which is smaller when the solar is more active. The simulation results are consistent with the data analysis. The variation of the response may be caused by the distribution of electron density below 250 km in the ionosphere. If the electron density is higher, the electromagnetic response will be lower.
The variation of the electric field is much more obvious than the magnetic field, which means finding earthquake anomalies from electric data of satellite may be more efficient than the magnetic field. Getting precise electron profiles below 250 km will be very helpful for studying the electromagnetic anomalies in the ionosphere related with earthquakes.
电离层等离子体主要由太阳极紫外(EUV)和X射线辐射造成上层大气的中性成分电离产生.因此，电离层具有明显的日变化，27天变化和太阳活动周期变化(熊年禄等，1999；Rich et al.，2003; Liu et al.，2006，2007)等.这些变化构成了电离层背景特征.一些研究者利用法国DEMETER卫星观测到了发生在地震活动区上空的电磁场扰动(Parrot et al.，2006; Zhang et al.，2009，2014)，它们有很宽的频带范围(Parrot，1994; Molchanov and Hayakawa，2008)，从Ultra Low Frequency(ULF，DC~15 Hz)到Extreme Low Frequency(ELF，15~30 Hz)甚至到几十千赫兹的Very Low Frequency(VLF，3~30 kHz)频段.地震前兆激发的电磁辐射在不同电离层背景下呈现怎样的电磁响应变化特征？了解这些特征可为我们判别地震电离层电磁前兆提供依据.由于地震电磁辐射的辐射能量及频率等信息目前尚无法确定，因此我们可以借助全球范围内的地面VLF潜艇通信台和通信台发射的电磁信号，它们的发射功率和发射频率都是已知的(表 1)，其信号进入电离层能够被卫星观测仪器记录.本文将利用澳大利亚NWC发射的大功率、高信噪比VLF信号，研究其在电离层中激发响应的时变特征，为地震导致的电离层电磁辐射异常判定提供背景参考.在众多发射站中NWC是目前南半球最大的VLF发射源，并且由于它处于中纬地区，也是对辐射带影响最大的人工发射源(Kulkarni et al.，2008; Li et al.，2014).
|
表 1 主要VLF发射站列表 Table 1 List of major VLF transmitter |
本文通过统计分析DEMETER卫星记录的NWC通信台激发的电磁场数据，研究DEMETER卫星运行期间(2005-2009年)的电磁响应的日变化、季节变化和年变化特征.采用传递矩阵解法(Nagano et al.，1975; Nygrén，1982; Yagitani et al.，1994)求解电离层中的电磁波传播方程，模拟VLF电磁波在不同电离层环境中传播的变化特征，与观测记录进行对比，给出不同时间的电离层背景对电磁波传播的影响，为寻找地震导致的电离层电磁辐射异常判定提供背景参考.
2 电磁响应时变特征统计分析法国DEMETER(Detection of Electro-Magnetic Emission Transmitted from Earthquake Regions)卫星自2004年8月运行至2010年12月，为准太阳同步极轨卫星，轨道倾角98.23°，运行高度710 km(2005 年12 月中旬改为660 km)，卫星重量130 kg. 卫星上的电磁场探测载荷有： 1)电场探测设备ICE，用四个电传感器测量电场的三分量(利用各传感器之间的电势差)，频率范围从DC/ULF(0~15 Hz)、ELF(15 Hz~1 kHz)、VLF(15 Hz~17.4 kHz)到HF(15 Hz~3.175 MHz)；2)磁场探测设备IMSC，由三个正交的磁探测线圈组成，频率范围为DC~ VLF.DEMETER卫星观测有两种模式：巡查(survey)模式和加密(burst)模式，其中巡查模式数据量小，用于记录不变纬度低于65°的所有地球区域的数据；加密模式用在地震多发区域，数据量大，具有更高的时间和频率分辨率.针对电磁场探测数据，详查模式(burst)保存地震多发区域的电场和磁场波形数据，采样率40 kHz；巡查模式(survey)保存全球的电场和磁场谱数据，谱分辨率约20 Hz，时间分辨率约为2 s.本文研究利用的是survey 模式下的DEMETER电场探测仪(ICE)和感应式磁力仪(IMSC)记录的VLF频段电场和磁场谱数据，数据从DEMETER卫星中心获得.
澳大利亚VLF通信台NWC，位于地理坐标21.816°S，114.166°E处，发射功率为1000 kW，发 射中心频率为19.8 kHz，发射带宽约为200 Hz.NWC 通信台为连续工作，在2007年6月关机，直到2008年1月底重新开机，此次维护时间约为8个月左右.
我们从DEMETER卫星观测的电磁场功率谱数据中提取以19.8 kHz为中心频率、200 Hz带宽内的功率谱数据，求得对应的均方根(RMS)电场和磁场强度，认为是NWC发射站信号在电离层中的响应，对获得的电场和磁场强度数据进行统计分析.图 1是一天中DEMETER卫星处在上升轨道时(地方时夜间)测量得到电场强度的全球分布图，其中红色五角星表示发射站位置.DEMETER卫星通过NWC发射站区域的时间为地方时10时和22时左右.数据以半轨道记录，存储在星载高容量存储器上，当卫星经过地面控制站时，每天和地面发生两次数据交换.
 |
图 1 DEMETER卫星一天所有上升轨道观测电场数据(19.8 k±200 Hz)投影图 Fig. 1 Projection of the electric field data(19.8 k±200 Hz)of all the up orbits of DEMETER satellite during one day |
Cohen和Inan(2012)利用DEMETER卫星(660 km高度)整个运行期间观测到的NWC发射站信号的电场数据，研究了电磁响应的空间分布特征，发现空间辐射图样呈现同心圆分布，且辐射的同心圆中心点在地面的投影相对于发射源所在的位置有沿着经向的偏移，位于发射源上空80 km处的磁力线在660 km高度的穿刺点位置，这是由于在80 km以上VLF电磁波主要沿着磁力线传播导致.因此本研究在统计分析发射源在电离层中激发的响应时，没有选取以发射源为中心的区域，而是选取的以发射源在660 km高度激发的辐射图样的中心为数据中心.
因此，本文选取同心圆的中心位置周边10°×10°地理范围内，2005年至2009年的DEMETER卫星电磁场数据进行统计分析，研究地面VLF辐射源在电离层中激发电磁响应的年变化、季节变化以及日变化规律.首先比较了2005年电场和磁场春季(8-10月)、夏季(11月-次年1月)、秋季(2-4 月)、冬季(5-7月)四个季节的季节变化，以及白天(10LT)与夜间(22LT)之间的日变化.
图 2为2005年NWC发射站在电离层中激发的电场的响应的日变化和季节变化，从图中可以发现电场强度有显著的日夜变化，而季节变化不显著.场强日变化可达40 dB，季节变化最大约7 dB.
 |
图 2 2005年太阳活动低年VLF发射源激发的电场强度不同时间统计分布(第一列为四个季节10LT，第二列为四个季节22LT) Fig. 2 Statistical distribution of electric field intensity excited by VLF transmitter in 2005 low solar activity(First column shows the results at 10LT and second column at 22LT in the four different seasons) |
图 2还可见在10LT，0~3 dB-μV/m(即1~1.4 μV/m)处概率值很大，这是因为在10LT，由于 电子密度很大，NWC信号除了在发射源上空附近 能收到一定强度的信号外，选取的研究区域的其他位置基本收不到信号，也就是说辐射源信号到达卫星高度后在水平方向的衰减很快.
图 3为2005年NWC发射站在电离层中激发的磁场的响应的日变化和季节变化，从图中可以发现磁场强度比电场强度变化小，但仍有显著的日夜变化，而季节变化不显著.10LT与22LT之间磁场强度变化可达20 dB，而四个季节的磁场强度变化约为3 dB.
 |
图 3 2005年太阳活动低年VLF发射源激发的磁场强度不同时间统计分布(第一列为四个季节10LT，第二列为四个季节22LT) Fig. 3 tatistical distribution of magnetic field intensity excited by VLF transmitter in 2005 low solar activity(First column shows the results at 10LT and second column at 22LT in the four different seasons) |
为了研究电磁响应的年变化特征，从图 2可知春季夜间电磁场响应最强，因此本研究利用2005年至2009年DEMETER卫星数据分析了春季夜间电磁场的年变化，结果如图 4所示，左列五幅图分别是 2005年至2009年电场统计结果，右列五幅为磁场统计结果.可以发现电场响应的年变化不明显，磁场有同样结果.2005年至2009年的电磁场强度年变化并不显著的可能原因是2005年至2009年太阳活动均处于较低年份，电离层电子密度变化不显著.电离层电子密度对电磁响应的影响将在本文第4节数值模拟中进一步分析.
 |
图 4 2005年至2009年春季夜间VLF发射源激发的电场强度统计分析 Fig. 4 Statistical distribution of electric field intensity excited by VLF transmitter from 2005 to 2009 |
值得注意的是2007年8月至10月的电场和磁场强度明显小于其他年份，这是由于2007年6月至2008年1月曾暂停发射8个月.
3 计算方法由于VLF频段波长很长，电离层参数在一个波长变化范围内变化剧烈，因此把电离层假设为水平分层介质.而传播矩阵法是计算低频电磁波在分层介质中传播的一种手段.VLF电磁波穿透低电离层，由于碰撞作用被强烈吸收，电子与离子及中性粒子的碰撞不可忽略.穿过电离层D区后地磁场对电 磁波传播起主导作用，因此模型计算时必须考虑地磁场作用导致的电离层各向异性.综上所述，实质是解决水平分层各向异性有耗介质的电磁波传播问题，该问题一直以来是电磁波传播计算中的难点.
计算选取的笛卡儿坐标系如图 5所示，y轴和z轴分别指向地磁北和垂直方向.x轴与y和z轴正 交.一个平面波入射只有z方向不均匀的各向异性电离层，入射面为x-z面，入射角为θi，斜向地磁场 H 0的方向余弦为(l，m，n).
 |
图 5 计算使用的笛卡儿坐标系，入射面为x-z平面，地磁场H0的方向余弦为(l，m，n) Fig. 5 Geometry of the problem. The plane of incident is the x-z plane. The direction of cosines of the Earth′s magnetic field H0 is (l，m，n) |
设场的时谐因子为 ejωt，则介质中电磁波满足下面的Maxwell方程组
其中 ε0，μ0，I分别表示自由空间的介电常数、磁导率和单位矩阵.对于VLF频段只考虑电子的影响，电离层的电极化率矩阵(Yeh and Liu，1972)为
其中 U=1+i
，ν 主要指电子与中性粒子、重离子的碰撞.
称为电子的等离子体角频率，这里 m、e和N 分别是电子的质量、电荷和体积密度.
称为电子的磁旋角频率.
假设电离层如图 6水平成层，在x和y方向均匀，只有z方向变化.则每层的入射波场各分量的时空因子为 exp[jωt－jk0(Sx+Cz)]，其中k0 是自由空间中的波数，S=sinθi，C=cosθi.
 |
图 6 水平分层介质中的波 Fig. 6 Waves in the horizontal stratified media |
利用Snell法则消去方程(1)中Ez和Hz分量，Maxwell方程化为更简洁的形式(Budden，1961)
解方程(3)时，首先估计第m层的T矩阵，得矩阵两特征向量对应上行波，作为初始条件向下利用传递矩阵法(Chew，1990)迭代至第1层得全波解.
其中Ks 称为各层的传播矩阵；A和B为任意常数，由 z=z1 处的入射波场决定.
当介质有耗时，(4)式的传递矩阵连乘时会出现溢出问题，本文采用Nagano提出的格兰施密特正交化方法(Kunz，1957)克服该数值溢出问题.
4 数值模拟结果计算时模型设80 km以下为自由空间；80~700 km为水平分层各向异性电离层，分层厚度为1.55 km；700 km以上认为是均匀各向异性半空间. 利用International Reference Ionosphere model(IRI模型)(Bilitza，2001)计算电子数密度在80~700 km的剖面.电子碰撞频率在80~300 km采用Cummer(2000)提出的模型；300~700 km库仑碰撞起主导作用，采用Helliwell(1965)的碰撞模型.利用International Geomagnetic Reference Field(IGRF地磁场模型)(Barton，1997)计算NWC发射站处地磁场强度45 μT，倾角为-55°，偏角设为0.1°.
4.1 NWC电磁响应日变化季节变化规律为解释NWC在电离层激发的电磁响应呈现的变化特征，本研究选取了NWC发射源2005年春分、秋分、夏至、冬至10LT及22LT的电子密度剖面(图 7)进行四个季节的模拟.
 |
图 7 计算使用的2005年四个季节电子密度剖面(实线表示10LT，虚线表示22LT) Fig. 7 Electron density profile used in the calculation of four seasons in 2005(Solid line present 10LT，and dashed line present 22LT) |
模拟结果显示电场强度日变化最大达30 dB以上，而季节变化最大约10 dB，夜间夏季和冬季的差别最大，夏秋几乎无差别，白天秋季与其他三个季节 差别均较大，磁场强度四个季节之间的夜间几乎无变化，白天最大差别约15 dB，这与DEMETER卫星的数据分析结果一致.
分析IRI模型的电子密度分布发现：250 km以下季节区别不大，而250 km以上季节差别较大.图 8为不同电子密度分布情况下，电磁场随传播高度的变化.计算结果表明，电磁波能量的传播衰减随季节变化不大，主要体现在夜间与白天的显著变化.这一模拟计算结果与前文DEMETER卫星数据分析结果一致，并且与250 km以下电子密度的时间变化特征一致，说明250 km以下电子密度分布对电磁波能量起着重要作用，电子密度越大，电磁波能量衰减越大.
 |
图 8 2005年电磁场随高度的变化 Fig. 8 The electric and magnetic field vary with altitude in 2005 |
数据分析中2005年至2009年年变化特征不显著的原因可能是2005年至2009年太阳活动都处于较低年份，因此为分析太阳活动对电磁响应的影响，我们对DEMETER卫星运行期间太阳活动最强年份2005年及太阳活动最低年份2009年冬夏两季夜间的电磁场也进行了模拟，结果发现，如图 9所示，2005年和2009年的电磁场场强变化不大.2009年的冬季电场及磁场均比2005年强，而夏季则比2005年弱，这与数据分析的结果也一致.
 |
图 9 太阳活动不同年份2005年(实线)和2009年(虚线)电磁场强度对比 Fig. 9 Comparison of electromagnetic field between different solar activities in 2005(solid line) and 2009(dashed line) |
从NOAA网站获取太阳射电通量F10.7指数 数据，给出了年均值变化趋势图，以及2001年和2005年至2010年F10.7指数的年变趋势，如图 10所示，可以发现，2005年与2009年实际上都处于太阳活动较低的阶段，为判断太阳活动对卫星高度电磁响应的影响，我们选择太阳活动最强的2001年作为对比，选取2001 年夏至冬至22LT的IRI电子浓度剖面，以NWC源作为输入进行数值模拟，与2009年电场和磁场结果对比，发现在太阳活动较强的2001年，电场及磁场的响应较2009年显著降低(图 11).因此，太阳活动对电离层电磁响应强度有 影响，并且太阳活动越强烈，电离层中的电磁响应越小.
 |
图 10 10.7cm 太阳射电通量年均值变化趋势及F10.7年变化趋势 Fig. 10 Annual change of average solar radio flux F10.7 |
 |
图 11 太阳活动不同年份2001年(实线)和2009年(虚线)电磁场强度对比 Fig. 11 Comparison of electromagnetic field between different solar activities in 2001(solid line) and 2009(dashed line) |
本文模拟了NWC发射源的电离层电磁响应时间变化特征，模拟结果与数据分析的结果一致，利用本研究建立的模型进行数值模拟，将有可能帮助判定地震相关的电离层异常.
图 8中2005年秋季(20050322)的仿真结果电场磁场均远远小于其他三个季节，这可能与我们选 取的IRI模型20050322这一天的电子浓度在200~250 km 明显大于其他三个季节电子浓度(图 7)有关，这与本文前面推断的250 km以下电子密度分布对电磁波能量起着重要作用的结论相吻合.因此，在判定地震电离层电磁异常的时候应更准确了解电离层电子密度的变化.还需指出的是本研究的模型是从电离层底部(50 km)到卫星高度，在以后的工作中还需考虑从地表到电离层底部的传播模型.
5.2 结论本文统计了DEMETER卫星在不同背景电离层下对NWC发射源的电磁响应，并进行了仿真模拟.得出以下结论：
(1)电场的日变化达40 dB和季节变化约7 dB，磁场日变化约15 dB，季节变化不显著，仿真模拟结果与数据分析结果一致.电场的变化比磁场的变化更明显，因此通过电场数据寻找地震电离层电磁异常可能更有效.
(2)250 km以下电子密度分布日变化显著，季节变化不大，而250 km以上电子密度季节变化显著.计算结果与DEMETER卫星数据分析结果均显示电场和磁场的日变化显著而季节变化不显著，这与250 km以下电子密度的时间变化特征一致，说明250 km以下电子密度分布对电磁波能量衰减起着更重要作用，电子密度越大，电磁波能量衰减越大.
(3)在寻找地震电离层电磁异常时，为了得到异常特征应该着重研究 250 km以下的电离层电子密度变化.
(4)太阳活动变化对VLF发射源的电离层电磁响应有一定的影响，太阳活动越强烈，卫星高度的电磁响应越低.
致谢 感谢法国DEMETER卫星任务中心提供的电磁场数据；日本金泽大学的Nagano教授在编程中提供了指导和帮助，作者在此表示感谢.感谢评审人的建议以及编辑的辛勤劳动.
| [1] | Barton C E. 1997. International Geomagnetic Reference Field: the seventh generation. J. Geomagn. Geoelec., 49(2): 123-148. |
| [2] | Bilitza D. 2001. International reference ionosphere 2000. Radio Sci., 36(2), doi: 10.1029/2000RS002432. |
| [3] | Budden K G. 1961. Radio Waves in the Ionosphere. London and New York: Cambridge Univ. Press, 389-390. |
| [4] | Chew W C. 1990. Waves and Fields in Inhomogeneous Media. New York: Van Nostrand Reinhold, 109-114. |
| [5] | Cohen M B, Inan U S. 2012. Terrestrial VLF transmitter injection into the magnetosphere. J. Geophys. Res., 117: A08310, doi: 10.1029/2012JA017992. |
| [6] | Cummer S A. 2000. Modeling electromagnetic propagation in the Earth-ionosphere waveguide. IEEE Transactions on Antennas and Propagation, 48(9): 1420-1429, doi: 10.1109/8.898776. |
| [7] | Helliwell R A. 1965. Whistlers and Related Ionospheric Phenomena. California: Stanford University Press, 63-76. |
| [8] | Kulkarni P, Inan U S, Bell T F, et al. 2008. Precipitation signatures of ground-based VLF transmitters. J. Geophys. Res., 113: A07214, doi: 10.1029/2007JA012569. |
| [9] | Kunz K S. 1957. Numerical Analysis. New York: McGraw-Hill, 236-238. |
| [10] | Li J D, Spasojevic M, Harid V, et al. 2014. Analysis of Magnetospheric ELF/VLF wave amplification from the Siple transmitter experiment. J. Geophys. Res. Space Physics, 119(3): 1837-1850, doi: 10.1002/2013JA019513. |
| [11] | Liu H X, Stolle C, Förster M, et al. 2007. Solar activity dependence of the electron density in the equatorial anomaly regions observed by CHAMP. J. Geophys. Res., 112(A11): A11311, doi: 10.1029/2007JA012616. |
| [12] | Liu L B, Wan W X, Ning B Q, et al. 2006. Solar activity variations of the ionospheric peak electron density. J. Geophys. Res., 111: A08304, doi: 10.1029/2006JA011598. |
| [13] | Molchanov O A, Hayakawa M. 2008. Seismo-Electromagnetics and Related Phenomena: History and Latest Results. Tokyo: Terra Sci. Pub. Co. |
| [14] | Nagano I, Mambo M, Hutatsuishi G. 1975. Numerical calculation of electromagnetic waves in an anisotropic multilayered medium. Radio Sci., 10(6): 611-617, doi: 10.1029/RS010i006p00611. |
| [15] | Nygrén T. 1982. A method of full wave analysis with improved stability. Planet. Space Sci., 30(4): 427-430, doi: 10.1016/0032-0633(82)90048-4. |
| [16] | Parrot M. 1994. Statistical study of ELF/VLF emissions recorded by a low-altitude satellite during seismic events. J. Geophys. Res., 99(A12): 23339-23347, doi: 10.1029/94JA02072. |
| [17] | Parrot M, Berthelier J J, Lebreton J P, et al. 2006. Examples of unusual ionospheric observations made by the DEMETER satellite over seismic regions. Physics and Chemistry of the Earth, Parts A/B/C, 31(4-9): 486-495. |
| [18] | Rich F J, Sultan P J, Burke W J. 2003. The 27-day variations of plasma densities and temperatures in the topside ionosphere. J. Geophys. Res., 108: 1297, doi: 10.1029/2002JA009731.A7. |
| [19] | Xiong N L, Tang C C, Li X J. 1999. Introduction of Ionospheric Physics (in Chinese). Wuhan: Wuhan University Press. |
| [20] | Yagitani S, Nagano I, Miyamura K, et al. 1994. Full wave calculation of ELF/VLF propagation from a dipole source located in the lower ionosphere. Radio Sci., 29(1): 39-54,doi: 10.1029/93RS01728. |
| [21] | Yeh K C, Liu C H. 1972. Theory of Ionospheric Waves. New York: Academic Press. |
| [22] | Zhang X, Shen X, Zhao S, et al. 2014. The characteristics of quasistatic electric field perturbations observed by DEMETER satellite before large earthquakes. Journal of Asian Earth Sciences, 79: 42-52. |
| [23] | Zhang X M, Qian J D, Ouyang X Y, et al. 2009. Ionospheric electromagnetic perturbations observed on DEMETER satellite before Chile M7.9 earthquake. Earthquake Science, 22(3): 251-255. |
| [24] | 熊年禄, 唐存琛, 李行健. 1999. 电离层物理概论. 武汉: 武汉大学出版社. |