王澤忠，司 遠，劉連光，張彩友，姚 輝，陳水耀

（1. 華北電力大學 高電壓與電磁兼容北京重點實驗室，北京102206；2. 華北電力大學 新能源電力系統(tǒng)國家重點實驗室，北京102206；3. 國網(wǎng)浙江省電力有限公司檢修分公司，浙江 杭州310000）

0 引言

電網(wǎng)地磁擾動（GMD）災害防治有設(shè)備防治、調(diào)度防治2種方法。設(shè)備防治方法受限于資金投入，無法大量使用［1］；調(diào)度防治方法需利用電網(wǎng)運行調(diào)度及時補償?shù)卮鸥袘?yīng)電流（GIC）引起的無功損耗［2］，因此研究計算精度高、速度快的GIC計算方法成為關(guān)鍵。

GMD 地電場是計算GIC 的關(guān)鍵，傳統(tǒng)的計算GMD 地電場的方法首先需建立大地電導率模型，均勻［3］、水平分層［4］及分層分區(qū)大地電導率模型［5］都是基于土壤電導率數(shù)據(jù)，在建模過程中都會遇到實測電導率數(shù)據(jù)不足的問題，因此需要一種既考慮地質(zhì)變量的隨機性又考慮其自相關(guān)性的散亂點插值方法。P. A. Fernberg 在研究地磁暴擾動下管地電位分布時未將重點放在建立大地電導率模型上，而是從具體的測點入手，提出基于各測點的視電阻率數(shù)據(jù)研究發(fā)生地磁暴時管道中感應(yīng)電流的思路［6］，整個過程省略了建立大地電導率模型的復雜過程。但對于計算GMD 地電場所需地磁場數(shù)據(jù)，只是使用了其假設(shè)值，而未使用地磁場的實測數(shù)據(jù)，所以該方法能否用于電網(wǎng)GIC的計算仍是需要研究的問題。

關(guān)于GMD 地電場的計算方法，平面波法［7］、分區(qū)平面波法［8］認為整個研究區(qū)域或某塊區(qū)域的GMD地電場是相同的，沒有充分考慮GMD 地電場的實際分布情況。有限元法［9］的基礎(chǔ)是準確的分層分區(qū)三維大地電導率模型。大地模型中存在尺度較小的突變區(qū)域，而有限元計算采用自適應(yīng)網(wǎng)格劃分，要求網(wǎng)格尺寸小于模型中最小的單元，使網(wǎng)格數(shù)量增加，所需存儲空間大，計算成本過高。此外，由電網(wǎng)GIC 計算中等效直流電壓源的計算方法［10-11］可知：關(guān)于地磁暴對電網(wǎng)的影響，真正需要關(guān)心的是有輸電線路區(qū)域的感應(yīng)地電場，即準確計算變壓器接地點的電位才是關(guān)鍵。而使用有限元法研究地電場時需要對整個區(qū)域進行計算，不能做到針對性地計算指定區(qū)域的感應(yīng)電場。

本文借鑒P. A. Fernberg等人的研究成果，對廣東500 kV 輸電線路進行微元化處理，將微元中點作為待估點，將實測視電阻率作為基礎(chǔ)數(shù)據(jù)，將地質(zhì)統(tǒng)計學中的空間-頻率變異函數(shù)［12］引入視電阻率估算，并結(jié)合2004 年11 月9 日地磁暴廣東肇慶地磁臺監(jiān)測到的GMD 數(shù)據(jù)直接計算微元GMD 地電場，省略了復雜的大地建模過程，利用微元地電場和輸電線路方向矢量做內(nèi)積運算得到等效電壓源，進而計算電網(wǎng)的GIC，以期為運行調(diào)度防治地磁暴災害提供基礎(chǔ)算法。

1 空間?頻率變異函數(shù)

視電阻率并不是指某巖層真實的電阻率，而是電場分量和磁場分量的比值，在大地電磁測深MT（MagnetoTelluric sounding）中能夠反映地下電性結(jié)構(gòu)［9］。利用變異函數(shù)對未知點的視電阻率進行估計，由于視電阻率與測量位置和場源頻率有關(guān)，故將空間變異函數(shù)拓展為空間-頻率變異函數(shù)。

1.1 假設(shè)條件

當區(qū)域化變量Z(xi，fi)滿足以下2個條件時為二階平穩(wěn)［12］。

（1）在研究區(qū)域內(nèi)，Z(xi，fi)的數(shù)學期望存在，且為常數(shù)，即：

其中，m 為常數(shù)；E[?]為數(shù)學期望；xi為某一點的空間位置；fi為某一點的頻率。

（2）Z(xi，fi)的協(xié)方差函數(shù)存在且相同，不依賴空間絕對位置，而依賴相對位置（只依賴滯后距hs、hω，而與xi、fi無關(guān)），即：

其中，Ccov{·，·}為協(xié)方差函數(shù)；hs為空間距離；hω為頻率間隔；Z(xi+hs，fi+hω)為與Z(xi，fi)空間距離相距hs、頻率相距hω的實測點；C(hs，hω)為僅與空間滯后距離和頻率間隔有關(guān)的協(xié)方差函數(shù)。

當Z(xi，fi)的 增 量Z(xi，fi)?Z(xi+hs，fi+hω)滿足以下2個條件時，稱Z(xi，fi)滿足本征假設(shè)。

（1）Z(xi，fi)的增量的數(shù)學期望存在，且為常數(shù)，即：

（2）增量Z(xi，fi)?Z(xi+hs，fi+hω)的方差函數(shù)（變異函數(shù)）存在且平穩(wěn)，即：

其中，Vvar(·)為方差函數(shù)；v(hs，hω)為空間-頻率變異函數(shù)。

1.2 實驗變異函數(shù)的估計

將視電阻率Z(xi，fi)作為區(qū)域化變量，其空間-頻率變異函數(shù)定義為：

其中，x 為視電阻率實測點的位置坐標；f 為視電阻率測量時的頻率。

利用所有實測點的視電阻率計算空間-頻率變異函數(shù)v(hs，hω)的估計值v?(hs，hω)，見式（6）。

其中，n(hs，hω)為實測點對數(shù)量。

1.3 變異函數(shù)理論模型

利用估計值分別繪制實驗空間、頻率變異函數(shù)的散點圖，然后采用式（7）所示的球形變異函數(shù)擬合變異函數(shù)均值數(shù)據(jù)點，得到理論變異函數(shù)曲線參數(shù)。

其中，c為偏基臺值；a為變程；c0為塊金值；c+c0為基臺值為空間／頻率變異函數(shù)。

將變異函數(shù)拓展到空間-頻域，常用的積和式空間-頻率變異函數(shù)v(hs，hω)如式（8）所示。

其中，Csω(0，0)、Cs(0)、Cω(0)為3 種變異函數(shù)的基臺值；k1—k3為系數(shù)。

2 基于變異函數(shù)的視電阻率估計算法

根據(jù)線路兩端變電站的經(jīng)緯度坐標，對線路進行微元化處理，將線路平均分為多個線段，得到每個微元中點的地理坐標，將其作為視電阻率的待估點。地理屬性既具有空間相關(guān)性，又具有隨機性。根據(jù)地理學第一定律及地質(zhì)統(tǒng)計學中的Kriging 法，利用所有已知點的數(shù)據(jù)加權(quán)估計微元中點的視電阻率，其中權(quán)重系數(shù)需滿足無偏性、估計誤差最小的要求。

2.1 無偏性

設(shè)某一視電阻率實測點x（ii=1，2，…，n，n 為實測點數(shù)量）在某一頻率fi下的值為Zi(xi，fi)，則待估點的視電阻率Z0(x，f)的預測結(jié)果定義為已知實測點視電阻率Zi(xi，fi)的加權(quán)和，如式（9）所示。

其中，λi為待定權(quán)重系數(shù)。

進而可得到：

2.2 估計誤差最小

為了使視電阻率估計誤差最小，對估計誤差進行分析。

其中，σ為標準差。

利用Lagrange 乘數(shù)法進一步整理式（12）可得式（13），具體推導過程見附錄A式（A1）—（A9）。

其中，vi，j(hs，hω)=v(Zi，Zj)=v(Zi?Zj)。

求解式（13），得到一組加權(quán)系數(shù)λ1、λ2、…、λn，將其代入式（9）估算待估點的視電阻率。

3 基于微元法的GMD地電場計算

根據(jù)每條線路兩端變電站的經(jīng)緯度坐標，對線路做微元化處理，得到每個微元中點的地理坐標，將這些點作為地電場的計算區(qū)域。將時域地磁場數(shù)據(jù)經(jīng)快速傅里葉變換（FFT）得到頻域數(shù)據(jù)，利用第2節(jié)的方法計算微元中點對應(yīng)地理坐標的視電阻率?；谝曤娮杪师?、波阻抗間的關(guān)系，根據(jù)式（14）計算微元中點GMD 地電場的南北向分量Exi、東西向分量Eyi，然后進行快速傅里葉逆變換（IFFT）即可得到每個線路微元中點的時域微元地電場分量。

其中，Hy、Hx分別為地磁場東西、南北向分量；μ0為真空的磁導率；ρi為微元i中點的視電阻率；ω為角頻率。

在電網(wǎng)GIC 的模型中，電網(wǎng)受地磁暴影響產(chǎn)生GIC 的原理相當于在不同的變電站間施加直流電壓源［13］。利用每條線路兩端變電站的經(jīng)緯度計算每條線路單位長度方向矢量，微元地電場Exi、Eyi與相應(yīng)的線路方向矢量做內(nèi)積并進行疊加得到等效電壓源V，可將等效電壓源的積分公式簡化為式（15）。

具體計算流程如圖1所示。

4 算例分析

圖1 計算流程Fig.1 Computation process

本文采用地磁臺地磁場測量數(shù)據(jù)計算GMD 地電場，結(jié)合廣東500 kV 電網(wǎng)結(jié)構(gòu)參數(shù)計算GIC，同時利用文獻［5］中的方法建立大地電導率模型并利用有限元法進行計算，然后對比分析2 種方法的計算結(jié)果。

4.1 基礎(chǔ)數(shù)據(jù)

廣東省地震局與國家地震局地質(zhì)研究所在廣東增城—深圳剖面的增城朱村、東莞謝崗、寶安龍崗、深圳大鵬等地開展了大地電磁測深研究，將頻率在10-4～10-1Hz范圍內(nèi)的實測視電阻率作為原始數(shù)據(jù)。

視電阻率實測點有14個，每個實測點對應(yīng)20個頻點，由于個別實測點缺失頻點數(shù)據(jù)，總的實測數(shù)據(jù)不少于14×20 個。利用ArcGIS 中的Geostatistical Analyst 工具箱進行數(shù)據(jù)檢驗，并對數(shù)據(jù)進行對數(shù)變換，處理后的數(shù)據(jù)圖形接近直線（見附錄B 圖B1），可見數(shù)據(jù)服從正態(tài)分布，即數(shù)學期望和方差都存在，滿足二階平穩(wěn)假設(shè)。

本文算例采用肇慶地磁臺在2004 年11 月9 日監(jiān)測到的GMD 數(shù)據(jù)計算廣東500 kV 電網(wǎng)的GIC。根據(jù)廣東電網(wǎng)數(shù)據(jù)和資料，廣東500 kV 電網(wǎng)共有54座變電站和62 條線路。各站點的經(jīng)緯度和線路的具體參數(shù)參見文獻［14］，根據(jù)數(shù)據(jù)構(gòu)建電網(wǎng)等效模型［15］從而計算GIC。

4.2 微元中點視電阻率計算

將每條線路分為100 個微元，得到各微元中點的地理坐標。用變異函數(shù)計算各微元對應(yīng)地理坐標的視電阻率，實驗空間、頻率變異函數(shù)散點圖分別見圖2（a）、（b），利用式（7）球形理論變異函數(shù)擬合所得理論空間、頻率變異函數(shù)曲線分別見圖2（c）、（d）。

根據(jù)實驗變異函數(shù)散點擬合得到的球形理論空間變異函數(shù)、球形理論頻率變異函數(shù)分別如式（16）和式（17）所示。

圖2 空間、頻率變異函數(shù)擬合結(jié)果Fig.2 Fitting results of spatial and frequency variation functions

圖3 空間-頻率變異函數(shù)的實驗模型Fig.3 Experimental model of spatial-frequency variation function

根據(jù)式（6）建立空間-頻率變異函數(shù)的實驗模型，結(jié)果如圖3所示。擬合得到的空間變異函數(shù)基臺值Cs(0)=2.73，頻率變異函數(shù)基臺值Cω(0)=4.1，實驗空間-頻率變異函數(shù)的基臺值Csω(0，0)=3.25，將其代入式（8）得到k1=0.32，k2=?0.311，k3=0.13，則理論空間-頻率變異函數(shù)如式（18）所示。

得到變異函數(shù)后，將微元中點作為預測點，計算權(quán)重系數(shù)，進而得到待估點的視電阻率。

4.3 微元GMD地電場計算

計算得到各微元中點的視電阻率后，將肇慶地磁臺的地磁場時域數(shù)據(jù)進行FFT 得到頻域數(shù)據(jù)，結(jié)合式（14）計算得到各微元的GMD 地電場，然后進行IFFT 得到時域下微元的GMD 地電場［16］。廣東地區(qū)2004年11月9日22:48:00時刻存在輸電線路區(qū)域的GMD 地電場見附錄B 圖B2。圖中每個矢量的方向代表微元GMD 地電場的方向，長短表示微元GMD地電場的大小，可以看出各微元的GMD 地電場是不同的，從而證明了在利用GMD 地電場沿線路積分求等效直流電壓源時，將某一區(qū)塊或者某一層的GMD地電場認為是同一值是不符合實際情況的。

計算得到肇慶站在22:48:00—22:52:00 時段內(nèi)的GMD 地電場分量Ex、Ey，結(jié)果分別見圖4（a）、（b）。根據(jù)式（15）計算等效電壓源，并結(jié)合廣東地區(qū)的電網(wǎng)參數(shù)計算得到嶺澳站的GIC，結(jié)果見圖4（c）。

圖4 地電場分量和GIC的計算結(jié)果Fig.4 Calculative results of geoelectric field components and GIC

4.4 對比分析

由于可查到肇慶地磁臺的地電場監(jiān)測數(shù)據(jù)、嶺澳站的GIC 監(jiān)測數(shù)據(jù)，選擇這2個站點進行2種算法的準確度比較。將計算所得肇慶站的300個GMD地電場數(shù)據(jù)、嶺澳站的300 個GIC 數(shù)據(jù)與實測值進行對比分析，結(jié)果見表1（表中數(shù)據(jù)樣本數(shù)量為300）。

表1 2種算法的計算結(jié)果與實測值對比Table 1 Comparison between calculative results of two algorithms and measured values

由表1 可知，本文所提算法的計算精度較高，但所提算法的計算結(jié)果與實測值之間仍然存在誤差，主要是因為以下2 點原因：①計算微元中點GMD 地電場時均采用肇慶地磁臺提供的地磁場數(shù)據(jù)，沒有考慮地磁場的實際分布；②只是將每條線路分為100個微元，各微元過長。如果將線路分成更多的微元，并綜合多個地磁臺的地磁場測量數(shù)據(jù)，通過插值得到不同微元對應(yīng)地理位置的地磁場數(shù)據(jù)，并將其應(yīng)用到GMD地電場計算中可進一步提高計算精度。

在相同的計算機上完成2 種算法用于對比其計算速度。當采用有限元法計算GMD 地電場時，大地電導率模型的南北邊界范圍為600 km，東西邊界范圍為800 km，深度為300 km，網(wǎng)格總數(shù)約為107個，整個有限元法占用內(nèi)存約為45 152 961 MByte，計算時間為9 385 s。此外，還需考慮在ANSYS 中進行復雜的分層分區(qū)大地電導率建模過程所用時間和內(nèi)存。相比于先建立分層分區(qū)大地電導率模型后采用有限元法計算GMD 地電場的方法，本文算法省略了復雜的大地建模過程，同時使用空間-頻率變異函數(shù)結(jié)合滑動加權(quán)平均法進行視電阻率散點插值，利用視電阻率與地磁場的關(guān)系代替復雜的有限元計算，不再以大地模型為導向，而是以輸電線路所在位置為出發(fā)點，進行針對性的計算，僅需計算存在輸電線路區(qū)域的GMD 地電場，即只需要計算6 200 個微元中點的地電場，同時考慮微元中點的視電阻率散點插值過程，共占用內(nèi)存130 426 MByte，計算時間為1 082 s?？梢?，即使考慮為了進一步提高計算精度而將線路分為更多的微元，但所分微元的數(shù)量相比于有限元的網(wǎng)格總數(shù)仍很少，所以基于視電阻率和微元法的GIC 計算方法可節(jié)省計算時間和內(nèi)存，適用于電網(wǎng)調(diào)度的實時計算。

5 結(jié)論

（1）本文首次將地質(zhì)統(tǒng)計學中的變異函數(shù)引入GMD 地電場計算中，從線路角度出發(fā)，將研究區(qū)域內(nèi)實測視電阻率作為原始數(shù)據(jù)，利用變異函數(shù)求解線路微元中點對應(yīng)的視電阻率。結(jié)合地磁臺地磁場的測量數(shù)據(jù)計算每個微元的地電場，GIC 和GMD 地電場的計算結(jié)果與實際監(jiān)測結(jié)果之間的擬合程度更高、更準確。

（2）使用地磁場實測數(shù)據(jù)和視電阻率實測數(shù)據(jù)，為計算GMD 地電場提供了一種非常簡單的方法，省略了復雜的大范圍大地電性結(jié)構(gòu)建模的過程。此外，本文所提算法只需研究存在輸電線路區(qū)域的GMD 地電場就可以計算不同站點之間的等效直流電壓源，提高了計算速度。

（3）各微元的地磁場數(shù)據(jù)使用的均為同一地磁臺的地磁場測量數(shù)據(jù)，仍存在一定的誤差。若能根據(jù)太陽風的參數(shù)、數(shù)據(jù)，結(jié)合全球三維磁層數(shù)值模型實時提供各微元的地磁場擾動數(shù)據(jù)，就可以用于進一步提高本文算法的精確度，有利于GIC 事故高風險點GIC預測及運行調(diào)度的GIC事故在線分析。

附錄見本刊網(wǎng)絡(luò)版（http：//www.epae.cn）。