基于廣義帕累托分布的中緯度地區地磁場極值預測

周寧馨，劉 青，查虹麗，馬龍雄

(西安科技大學 電氣與控制工程學院，陜西 西安 710054)

0 引 言

極端地磁暴是一種客觀存在的自然災害。磁暴期間的感應地電場在電網中誘發的地磁感應電流(GIC)，對電網的威脅已引起普遍關注[1-3]。

目前國內關于感應地電場和地磁感應電流的計算多采用固定的地磁場變化率，或者針對某一次特定的地磁暴事件進行評估分析[4-7]。實際上，自然災害對實際工程造成的潛在風險進行評估時，需要考慮的并非已經觀測到的一般事件的概率，而是極端事件可能發生的微小概率，因此往往采用百年重現水平值作為評估的輸入條件[8-9]。

BEGGAN等人將地磁場極值估計的結果與英國大地模型和高壓輸電網絡相結合，預測了100 a一遇和200 a一遇的極端地磁暴場景產生的變壓器節點GIC[10]。這項研究采用的極值數據來源于THOMSON等人對歐洲地磁場重現水平的估計，他們利用極值理論建立地磁場水平分量的分鐘值及其變化率的廣義帕累托分布(generalized pareto distribution，GPD)模型，得到高緯度地區100 a一遇和200 a一遇地磁場可能出現的極值范圍[11]。

地磁暴作為一種全球性、幾乎同時發生變化的地磁場擾動，從赤道到極區均可觀測到顯著變化現象[12]。但不同緯度地區地磁場所受到磁暴影響的機理有所差異[13-15]，因此需要結合中國高壓輸電網絡主要處于中緯度范圍的地理特點，選取中緯度地區的地磁臺觀測數據進行分析。

國內利用極值理論對極端天氣進行重現水平的分析應用廣泛，但在地磁暴領域的研究目前相對較少，吳偉麗等人利用磁場數據和復鏡像法構造地電場統計樣本，通過建立地電場的GPD模型估算出喀什地區50 a一遇和100 a一遇的感應地電場強度[16]。

文中以蘭州地磁臺的觀測數據為基礎，利用基于廣義帕累托分布的超閾值模型(peak over threshold，POT)對地磁場的水平分量及其分鐘變化率進行擬合，通過極大似然估計方法和輪廓似然函數得到T年一遇的重現值結果，為施加地磁場激勵的幅值選取提供理論參考。

1 地磁暴誘發的感應地電場

現有地磁暴感應地電場計算方法中應用最廣泛的為平面波理論，即假設空間電流源為電流密度恒定，距離地面足夠遠的無限大面電流。若大地電導率均勻，則根據麥克斯韋方程和法拉第電磁感應定律，通過地磁臺觀測的地磁場水平分量估算出感應地電場為

(1)

(2)

式中Ex，Ey分別為南北及東西向感應地電場；μ0為真空磁導率；σ為大地電導率；gx(u)為南北方向磁場變化率；gy(u)為東西方向磁場變化率。

由于地磁臺數據為離散序列，在編程計算時需要將公式離散化。設地磁臺記錄水平分量B數據的采樣周期表示為Δ=Tn-Tn-1，在采樣間隔[Tn-1,Tn]內有g(t)=(Bn-Bn-1)/Δ。設t=TN，代入式(1)并進行離散化處理得

(3)

將g(t)=(Bn-Bn-1)/Δ代入上式得

E(TN)=

(4)

其中bn=Bn-Bn-1為地磁場分量的一階差分。通過式(4)可以看出影響磁暴感應地電場的關鍵因素是地磁場水平分量變化率。

以2006年12月的地磁暴事件為例，基于蘭州地磁臺(LZH)磁暴監測數據計算感應地電場。圖1給出了磁暴期間Bx和dBx/dt隨時間的變化曲線。設置均勻大地電導率為0.001 S/m，由此計算的東西方向感應地電場Ey如圖2所示。

圖1 磁暴期間地磁場水平分量及其變化率

圖2 磁暴期間感應地電場

對比可知，相比較于地磁場分量Bx，地電場幅值出現最大值的時刻與地磁場變化率dBx/dt更相關，證明地磁場變化率是地磁活動影響電力系統過程中最主要的物理量，故文中選取地磁場水平分量BH和變化率dBH/dt進行概率模型的擬合。

2 地磁場的極值預測模型

2.1 極值理論與應用

鑒于地磁場的觀測數據有獨立的分鐘值樣本，文獻[8]提出在地磁活動的統計分析中廣義帕累托分布(GPD)最為適合。這種基于廣義帕累托分布對超過某一充分大的閾值的所有觀測數據進行擬合的模型稱為超閾值(POT)模型。

假定地磁場水平分量及水平分量變化率的隨機變量列x1，x2，…，xn服從GPD分布，其表達式為

(5)

式中μ為位置參數；σ為尺度參數(σ>0)；ξ為形狀參數。位置參數μ將選定為閾值u。隨機變量列x1，x2，…，xn是否服從GPD分布還需要進行模型檢驗。

通過將地磁臺的觀測數據集與GPD函數進行擬合建立POT模型，可以評估觀測數據集以外更為極端的地磁活動發生概率[17]。通常采用極大似然估計法確定模型中的參數，也可利用輪廓似然函數得到更為精確參數的估計范圍[18]。

在對參數進行估計后，通過概率分布函數中分位數p的定義，可以得到某一確定重現期T下磁暴極值的估算公式

(6)

式中N/Nμ為樣本總數與樣本中超出量個數之比；xp為重現期為T=1/(1-p)所對應重現水平。

利用P-P圖(概率圖)與Q-Q圖(分位數圖)進行模型檢驗，理論上當某種分布函數與樣本實際分布越貼近時，P-P圖和Q-Q圖應該更近似為直線[19]。通過分析模型檢驗結果可以確定選取的概率分布函數是否符合實際數據的真實分布情況。

2.2 地磁數據資料

選取我國觀測狀態正常、保存數據記錄年份較長的蘭州地磁臺(LZH)進行分析[20-23]。結合我國超、特高壓輸電網絡主要處于中緯度范圍的地理特點，同樣選取國際上其他4個中緯地磁臺數據進行GPD擬合，并將重現值估計結果對比分析。地磁臺地理坐標和數據記錄時間見表1。

表1 臺站列表

2.3 模型擬合與檢驗

利用R語言進行數據處理，建立LZH地磁臺2001—2019年的地磁場水平分量分鐘值BH的POT模型。擬合過程中需要選取某一確定閾值，閾值太大將只有少數幾個超出量，導致估計量方差過大；而閾值太小則使得超出量與GPD分布出現較大偏差，估計量成為有偏估計[24]。

圖3給出了形狀參數和修正尺度參數隨閾值變化的趨勢，修正尺度參數σ*=σ-ξμ。若閾值適當，則超出量服從GPD分布，參數的估計值應該保持在變化較小的范圍內。從圖3可以看出在[31 050，31 055]的區間內，ξ和σ*的估計量都比較穩定。為了保證模型的準確性，應在參數估計量呈平穩趨勢的區間內盡量選取較大的閾值，因此選擇μ=31 055，此時閾值門檻的篩選比例為99.97%，與文獻[8]中確定的篩選標準相近，同時超閾值的樣本量不足1%，可以較好的滿足GPD尾部建模條件。

通過求解GPD分布的對數似然函數的最大值，得到尺寸參數σ極大似然估計為9.295 6，標準誤差1.239 0；形狀參數ξ極大似然估計為-0.291 9，標準誤差0.085 8。由估計值和標準誤差可知σ和ξ置信度近似為95%的置信區間分別為[6.87，11.72]和[-0.46，-0.12]。由于形狀參數ξ決定分布的類型，圖4利用ξ的輪廓對數似然函數得到更高精度的95%置信區間為[-0.44，-0.09]，與之前得到的區間相比范圍稍有縮小。ξ所在置信區間為負意味著對應一個有界分布，這與地磁場水平分量具有上限的實際觀測情況相符合。

圖4 形狀參數ξ的輪廓似然估計

圖5為利用GPD分布擬合地磁場水平分量的POT模型診斷圖。圖5(a)和圖5(b)分別為P-P圖和Q-Q圖，可以明顯看出樣本數據點均落在一條直線。圖5(c)為重現水平估計值隨著重現期變化的趨勢及其95%置信區間，由于ξ的估計值和考慮置信區間后均為負值，相應分布具有有限上界，故重現水平漸進地趨于某個有限值。最后，圖5(d)的概率密度曲線也與樣本數據直方圖相吻合。診斷圖表明POT模型適用于地磁數據的擬合，而且在進行模型數值外推時具有穩定表現。

圖5 地磁場水平分量BH的POT模型擬合診斷

在分析地磁場水平分量變化率時同樣利用式(5)擬合建立POT模型，數據來源于BH進行差分處理后的分鐘值變化率dBH/dt。結合MRL曲線和參數隨閾值變化情況，選取的閾值門檻為99.5%。對GPD分布中參數進行估計，得到尺寸參數σ極大似然估計為6.311 6，標準誤差2.350 3；形狀參數ξ極大似然估計為0.071 4，標準誤差0.305 7。通過計算得到形狀參數ξ的95%置信區間為[-0.53，0.67]，利用輪廓似然曲線得到更精確的置信區間為[-0.28，0.64]。

擬合診斷圖如圖6所示。圖6(a)和圖6(b)分別為P-P圖和Q-Q圖，可以看出樣本點幾乎均勻地落在參考直線附近，表明POT模型適用于對樣本數據的擬合。圖6(c)的重現水平曲線近似為線性，考慮到ξ的估計值和大部分置信區間都為正，相應分布的沒有上界，這表明極端地磁活動不被某些最大值所限制；同時服從Gumbel或Fraechet尾部分布，意味著其概率將隨著重現水平的增大而呈指數或多項式下降趨勢，在更長的重現期才可能會出現更高的極值。圖6(d)的概率密度曲線與樣本實際分布的直方圖較為吻合，綜合4個診斷圖認為POT模型同樣適用于對地磁場水平分量分鐘值變化率dBH/dt的擬合。

圖6 地磁場水平分量變化率dBH/dt的POT模型擬合診斷

3 重現水平估計

在建立地磁場水平分量BH和地磁場水平分量變化率dBH/dt的POT模型后，根據式(6)計算出不同重現期下的輪廓對數似然曲線，如圖7和圖8所示，包括50 a一遇、100 a一遇和200 a一遇極值的重現水平估計。圖中橫坐標為重現水平范圍，縱坐標為輪廓對數似然函數，其極大值對應橫坐標即為重現水平估計值，用于參考的水平直線指出了該重現水平的95%置信區間。

圖7 地磁場水平分量BH重現水平的輪廓似然估計

圖8 地磁場變化率dBH/dt重現水平的輪廓似然估計

可以看出隨著重現期的增加，地磁場水平分量BH和地磁場水平分量變化率dBH/dt的輪廓對數似然曲線均表現出更強的不對稱性。由于重現期越長，觀測數據所能提供的信息就越少，得到重現水平估計值的95%置信區間就越寬。

表2給出在不同重現期下，地磁場水平分量分鐘變化率dBH/dt的估計上限和以最大觀測值為基準的放大倍數，表中第1行為20 a觀測數據中的最大值。可以看出隨著重現期的增長，dBH/dt呈現出較為緩和的增長趨勢，而非簡單的線性關系。相比較于采用統一倍數標準對觀測歷史最大值進行放大[2]，POT模型可以為極端磁暴活動提供更為準確的預測極值參考。

表2 T年重現期的dBH/dt重現水平上限值

考慮到經緯度的差異，利用GPD分布對全球范圍內中緯度地區的5個地磁站臺觀測數據進行擬合，給出了最大觀測值與所建立的POT模型中100 a一遇和200 a一遇估計結果的對比。圖9和圖10分別為5個不同站臺的地磁場水平分量BH和地磁場水平分量變化率dBH/dt的極端磁暴場景估計值，括號內為站臺所在緯度。

圖9 中緯度地區地磁場水平分量BH估計值

圖10 中緯度地區地磁場水平分量變化率dBH/dt估計值

從圖10中可以看出，隨著重現期的增加，地磁場變化率dBH/dt的重現水平有較大幅度的增加。此外，緯度相近的地磁站臺SPT與SUA的BH和dBH/dt重現水平均有較大差異，表明dBH/dt不僅與所處緯度有關。綜合全球范圍中緯度各個地磁臺的極值預測情況，給出中緯度地區100 a一遇的磁暴期間地磁場分鐘變化率約為200～500 nT/min，200 a一遇的磁暴期間地磁場分鐘變化率約為200～800 nT/min。

4 結 論

1)基于廣義帕累托分布的POT模型適用于磁暴期間地磁場觀測數據的擬合，可以通過該模型估計地磁場的極值水平。

2)POT模型中形狀參數ξ最為重要，其取值范圍影響概率分布特征及重現水平曲線的變化趨勢。

3)地磁場變化率的重現水平隨著重現期增加呈非線性增大趨勢，不能簡單地采用統一放大倍數構建極端磁暴場景。