基于視電阻率數據的電網地磁感應電流計算方法

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

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

0 引言

電網地磁擾動（GMD）災害防治有設備防治、調度防治2種方法。設備防治方法受限于資金投入，無法大量使用［1］；調度防治方法需利用電網運行調度及時補償地磁感應電流（GIC）引起的無功損耗［2］，因此研究計算精度高、速度快的GIC計算方法成為關鍵。

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

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

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

1 空間?頻率變異函數

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

1.1 假設條件

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

（1）在研究區域內，Z(xi，fi)的數學期望存在，且為常數，即：

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

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

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

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

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

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

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

1.2 實驗變異函數的估計

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

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

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

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

1.3 變異函數理論模型

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

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

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

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

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

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

2.1 無偏性

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

其中，λi為待定權重系數。

進而可得到：

2.2 估計誤差最小

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

其中，σ為標準差。

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

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

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

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

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

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

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

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

4 算例分析

圖1 計算流程Fig.1 Computation process

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

4.1 基礎數據

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

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

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

4.2 微元中點視電阻率計算

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

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

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

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

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

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

4.3 微元GMD地電場計算

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

計算得到肇慶站在22:48:00—22:52:00 時段內的GMD 地電場分量Ex、Ey，結果分別見圖4（a）、（b）。根據式（15）計算等效電壓源，并結合廣東地區的電網參數計算得到嶺澳站的GIC，結果見圖4（c）。

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

4.4 對比分析

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

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

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

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

5 結論

（1）本文首次將地質統計學中的變異函數引入GMD 地電場計算中，從線路角度出發，將研究區域內實測視電阻率作為原始數據，利用變異函數求解線路微元中點對應的視電阻率。結合地磁臺地磁場的測量數據計算每個微元的地電場，GIC 和GMD 地電場的計算結果與實際監測結果之間的擬合程度更高、更準確。

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

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

附錄見本刊網絡版（http：//www.epae.cn）。