王澤忠，司 遠(yuǎn)，劉連光，張彩友，姚 輝，陳水耀

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

0 引言

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

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

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

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

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

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

1.1 假設(shè)條件

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

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

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

（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ω的實測點(diǎn)；C(hs，hω)為僅與空間滯后距離和頻率間隔有關(guān)的協(xié)方差函數(shù)。

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

（1）Z(xi，fi)的增量的數(shù)學(xué)期望存在，且為常數(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 為視電阻率實測點(diǎn)的位置坐標(biāo)；f 為視電阻率測量時的頻率。

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

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

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

利用估計值分別繪制實驗空間、頻率變異函數(shù)的散點(diǎn)圖，然后采用式（7）所示的球形變異函數(shù)擬合變異函數(shù)均值數(shù)據(jù)點(diǎn)，得到理論變異函數(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)緯度坐標(biāo)，對線路進(jìn)行微元化處理，將線路平均分為多個線段，得到每個微元中點(diǎn)的地理坐標(biāo)，將其作為視電阻率的待估點(diǎn)。地理屬性既具有空間相關(guān)性，又具有隨機(jī)性。根據(jù)地理學(xué)第一定律及地質(zhì)統(tǒng)計學(xué)中的Kriging 法，利用所有已知點(diǎn)的數(shù)據(jù)加權(quán)估計微元中點(diǎn)的視電阻率，其中權(quán)重系數(shù)需滿足無偏性、估計誤差最小的要求。

2.1 無偏性

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

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

進(jìn)而可得到：

2.2 估計誤差最小

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

其中，σ為標(biāo)準(zhǔn)差。

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

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

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

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

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

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

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

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

4 算例分析

圖1 計算流程Fig.1 Computation process

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

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

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

視電阻率實測點(diǎn)有14個，每個實測點(diǎn)對應(yīng)20個頻點(diǎn)，由于個別實測點(diǎn)缺失頻點(diǎn)數(shù)據(jù)，總的實測數(shù)據(jù)不少于14×20 個。利用ArcGIS 中的Geostatistical Analyst 工具箱進(jìn)行數(shù)據(jù)檢驗，并對數(shù)據(jù)進(jìn)行對數(shù)變換，處理后的數(shù)據(jù)圖形接近直線（見附錄B 圖B1），可見數(shù)據(jù)服從正態(tài)分布，即數(shù)學(xué)期望和方差都存在，滿足二階平穩(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 條線路。各站點(diǎn)的經(jīng)緯度和線路的具體參數(shù)參見文獻(xiàn)［14］，根據(jù)數(shù)據(jù)構(gòu)建電網(wǎng)等效模型［15］從而計算GIC。

4.2 微元中點(diǎn)視電阻率計算

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

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

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

圖3 空間-頻率變異函數(shù)的實驗?zāi)Ｐ虵ig.3 Experimental model of spatial-frequency variation function

根據(jù)式（6）建立空間-頻率變異函數(shù)的實驗?zāi)Ｐ?，結(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ù)后，將微元中點(diǎn)作為預(yù)測點(diǎn)，計算權(quán)重系數(shù)，進(jìn)而得到待估點(diǎn)的視電阻率。

4.3 微元GMD地電場計算

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

計算得到肇慶站在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個站點(diǎn)進(jìn)行2種算法的準(zhǔn)確度比較。將計算所得肇慶站的300個GMD地電場數(shù)據(jù)、嶺澳站的300 個GIC 數(shù)據(jù)與實測值進(jìn)行對比分析，結(jié)果見表1（表中數(shù)據(jù)樣本數(shù)量為300）。

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

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

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

5 結(jié)論

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

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

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

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