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

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

0 引 言

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

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

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

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

國內(nèi)利用極值理論對極端天氣進行重現(xiàn)水平的分析應(yīng)用廣泛，但在地磁暴領(lǐng)域的研究目前相對較少，吳偉麗等人利用磁場數(shù)據(jù)和復(fù)鏡像法構(gòu)造地電場統(tǒng)計樣本，通過建立地電場的GPD模型估算出喀什地區(qū)50 a一遇和100 a一遇的感應(yīng)地電場強度[16]。

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

1 地磁暴誘發(fā)的感應(yīng)地電場

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

(1)

(2)

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

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

(3)

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

E(TN)=

(4)

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

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

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

圖2 磁暴期間感應(yīng)地電場

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

2 地磁場的極值預(yù)測模型

2.1 極值理論與應(yīng)用

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

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

(5)

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

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

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

(6)

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

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

2.2 地磁數(shù)據(jù)資料

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

表1 臺站列表

2.3 模型擬合與檢驗

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

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

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

圖4 形狀參數(shù)ξ的輪廓似然估計

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

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

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

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

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

3 重現(xiàn)水平估計

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

圖7 地磁場水平分量BH重現(xiàn)水平的輪廓似然估計

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

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

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

表2 T年重現(xiàn)期的dBH/dt重現(xiàn)水平上限值

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

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

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

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

4 結(jié) 論

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

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

3)地磁場變化率的重現(xiàn)水平隨著重現(xiàn)期增加呈非線性增大趨勢，不能簡單地采用統(tǒng)一放大倍數(shù)構(gòu)建極端磁暴場景。