華北地區大地電性結構三維建模及磁暴感應地電場有限元計算

王澤忠 董 博 劉春明 劉麗平 劉連光

（1. 華北電力大學高電壓與電磁兼容北京市重點實驗室 北京 102206 2. 華北電力大學電氣與電子工程學院 北京 102206 3. 中國電力科學研究院 北京 100192）

1 引言

太陽風暴攜帶的帶電粒子會導致地球磁場發生全球性劇烈變化，并在具有導電性的地質結構中感應出地電流。地電流在大地導體中流通會使地面不同位置出現電位差。在地面電位差的作用下，地磁感應電流（GIC）通過交流電網的中性點、油氣管網的管道壁及高鐵系統的鐵軌等接地點與大地構成流通回路[1]，引起變壓器直流偏磁、管道外壁腐蝕、破壞陰極保護系統和鐵路信號故障等一系列安全問題。而地面電位差是由地電場引起，因此準確求解地電場對于磁暴侵害技術系統的防護和治理至關重要[2]。

確定地電場的分布不僅要了解地磁場的時變情況，還要考慮大地電性結構的實際分布[2,3]。我國所處的中低緯地區地磁場變化頻率一般位于 0.01～0.000 1Hz之間[4,5]。若大地電導率只隨深度分層變化，可用地面波阻抗表征電性結構，與地磁場相乘求得地電場。這種方法稱作平面波法[6]。倘若輸電線路、油氣管網或鐵路網沿海岸線架設，或途經內陸不同地質結構交界面、大地斷層等電導率突變顯著的地區，則需要考慮海水與陸地之間、不同地質結構之間電導率的橫向變化。進一步地，對于我國特高壓長距離輸電線路，其空間尺度達到幾百公里乃至上千公里[7,8]，所覆蓋的區域包含多種地質構造，大地電導率是空間三維分布的。文獻[9,10]建立了包含海陸分界面的二維大地電導率模型，但沒有采用電導率的實測數據；文獻[11-13]是將北美、歐洲的大范圍區域分解成多個獨立的一維結構，在每個一維結構內用平面波法求得地電場，再疊加起來獲得全區域的地電場分布。這種方法沒有考慮電導率的橫向差異對地電場的影響。

本文采用我國大地電磁測深數據構建華北地區大范圍三維電導率模型。以磁暴時地磁臺測得的地磁場變化數據作為邊界條件，應用有限元法求解三維感應地電場，給出了不同電導率結構影響地電流分布的規律。應用地面電場的計算結果，考察了規劃中的“三華”特高壓同步電網“東縱”和“北橫”部分線路中選定變電站間的電位差分布情況，為進一步制定特高壓電網防范磁暴侵害的技術方案、優化電網設計及臺站選址提供依據。

2 基本方程與建模方法

2.1 基本方程

由于地磁場變化的頻率很低，與傳導電流相比，位移電流可以忽略。研究感應地電場時假定大地導體為各向同性結構，磁導率為μ0。采用A-φ位表示的伽遼金加權余量方程為[14]

式中，V為整個求解域；S為場域邊界；Vc為導體區；標量權函數w等于基函數，矢量權函數W代表w與直角坐標系的單位矢量ex、ey和ez相乘得到的3個矢量權函數[15]；K定義為

其中，en為邊界面的外法向。則K的大小等于邊界處磁場強度Bμ0的切向分量。假設各單元上的形狀函數為N（ii=1,2,…,np，i為np單元中的節點總數），則A和φ在各單元上的插值函數為

將形狀函數在全部單元上疊加得到基函數序列，與節點上的A和φ組成近似函數，代入式（1）和式（2），離散成代數方程組，求解全部節點上的A和φ，進一步可求得磁感應強度、電場強度和電流密度。

2.2 建模方法

圖1 華北地區三維電導率模型范圍與“三華”特高壓電網部分規劃線路示意圖Fig.1 Sketch of three-dimensional conductivity structure of North China and two selected UHV transmission lines

根據大地電磁測深法觀測和解釋的華北地區地質構造成果[16-19]，不考慮巖石磁性和地面系統中導磁材料的影響，建立三維電導率模型，如圖1所示。圖中虛線與坐標軸圍成的區域為建模區域。圖2給出了地面附近的電性分布圖（實測數據為電阻率）。通過與我國“三華”特高壓同步電網規劃圖對比發現，“東縱”線路的錫盟至徐州段、“北橫”線路的蒙西至濟南段位于模型覆蓋的范圍內，兩條線路的走向及擬建變電站的大致選址位置也在圖1中標出。

圖2 模型地面附近電導率分布圖（以電阻率形式給出）Fig.2 The conductivity near the Earth’s surface at North China shown as resistivity

建模時忽略地球曲率影響，選取直角坐標系，模型的南北范圍為888km，東西范圍約為962.5km。考慮到趨膚效應，取模型下邊界為地下300km處，認為變化的電磁場在此處衰減至零。模型上邊界選為地面，邊界條件為模型的四周邊界選取的離關心區域足夠遠，感應電流平行或垂直于邊界面。前者對應磁場垂直邊界面條件，后者對應磁場平行邊界面條件。模型內部場的分布可通過式（1）與式（2）求得。

3 計算數據選取

圖3 2004年11月7日地磁場水平分量的記錄結果Fig.3 The horizontal components of geomagnetic field observed from one observatory in Nov.7th,2004

磁暴時地磁場的變化情況可由分布于各地的地磁臺站監測數據獲得。本文選取2004年11月7日發生磁暴時，離建模區域最近的地磁臺記錄的地磁場水平分量結果，如圖 3所示，數據采樣間隔為 1s。用文獻[20]的方法估算出，在11月7日的10：51：40～10：54：00之間的變化地磁場可能產生較大地電場。因此計算時選取這 140s的數據，相應的地磁場南北分量Bx、東西分量By相對于起始時刻的變化情況以內插圖的形式繪于圖中。

由于Bx和By對應的K的方向不同，四周邊界條件也不同。因此在同一個模型上需要施加各自的邊界條件分別計算，最后將二者結果疊加得到總的結果。

電磁場的唯一性定理表明[21,22]，場域內的場分布由區域內的基本方程和區域邊界上的邊界條件唯一確定。式（1）和式（2）表明，場域內的方程為渦流場方程，其中的電導率是通過實測獲得的。模型的邊界條件中，地面處的邊界條件是實際記錄的地磁場變化情況；底部邊界是變化電磁場所能達到的最遠邊界，考慮趨膚效應后場在此處衰減為零；四周邊界條件與電流和場在邊界的行為吻合。計算采用的有限元法為成熟的數值方法，采用六面體八節點單元，網格剖分足夠密。以上這些可以保證求解結果的可靠性[23]。

4 結果與分析

選取圖1中所示的兩條線路走向為電場的積分路徑，選徐州站的電位為參考電位，各變電站相對于參考點的電位差分布如圖4所示。曲線標號由小到大代表沿積分路徑上各站距離參考站的由近到遠。

圖4 兩條特高壓線路上各站相對于徐州站的電位差變化情況Fig.4 The voltage differences of substations at two UHV transmission lines refer to Xuzhou substation

由圖4可見，電位差最大值出現在距參考點約950km的錫盟站，大小約為220V。最大值出現的時刻為第67s。而距參考點最遠的蒙西站（約1 220km）并未因距離最遠而出現最大電位差。造成這種現象的一個原因是由于所選時段地磁場東西分量比南北分量變化更劇烈（見圖3），進而對南北走向的輸電線路影響更嚴重。

進一步考察晉北站的電位分布情況。晉北站距參考點約770km，與蒙西站相距450km，然而對比圖4中的5號和7號曲線可見，兩站的電位差分布大致相同。圖2顯示蒙西站所處位置的大地電導率較大，而晉北站所在地區的大地電導率較小，因此大地電性結構的差異也是影響變電站電位分布的一個原因。

圖 5給出了電位差達到最大時（第 67s）及其前后一些時刻地面電場的分布情況。在每個時刻圖中標出了電場最大值和最小值所在的位置以及大小。

圖5 電位差最大時刻及前后時刻華北地區地面電場分布情況Fig.5 The geoelectric field vectors at the moment and around the moment of maximum voltage difference

由圖5可見，磁暴時的地面電場在不同時刻的大小和方向都是變化的，且同一時刻的最大值和最小值差別很大。地面電場的最大值為0.87V/km，出現在E112°、N39°附近。將各時刻最大值的位置與圖2對比可見，在電導率比較小，并且周邊電導率與之差異較大的地區容易產生較大的地電場。這是因為感應電流在不同電導率結構之間流通時，為了保證電流連續性，電場強度在不同電導率交界面處出現突變。而電導率的差異越大，電場強度突變程度也越大。可見電導率的橫向差異會使地面電場的分布呈現空間不均勻性。

圖6給出了第67s時，地下不同深度的電流密度分布情況。從圖中可見，由于趨膚效應的存在，感應電流在模型上部較大，下部較小。地下200km處的最大感應電流大約是地下0.3km處最大感應電流的 1.1%。因此模型下邊界選為 300km可以確保變化電磁場在此處衰減至零，從而保證計算結果的準確可靠。如果地磁場的變化頻率更低，模型下邊界要選取得更深一些。

5 結論

（1）首次采用三維渦流場模型分析了復雜分布的大地電性結構對磁暴時地面電場的影響。計算結果表明，電導率的不連續使得電場強度出現突變，導致地面電場呈現空間的不均勻分布，并進一步影響位于不同地質結構的變電站的電位差分布。

（2）透入深度決定了模型邊界范圍的選取。磁暴時地磁場的變化頻率很低，電磁場會透入到地下幾百公里。從防范磁暴災害的角度出發，在這個深度范圍內的大地電性構造都要加以考慮。

圖6 電位差最大時刻不同深度處感應電流分布情況Fig.6 The induced currents at the moment of maximum voltage difference

（3）建模的電導率數據來自于大地電磁測深的探測結果，計算的地磁場數據來自于地磁臺的實際記錄，計算結果真實可靠，相應結論為準確評估特高壓電網受磁暴侵害的影響奠定了基礎，為電網規劃建設及變電站選址提供了參考依據。

[1] Pirjola R. Geomagnetically induced currents during magnetic storms[J]. IEEE Transactions on Plasma Science,2000,28(6): 1867-1873.

[2] 劉連光,吳偉麗. 磁暴影響電力系統安全風險評估思路與理論框架[J]. 中國電機工程學報,2014,34(10): 1583-1591.

Liu Lianguang,Wu Weili. Security risk assessment ideas and theoretical framework for power system considering geomagnetic storm[J]. Proceedings of the CSEE,2014,34(10): 1583-1591.

[3] Zheng K,Trichtchenko L,Pirjola R,et al. Effects of geophysical parameters on GIC Illustrated by benchmark network modeling[J]. IEEE Transactions on Power Delivery,2013,28(2): 1183-1191.

[4] Kappenman J. Storm sudden commencement events and the associated geomagnetically induced current risks to ground-based systems at low-latitude and midlatitude locations[J]. Space Weather,2003,1(3):1016-1031.

[5] 劉春明. 中低緯電網地磁感應電流及其評估方法研究[D]. 北京: 華北電力大學,2009.

[6] Marti L,Rezaei Zare A,Boteler D. Calculation of induced electric field during a geomagnetic storm using recursive convolution[J]. IEEE Transactions on Power Delivery,2014,29(2): 802-807.

[7] Liu C M,Liu L G,Pirjola R. Geomagnetically induced currents in the high-voltage power grid in china[J].IEEE Transactions on Power Delivery,2009,24(4):2368-2374.

[8] Liu C M,Li Y L,Pirjola R. Observations and modeling of GIC in the Chinese large-scale highvoltage power networks[J]. Journal of Space weather and Space Climate,2014,4(A03): 1-7.

[9] Gilbert J L. Modelling the effect of the ocean-land interface on induced electric fields during geomagnetic storms[J]. Space Weather,2005,(S04A03): 1-9.

[10] Pirjola R. Practical model applicable to investigating the coast effect on the geoelectric field in connection with studies of geomagnetically induced currents[J].Advances in Applied Physics,2013,1(1): 9-28.

[11] Wei L H,Homeier N,Gannon J L. Surface electric fields for north America during historical geomagnetic storms[J]. Space Weather,2013,11(8): 451-462.

[12] Marti L,Yiu C,Rezaei Zare A,et al. Simulation of geomagnetically induced currents with piecewise layered-earth models[J]. IEEE Transactions on Power Delivery,2014,29(4): 1886-1893.

[13] Viljanen A,Pirjola R,Wik M,et al. Continental scale modelling of geomagnetically induced currents[J].Journal of Space weather and Space Climate,2012,2(A17): 1-11.

[14] 周克定. 工程電磁場數值計算理論方法及應用[M].北京: 高等教育出版社,1994.

[15] Bíró O,Preis K. On the use of magnetic vector potential in the finite element analysis of threedimensional eddy currents[J]. IEEE Transactions on Magnetics,1989,25(4): 3145-3159.

[16] 董樹文,李廷棟,陳宣華,等. 我國深部探測技術與實驗研究進展綜述[J]. 地球物理學報,2012,55(12): 3884-3901.

Dong Shuwen,Li Tingdong,Chen Xuanhua,et al.Progress of deep exploration in mainland China: A review[J]. Chinese Journal of Geophysics,2012,55(12): 3884-3901.

[17] 董樹文,李廷棟. SinoProbe—中國深部探測實驗[J].地質學報,2009,83(7): 895-909.

Dong Shuwen,Li Tingdong. SinoProbe: the exploration of the deep interior beneath the Chinese continent[J]. Acta Geologica Sinica,2009,83(7): 895-909.

[18] 魏文博,葉高峰,金勝,等. 華北地區東部巖石圈到電性結構研究—減薄的華北巖石圈特點[J]. 地學前緣,2008,15(4): 204-216.

Wei wenbo,Ye Gaofeng,Jin Sheng,et al. Geoelectric structure of lithosphere beneath eastern North China:features of a thinned lithosphere from magnetotelluric soundings[J]. Earth Science Frontiers,2008,15(4):204-216.

[19] 葉高峰,魏文博,鄧明,等. 青藏及華北陣列式區域大地電磁場標準觀測網建設方法與實驗[J]. 地質學報,2010,84(6): 801-807.

Ye Gaofeng,Wei wenbo,Deng Ming,et al. Construction methods and experiments of magnetotelluric standard network at the Tibetan plateau and north China[J]. Acta Geologica Sinica,2010,84(6): 801-807.

[20] 余華兵. 針對電網地磁感應電流的地磁場變化特征提取[D]. 北京: 華北電力大學,2014.

[21] 盛劍霓等. 工程電磁場數值分析[M]. 西安: 西安交通大學出版社,1991.

[22] 楊儒貴. 電磁定理和原理及其應用[M]. 成都: 西南交通大學出版社,2002.

[23] Dong B,Danskin D,Pirjola R,et al. Evaluating the applicability of the finite element method for modelling of geoelectric fields[J]. Annales Geophysicae,2013,31(10): 1689-1698.