地磁感应地电场时频变换快速计算方法

倪政泽， 王泽忠， 李宇妍， 司远

(华北电力大学高电压与电磁兼容北京市重点实验室， 北京 102206)

地磁暴是指由于太阳活动而造成地球磁场全球性剧烈扰动的现象。地球磁场的剧烈扰动将在地面中感应出电场，该感应电场会在由输电线路，中性点接地的变压器及大地构成的回路中产生地磁感应电流(geomagnetically induced current，GIC)，GIC频率在0.01～0.000 1 Hz范围内变化，这种准直流电流对变压器等电网设备会造成严重威胁[1-4]。同时地磁场变化迅速，如何准确快速地求取磁暴发生时的地面感应地电场是防治磁暴灾害的关键。

传统的地电场计算方法需要根据不同地区土壤的电性结构确定大地电导率模型，文献[5]提出建立均匀大地模型，没有充分考虑地电场的实际分布情况。文献[6]提出建立水平分层大地电导率模型，利用平面波法[7]计算地电场，对于中国特高压长距离输电线路，其空间尺度达到几百公里甚至上千公里，地质结构复杂[8]，各地区土壤结构差异较大，难以建立准确模型。文献[9-10]建立了包含海陆分界面的二维大地电导率模型，但没有电导率的实测数据，文献[11-12]提出建立三维大地分层分区电导率模型，利用地磁台站记录的地磁场变化量作为地面边界条件，采用有限元法计算地电场，但自适应剖分网格过多，计算时间过长。无法准确得出各地大地模型并快速计算成为传统地电场方法的局限性。文献[13]提出时谐拟合算法将时变的地磁场转化为频域相量法计算，数据拟合得到的仅为单一频率，与实际相差较大。文献[14]提出利用视电阻率数据计算磁暴感应地电场频域幅值，无法准确选取包含磁暴最大幅值的脉冲，计算结果准确性受数据影响较大。

地磁测深是利用天然交变电磁场研究地球电性结构的一种物理勘探方法，具有成本低，工作方便，不受高阻层屏蔽的优点。在地下岩石电性分布不均匀(有两种或两种以上导电性不同的岩石或矿石)或地表起伏不平的情况下，仍按测定均匀水平大地电阻率的方法和计算公式求得的电阻率称之为视电阻率，其综合反映了大地的电性结构，是地磁测深法的重要数据。现利用2004年11月7日地磁台实测秒级数据与视电阻率数据在频域计算得到频域电场值，再利用傅里叶反变换得到时域电场值。利用这种方法对地磁暴发生时刻阿尔山地区，阿拉善盟，大别山地区，广东东莞谢岗及增城朱村的感应地电场最大值进行对比。该计算方法简化了地面感应地电场的计算，为工程应用提供了有效算法，在地磁暴灾害防控中有着重要意义。

1 视电阻率在地电场算法中的应用

在大地电磁测深方法中，考虑到应用的观测频率范围以及构成地壳浅部介质的电阻率范围，在大地介质中可忽略位移电流对场分布的影响[15]，于是谐变场的Maxwell方程组表示为

(1)

(2)

(3)

(4)

式中：H为地磁场强度向量；E为地电场强度向量；ω为角频率；μ为磁导率；γ为电导率；i为虚数单位。

如图1所示建立笛卡尔坐标系。地磁暴对地电场的影响可看作是平面电磁波垂直入射大地介质(即OZ方向)中产生的电磁场，同时对于电磁测深方法，视电阻率ρ的测定以均匀水平大地为前提。

O为观测点(坐标原点)；OXY为水平面；OX指向地磁北；OY指向地磁东；OZ为铅直向下；N为地磁北方向；E为地磁东方向图1 地磁场坐标系的定义Fig.1 Definition of geomagnetic coordinate system

根据均匀介质中的平面波理论[16]，有

(5)

(6)

式中：Zxy为TE模式下的波阻抗；Zyx为TM模式下的波阻抗；Ex和Ey分别为地电场南北分量和东西分量；Hx和Hy分别为地磁场南北分量和东西分量；ρxy为TE模式下的视电阻率；ρyx为TM模式下的视电阻率；i为虚数单位；μ为磁导率，大地磁导率近似为μ0。

若只考虑视电阻率与对应波阻抗的幅值，有以下关系：

(7)

(8)

目前已有通过地磁测深得到的视电阻率数据，可通过式(7)和式(8)得到对应阻抗Zxy和Zyx。将时域地磁场数据经快速傅里叶变换(FFT)后得到频域磁场数据，利用式(5)和式(6)的关系计算地电场Ex和Ey的频域数值，然后进行快速傅里叶反变换(IFFT)即可得到地电场Ex、Ey的时域数值。该计算方法以视电阻率数据作为支撑，将时域中多点求解问题转化为频域中的单点求解问题。

1.1 磁场数据的来源

目前收集到地磁台磁暴数据有分钟级和秒级地磁场三分量(即地磁场总强度F、磁偏角D和地磁场水平分量H)数据，为了得到地磁南北和东西分量变化率脉冲幅值，现使用秒级数据。所收集到的数据中存在一些因设备故障而导致的错误数据，剔除错误数据后，能够使用的据包括2004年11月7—9日黑龙江德都(DED)、广东广州(GZH)、甘肃嘉峪关(JYG)、海南琼中(QGZ)、湖北武汉(WHN)、新疆乌鲁木齐(WMQ)和新疆喀什(KSH)7个地磁台的测量数据。各台站的地理坐标如表1所示。

表1 地磁台站地理坐标Table 1 Geographic coordinates of geomagnetic stations

1.2 视电阻率数据的选取与处理

选取的阿尔山火山区大地电磁测深数据由文献[17]中获得，观测了一条西北向测线上的7个大地电磁测深点，选取其中编号3、4两个测点；阿拉善盟的数据由文献[18]中获得，其观测了东西向测线共65个测点，选取其中编号43、58两个测点；大别山的数据由文献[19]中获得，其观测了自西向东3条大地电磁测深侧线，测线LA沿NE26°方向，测线LB沿着NE28°方向，选取测点LA7、测点LB7两个测点；广东增城朱村，东莞谢岗的数据由广东省地震局与国家地震局地质研究所开展的大地电磁测深研究中获得，其电阻率幅值曲线如图2所示。

1.3 磁暴期间磁场数据的选取与处理

为了提高计算的准确度和精度，根据选取同一纬度范围内，距离最近的地磁台站作为地磁场数据来源的准则，选取德都地磁台(DED)，嘉峪关地磁台(JYG)，武汉地磁台(WHN)，广州地磁台(GZH)的秒级磁场数据分别作为阿尔山地区，阿拉善盟，大别山造山带，广东增城朱村与东莞谢岗的磁场数据来源。

在2004年11月7日18：28：08的磁暴急始时刻，监测广东岭澳核电站内的中性点电流峰值达到了-47.2 A(负号表示电流从大地流入中性点)[20]，故选取2004年11月7日17：00—21：00各地磁台磁场数据作为原始数据进行分析，各地磁台该时段的磁场数据如图3所示，纵轴单位为nT，1T=1×109nT。

图2 测点视电阻率幅值曲线Fig.2 Apparent resistivity amplitude curve of measuring point

地面感应电场与地磁场存在半阶导数关系[21]，即

(9)

(10)

式中：Ex(ω)、Ey(ω)为感应地电场频域的x、y分量；Bx(ω)、By(ω)地磁场频域的x、y分量；j为虚数。其对应的时域内地电场分量和对应地磁场分量的关系式：

(11)

(12)

式中：u为变量名；t为待求电场的时间；gx(u)=dBx/dt、gy(u)=dBy/dt为地磁场时域的x、y分量变化率；μ0为真空磁导率；γ为大地电导率。将式(11)和式(12)归一离散化后，得

(13)

式(13)中：bn=Bn-Bn-1为地磁场分量的一阶差分；M为一足够大的数；TN为待求电场的第N个时刻；Tn为第n个采样点的时间。由式(13)即可由磁场数据的时间序列计算出地面电场相对值的时间序列，2004年11月7日17：00—21：00各地台感应地电场相对值如图4所示。

图3 各地磁台磁场幅值Fig.3 Magnetic field amplitude of geomagnetic station

图4 感应地电场相对值Fig.4 Relative value of induced geoelectric field

由图4可知，在2004年11月7日18：00：00—19：00：00时段存在感应地电场最大值，与GIC最大值出现时刻对应。在评估及预防地磁暴风险灾害时，其中最关心的指标之一就是感应地电场的最大瞬时值，故选取了包含最大值脉冲的18：24：57—18：42：00时段计算感应地电场的准确值。各地磁台该时段的磁场数据如图5所示。

2 基于电磁测深的感应地电场计算方法

2.1 磁暴期间磁场数据的选取与处理

已有4个地磁台包含感应地电场最大值脉冲的磁场时域数据，经过快速傅里叶变换(FFT)可得到各地磁台频域磁场数据，如图6所示。

图5 各地磁台18：24：57—18：42：00时段磁场幅值线Fig.5 Magnetic field amplitude line of local magnetic stations from 18：24：57 to 18：42：00

f为频率图6 各地磁台频域磁场幅值Fig.6 Frequency domain magnetic field amplitude of local magnetic stations

由图6可以看出，所选取的地磁场数据的主要频率范围在GIC频率范围内，即10-4～10-2Hz，第一个频率点对应为1/1 024 Hz。

2.2 感应地电场频域特性

采集2.1节得到的测点磁场频率范围，将测点视电阻率在磁场包含频率上作线性插值后计算得到对应阻抗，将阻抗与对应幅值相乘可得到感应地电场频域幅值。相位关系依照式(5)和式(6)，未作图展示。不用考虑直流分量，因为直流对应的阻抗恒为0，不影响计算结果。

同一地区不同测点出现最大幅值对应的频率基本相同，如图7所示。阿尔山地区感应地电场Y方向分量最大值为6.836×10-3Hz对应的0.121 9 V/km，X方向分量最大值为3.906×10-3Hz对应的3.467×10-2V/km；阿拉善盟感应地电场Y方向分量最大值为1.953×10-3Hz对应的1.245×10-2V/km，X方向分量最大值为1.953×10-3Hz对应的2.64×10-2V/km；大别山造山带感应地电场Y方向分量最大值为1.953×10-3Hz对应的2.824×10-2V/km，X方向分量同地区感应地电场频域最大幅值差别较大。

图7 各测点频域特性Fig.7 Frequency domain electric field amplitude of each measuring point

2.3 感应地电场时域特性

对2.2节中得到的各测点频域幅值进行快速傅里叶反变换可得到地电场幅值如图8所示。

2.4 时域感应地电场波形处理

由于在傅里叶变换中对数据进行了矩形截断，会产生截断误差，即数据两端会出现“尖刺”。对数据两端分别截去5%的数据点，处理后重新作图，结果如图9所示。处理后的时域特性消除了数据的截断误差，时段内数据变化更加清晰，更准确地反映了感应地电场时域特性。

图8 各测点地电场时域特性Fig.8 Time domain electric field amplitude of each measuring point

图9 波形处理后各测点时域特性Fig.9 Time domain electric field amplitude after correction of each measuring point

3 结论

(1)地磁暴对各地地电场的影响不均匀，与当地的地质电性结构有很大的关系。在计算广东增城朱村、东莞谢岗的感应地电场时，其都使用了同一纬度上距离最近的广州GZH地磁台地磁数据，但增城朱村南北方向分量幅值最大值为0.932 0 V/km，而东莞谢岗南北方向分量幅值方向最大值仅为0.062 1 V/km，相差十倍以上。

(2)磁暴发生期间，所有测点的南北方向分量幅值在0.031 5～0.932 0 V/km范围内，东西方向分量幅值在0.015 6～0.109 3 V/km范围内。可见感应地电场南北方向分量幅值普遍大于东西方向分量幅值，说明南北走向的电网将产生更大的GIC，更易收到地磁暴的侵害，应作为主要的关注对象。

(3)采用本文所提出的基于地磁测深的感应地电场计算方法，仅需要1 024 s的磁场数据就可以准确得出关心时段内的感应地电场最大值，对于本文所选时段，用512 s的磁场数据也可以得到准确的计算结果。仅利用包含最大值的脉冲计算，可以在不影响准确度的前提下，极大缩短磁暴发生时感应地电场的计算时间，减少计算所需内存，对应对地磁暴灾害有着重要意义。