地磁感應(yīng)地電場時頻變換快速計算方法

倪政澤， 王澤忠， 李宇妍， 司遠(yuǎn)

(華北電力大學(xué)高電壓與電磁兼容北京市重點實驗室， 北京 102206)

地磁暴是指由于太陽活動而造成地球磁場全球性劇烈擾動的現(xiàn)象。地球磁場的劇烈擾動將在地面中感應(yīng)出電場，該感應(yīng)電場會在由輸電線路，中性點接地的變壓器及大地構(gòu)成的回路中產(chǎn)生地磁感應(yīng)電流(geomagnetically induced current，GIC)，GIC頻率在0.01～0.000 1 Hz范圍內(nèi)變化，這種準(zhǔn)直流電流對變壓器等電網(wǎng)設(shè)備會造成嚴(yán)重威脅[1-4]。同時地磁場變化迅速，如何準(zhǔn)確快速地求取磁暴發(fā)生時的地面感應(yīng)地電場是防治磁暴災(zāi)害的關(guān)鍵。

傳統(tǒng)的地電場計算方法需要根據(jù)不同地區(qū)土壤的電性結(jié)構(gòu)確定大地電導(dǎo)率模型，文獻(xiàn)[5]提出建立均勻大地模型，沒有充分考慮地電場的實際分布情況。文獻(xiàn)[6]提出建立水平分層大地電導(dǎo)率模型，利用平面波法[7]計算地電場，對于中國特高壓長距離輸電線路，其空間尺度達(dá)到幾百公里甚至上千公里，地質(zhì)結(jié)構(gòu)復(fù)雜[8]，各地區(qū)土壤結(jié)構(gòu)差異較大，難以建立準(zhǔn)確模型。文獻(xiàn)[9-10]建立了包含海陸分界面的二維大地電導(dǎo)率模型，但沒有電導(dǎo)率的實測數(shù)據(jù)，文獻(xiàn)[11-12]提出建立三維大地分層分區(qū)電導(dǎo)率模型，利用地磁臺站記錄的地磁場變化量作為地面邊界條件，采用有限元法計算地電場，但自適應(yīng)剖分網(wǎng)格過多，計算時間過長。無法準(zhǔn)確得出各地大地模型并快速計算成為傳統(tǒng)地電場方法的局限性。文獻(xiàn)[13]提出時諧擬合算法將時變的地磁場轉(zhuǎn)化為頻域相量法計算，數(shù)據(jù)擬合得到的僅為單一頻率，與實際相差較大。文獻(xiàn)[14]提出利用視電阻率數(shù)據(jù)計算磁暴感應(yīng)地電場頻域幅值，無法準(zhǔn)確選取包含磁暴最大幅值的脈沖，計算結(jié)果準(zhǔn)確性受數(shù)據(jù)影響較大。

地磁測深是利用天然交變電磁場研究地球電性結(jié)構(gòu)的一種物理勘探方法，具有成本低，工作方便，不受高阻層屏蔽的優(yōu)點。在地下巖石電性分布不均勻(有兩種或兩種以上導(dǎo)電性不同的巖石或礦石)或地表起伏不平的情況下，仍按測定均勻水平大地電阻率的方法和計算公式求得的電阻率稱之為視電阻率，其綜合反映了大地的電性結(jié)構(gòu)，是地磁測深法的重要數(shù)據(jù)。現(xiàn)利用2004年11月7日地磁臺實測秒級數(shù)據(jù)與視電阻率數(shù)據(jù)在頻域計算得到頻域電場值，再利用傅里葉反變換得到時域電場值。利用這種方法對地磁暴發(fā)生時刻阿爾山地區(qū)，阿拉善盟，大別山地區(qū)，廣東東莞謝崗及增城朱村的感應(yīng)地電場最大值進行對比。該計算方法簡化了地面感應(yīng)地電場的計算，為工程應(yīng)用提供了有效算法，在地磁暴災(zāi)害防控中有著重要意義。

1 視電阻率在地電場算法中的應(yīng)用

在大地電磁測深方法中，考慮到應(yīng)用的觀測頻率范圍以及構(gòu)成地殼淺部介質(zhì)的電阻率范圍，在大地介質(zhì)中可忽略位移電流對場分布的影響[15]，于是諧變場的Maxwell方程組表示為

(1)

(2)

(3)

(4)

式中：H為地磁場強度向量；E為地電場強度向量；ω為角頻率；μ為磁導(dǎo)率；γ為電導(dǎo)率；i為虛數(shù)單位。

如圖1所示建立笛卡爾坐標(biāo)系。地磁暴對地電場的影響可看作是平面電磁波垂直入射大地介質(zhì)(即OZ方向)中產(chǎn)生的電磁場，同時對于電磁測深方法，視電阻率ρ的測定以均勻水平大地為前提。

O為觀測點(坐標(biāo)原點)；OXY為水平面；OX指向地磁北；OY指向地磁東；OZ為鉛直向下；N為地磁北方向；E為地磁東方向圖1 地磁場坐標(biāo)系的定義Fig.1 Definition of geomagnetic coordinate system

根據(jù)均勻介質(zhì)中的平面波理論[16]，有

(5)

(6)

式中：Zxy為TE模式下的波阻抗；Zyx為TM模式下的波阻抗；Ex和Ey分別為地電場南北分量和東西分量；Hx和Hy分別為地磁場南北分量和東西分量；ρxy為TE模式下的視電阻率；ρyx為TM模式下的視電阻率；i為虛數(shù)單位；μ為磁導(dǎo)率，大地磁導(dǎo)率近似為μ0。

若只考慮視電阻率與對應(yīng)波阻抗的幅值，有以下關(guān)系：

(7)

(8)

目前已有通過地磁測深得到的視電阻率數(shù)據(jù)，可通過式(7)和式(8)得到對應(yīng)阻抗Zxy和Zyx。將時域地磁場數(shù)據(jù)經(jīng)快速傅里葉變換(FFT)后得到頻域磁場數(shù)據(jù)，利用式(5)和式(6)的關(guān)系計算地電場Ex和Ey的頻域數(shù)值，然后進行快速傅里葉反變換(IFFT)即可得到地電場Ex、Ey的時域數(shù)值。該計算方法以視電阻率數(shù)據(jù)作為支撐，將時域中多點求解問題轉(zhuǎn)化為頻域中的單點求解問題。

1.1 磁場數(shù)據(jù)的來源

目前收集到地磁臺磁暴數(shù)據(jù)有分鐘級和秒級地磁場三分量(即地磁場總強度F、磁偏角D和地磁場水平分量H)數(shù)據(jù)，為了得到地磁南北和東西分量變化率脈沖幅值，現(xiàn)使用秒級數(shù)據(jù)。所收集到的數(shù)據(jù)中存在一些因設(shè)備故障而導(dǎo)致的錯誤數(shù)據(jù)，剔除錯誤數(shù)據(jù)后，能夠使用的據(jù)包括2004年11月7—9日黑龍江德都(DED)、廣東廣州(GZH)、甘肅嘉峪關(guān)(JYG)、海南瓊中(QGZ)、湖北武漢(WHN)、新疆烏魯木齊(WMQ)和新疆喀什(KSH)7個地磁臺的測量數(shù)據(jù)。各臺站的地理坐標(biāo)如表1所示。

表1 地磁臺站地理坐標(biāo)Table 1 Geographic coordinates of geomagnetic stations

1.2 視電阻率數(shù)據(jù)的選取與處理

選取的阿爾山火山區(qū)大地電磁測深數(shù)據(jù)由文獻(xiàn)[17]中獲得，觀測了一條西北向測線上的7個大地電磁測深點，選取其中編號3、4兩個測點；阿拉善盟的數(shù)據(jù)由文獻(xiàn)[18]中獲得，其觀測了東西向測線共65個測點，選取其中編號43、58兩個測點；大別山的數(shù)據(jù)由文獻(xiàn)[19]中獲得，其觀測了自西向東3條大地電磁測深側(cè)線，測線LA沿NE26°方向，測線LB沿著NE28°方向，選取測點LA7、測點LB7兩個測點；廣東增城朱村，東莞謝崗的數(shù)據(jù)由廣東省地震局與國家地震局地質(zhì)研究所開展的大地電磁測深研究中獲得，其電阻率幅值曲線如圖2所示。

1.3 磁暴期間磁場數(shù)據(jù)的選取與處理

為了提高計算的準(zhǔn)確度和精度，根據(jù)選取同一緯度范圍內(nèi)，距離最近的地磁臺站作為地磁場數(shù)據(jù)來源的準(zhǔn)則，選取德都地磁臺(DED)，嘉峪關(guān)地磁臺(JYG)，武漢地磁臺(WHN)，廣州地磁臺(GZH)的秒級磁場數(shù)據(jù)分別作為阿爾山地區(qū)，阿拉善盟，大別山造山帶，廣東增城朱村與東莞謝崗的磁場數(shù)據(jù)來源。

在2004年11月7日18：28：08的磁暴急始時刻，監(jiān)測廣東嶺澳核電站內(nèi)的中性點電流峰值達(dá)到了-47.2 A(負(fù)號表示電流從大地流入中性點)[20]，故選取2004年11月7日17：00—21：00各地磁臺磁場數(shù)據(jù)作為原始數(shù)據(jù)進行分析，各地磁臺該時段的磁場數(shù)據(jù)如圖3所示，縱軸單位為nT，1T=1×109nT。

圖2 測點視電阻率幅值曲線Fig.2 Apparent resistivity amplitude curve of measuring point

地面感應(yīng)電場與地磁場存在半階導(dǎo)數(shù)關(guān)系[21]，即

(9)

(10)

式中：Ex(ω)、Ey(ω)為感應(yīng)地電場頻域的x、y分量；Bx(ω)、By(ω)地磁場頻域的x、y分量；j為虛數(shù)。其對應(yīng)的時域內(nèi)地電場分量和對應(yīng)地磁場分量的關(guān)系式：

(11)

(12)

式中：u為變量名；t為待求電場的時間；gx(u)=dBx/dt、gy(u)=dBy/dt為地磁場時域的x、y分量變化率；μ0為真空磁導(dǎo)率；γ為大地電導(dǎo)率。將式(11)和式(12)歸一離散化后，得

(13)

式(13)中：bn=Bn-Bn-1為地磁場分量的一階差分；M為一足夠大的數(shù)；TN為待求電場的第N個時刻；Tn為第n個采樣點的時間。由式(13)即可由磁場數(shù)據(jù)的時間序列計算出地面電場相對值的時間序列，2004年11月7日17：00—21：00各地臺感應(yīng)地電場相對值如圖4所示。

圖3 各地磁臺磁場幅值Fig.3 Magnetic field amplitude of geomagnetic station

圖4 感應(yīng)地電場相對值Fig.4 Relative value of induced geoelectric field

由圖4可知，在2004年11月7日18：00：00—19：00：00時段存在感應(yīng)地電場最大值，與GIC最大值出現(xiàn)時刻對應(yīng)。在評估及預(yù)防地磁暴風(fēng)險災(zāi)害時，其中最關(guān)心的指標(biāo)之一就是感應(yīng)地電場的最大瞬時值，故選取了包含最大值脈沖的18：24：57—18：42：00時段計算感應(yīng)地電場的準(zhǔn)確值。各地磁臺該時段的磁場數(shù)據(jù)如圖5所示。

2 基于電磁測深的感應(yīng)地電場計算方法

2.1 磁暴期間磁場數(shù)據(jù)的選取與處理

已有4個地磁臺包含感應(yīng)地電場最大值脈沖的磁場時域數(shù)據(jù)，經(jīng)過快速傅里葉變換(FFT)可得到各地磁臺頻域磁場數(shù)據(jù)，如圖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可以看出，所選取的地磁場數(shù)據(jù)的主要頻率范圍在GIC頻率范圍內(nèi)，即10-4～10-2Hz，第一個頻率點對應(yīng)為1/1 024 Hz。

2.2 感應(yīng)地電場頻域特性

采集2.1節(jié)得到的測點磁場頻率范圍，將測點視電阻率在磁場包含頻率上作線性插值后計算得到對應(yīng)阻抗，將阻抗與對應(yīng)幅值相乘可得到感應(yīng)地電場頻域幅值。相位關(guān)系依照式(5)和式(6)，未作圖展示。不用考慮直流分量，因為直流對應(yīng)的阻抗恒為0，不影響計算結(jié)果。

同一地區(qū)不同測點出現(xiàn)最大幅值對應(yīng)的頻率基本相同，如圖7所示。阿爾山地區(qū)感應(yīng)地電場Y方向分量最大值為6.836×10-3Hz對應(yīng)的0.121 9 V/km，X方向分量最大值為3.906×10-3Hz對應(yīng)的3.467×10-2V/km；阿拉善盟感應(yīng)地電場Y方向分量最大值為1.953×10-3Hz對應(yīng)的1.245×10-2V/km，X方向分量最大值為1.953×10-3Hz對應(yīng)的2.64×10-2V/km；大別山造山帶感應(yīng)地電場Y方向分量最大值為1.953×10-3Hz對應(yīng)的2.824×10-2V/km，X方向分量同地區(qū)感應(yīng)地電場頻域最大幅值差別較大。

圖7 各測點頻域特性Fig.7 Frequency domain electric field amplitude of each measuring point

2.3 感應(yīng)地電場時域特性

對2.2節(jié)中得到的各測點頻域幅值進行快速傅里葉反變換可得到地電場幅值如圖8所示。

2.4 時域感應(yīng)地電場波形處理

由于在傅里葉變換中對數(shù)據(jù)進行了矩形截斷，會產(chǎn)生截斷誤差，即數(shù)據(jù)兩端會出現(xiàn)“尖刺”。對數(shù)據(jù)兩端分別截去5%的數(shù)據(jù)點，處理后重新作圖，結(jié)果如圖9所示。處理后的時域特性消除了數(shù)據(jù)的截斷誤差，時段內(nèi)數(shù)據(jù)變化更加清晰，更準(zhǔn)確地反映了感應(yīng)地電場時域特性。

圖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 結(jié)論

(1)地磁暴對各地地電場的影響不均勻，與當(dāng)?shù)氐牡刭|(zhì)電性結(jié)構(gòu)有很大的關(guān)系。在計算廣東增城朱村、東莞謝崗的感應(yīng)地電場時，其都使用了同一緯度上距離最近的廣州GZH地磁臺地磁數(shù)據(jù)，但增城朱村南北方向分量幅值最大值為0.932 0 V/km，而東莞謝崗南北方向分量幅值方向最大值僅為0.062 1 V/km，相差十倍以上。

(2)磁暴發(fā)生期間，所有測點的南北方向分量幅值在0.031 5～0.932 0 V/km范圍內(nèi)，東西方向分量幅值在0.015 6～0.109 3 V/km范圍內(nèi)。可見感應(yīng)地電場南北方向分量幅值普遍大于東西方向分量幅值，說明南北走向的電網(wǎng)將產(chǎn)生更大的GIC，更易收到地磁暴的侵害，應(yīng)作為主要的關(guān)注對象。

(3)采用本文所提出的基于地磁測深的感應(yīng)地電場計算方法，僅需要1 024 s的磁場數(shù)據(jù)就可以準(zhǔn)確得出關(guān)心時段內(nèi)的感應(yīng)地電場最大值，對于本文所選時段，用512 s的磁場數(shù)據(jù)也可以得到準(zhǔn)確的計算結(jié)果。僅利用包含最大值的脈沖計算，可以在不影響準(zhǔn)確度的前提下，極大縮短磁暴發(fā)生時感應(yīng)地電場的計算時間，減少計算所需內(nèi)存，對應(yīng)對地磁暴災(zāi)害有著重要意義。