基于广义帕累托分布的中纬度地区地磁场极值预测

周宁馨 刘青 查虹丽 马龙雄

摘 要：极端地磁活动会对电力系统的正常运行造成影响，研究地磁场量的极值水平可以为量化电网中的GIC受地磁暴影响程度提供理论支撑。选取兰州地磁台的观测数据为典型算例，利用基于广义帕累托分布（GPD）的超阈值模型（POT）对磁暴期间地磁场的水平分量及其分钟变化率进行拟合，并结合概率图（P-P图）和分位数图（Q-Q图）分析模型检验结果。利用轮廓似然估计方法得到50 a，100 a和200 a一遇的地磁场水平分量及其分钟变化率的重现水平，并对重现水平及其95%置信区间随着重现期增长的变化趋势进行讨论。通过对中纬度地区的5个地磁台重现值结果的综合分析，给出了100 a一遇的磁暴期间地磁场分钟变化率约为200～500 nT/min，200 a一遇的磁暴期间地磁场分钟变化率约为200～800 nT/min水平。模型估计的结果表明POT模型适用于地磁场量的极值分析，且地磁场变化率的重现水平随着重现期增加呈非线性增大趋势，给出的中纬度地区地磁场变化率极值范围，可为计算电网中地磁感应电流和分析磁暴灾害提供理论支撑。

关键词：地磁暴;广义帕累托分布;超阈值模型;预测

中图分类号：P 318

文献标志码：A

文章编号：1672-9315（2021）03-0524-07

DOI：10.13800/j.cnki.xakjdxxb.2021.0318

Abstract：An examination of the geomagnetic extreme value can provide theoretical support for quantifying the GIC in the power grid，and peak over threshold model（POT）based on the generalized Pareto distribution（GPD）has found wide applications.The observation data of the Lanzhou station was selected as an example，and the POT model was used to fit the data with the results checked.Also the profile likelihood estimation was used to obtain the return level under different return periods.Among five typical midlatitude observatories，horizontal geomagnetic field changes may reach 200～500 nT/min in one magnetic storm once every 100 years，and for return periods of 200 years the equivalent figure is 200～800 nT/min.The results show that the return level of the rate of change of the geomagnetic field increases nonlinearly with the increasing of the return period.Key words：geomagnetic storms;generalized Pareto distribution;peak over threshold model;prediction

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 地磁暴诱发的感应地电场

现有地磁暴感应地电场计算方法中应用最广泛的为平面波理论，即假设空间电流源为电流密度恒定，距离地面足够远的无限大面电流。若大地电导率均匀，则根据麦克斯韦方程和法拉第电磁感应定律，通过地磁台观测的地磁场水平分量估算出感应地电场为

其中bn=Bn-Bn-1为地磁场分量的一阶差分。通过式（4）可以看出影响磁暴感应地电场的关键因素是地磁场水平分量变化率。

以2006年12月的地磁暴事件为例，基于兰州地磁台（LZH）磁暴监测数据计算感应地电场。图1给出了磁暴期间Bx和dBx/dt随时间的变化曲线。设置均匀大地电导率为0.001 S/m，由此计算的东西方向感应地电场Ey如图2所示。

对比可知，相比较于地磁场分量Bx，地电场幅值出现最大值的时刻与地磁场变化率dBx/dt更相关，证明地磁场变化率是地磁活动影响电力系统过程中最主要的物理量，故文中选取地磁场水平分量BH和变化率dBH/dt进行概率模型的拟合。

2 地磁场的极值预测模型

2.1 极值理论与应用

鉴于地磁场的观测数据有独立的分钟值样本，文献[8]提出在地磁活动的统计分析中广义帕累托分布（GPD）最为适合。这种基于广义帕累托分布对超过某一充分大的阈值的所有观测数据进行拟合的模型称为超阈值（POT）模型。

假定地磁场水平分量及水平分量变化率的随机变量列x1，x2，…，xn服从GPD分布，其表达式为

式中 μ为位置参数;σ为尺度参数（σ>0）;ξ为形状参数。位置参数μ将选定为阈值u。随机变量列x1，x2，…，xn是否服从GPD分布还需要进行模型检验。

通过将地磁台的观测数据集与GPD函数进行拟合建立POT模型，可以评估观测数据集以外更为极端的地磁活动发生概率[17]。通常采用极大似然估计法确定模型中的参数，也可利用轮廓似然函数得到更为精确参数的估计范围[18]。

在对参数进行估计后，通过概率分布函数中分位数p的定义，可以得到某一确定重现期T下磁暴极值的估算公式

式中 N/Nμ为样本总数与样本中超出量个数之比;xp为重现期为T=1/（1-p）所对应重现水平。

利用P-P图（概率图）与Q-Q图（分位数图）进行模型检验，理论上当某种分布函数与样本实际分布越贴近时，P-P图和Q-Q图应该更近似为直线[19]。通过分析模型检验结果可以确定选取的概率分布函数是否符合实际数据的真实分布情况。

2.2 地磁数据资料

选取我国观测状态正常、保存数据记录年份较长的兰州地磁台（LZH）进行分析[20-23]。结合我国超、特高压输电网络主要处于中纬度范围的地理特点，同样选取国际上其他4个中纬地磁台数据进行GPD拟合，并将重现值估计结果对比分析。地磁台地理坐标和数据记录时间见表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]，与之前得到的区间相比范围稍有缩小。ξ所在置信区间为负意味着对应一个有界分布，这与地磁场水平分量具有上限的实际观测情况相符合。

图5为利用GPD分布拟合地磁场水平分量的POT模型诊断图。图5（a）和图5（b）分别为P-P图和Q-Q图，可以明显看出样本数据点均落在一条直线。图5（c）为重现水平估计值随着重现期变化的趋势及其95%置信区间，由于ξ的估计值和考虑置信区间后均为负值，相应分布具有有限上界，故重现水平渐进地趋于某个有限值。最后，图5（d）的概率密度曲线也与样本数据直方图相吻合。诊断图表明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的擬合。

3 重现水平估计

在建立地磁场水平分量BH和地磁场水平分量变化率dBH/dt的POT模型后，根据式（6）计算出不同重现期下的轮廓对数似然曲线，如图7和图8所示，包括50 a一遇、100 a一遇和200 a一遇极值的重现水平估计。图中横坐标为重现水平范围，纵坐标为轮廓对数似然函数，其极大值对应横坐标即为重现水平估计值，用于参考的水平直线指出了该重现水平的95%置信区间。

可以看出随着重现期的增加，地磁场水平分量BH和地磁场水平分量变化率dBH/dt的轮廓对数似然曲线均表现出更强的不对称性。由于重现期越长，观测数据所能提供的信息就越少，得到重现水平估计值的95%置信区间就越宽。

表2给出在不同重现期下，地磁场水平分量分钟变化率dBH/dt的估计上限和以最大观测值为基准的放大倍数，表中第1行为20 a观测数据中的最大值。可以看出随着重现期的增长，dBH/dt呈现出较为缓和的增长趋势，而非简单的线性关系。相比较于采用统一倍数标准对观测历史最大值进行放大[2]，POT模型可以为极端磁暴活动提供更为准确的预测极值参考。

考虑到经纬度的差异，利用GPD分布对全球范围内中纬度地区的5个地磁站台观测数据进行拟合，给出了最大观测值与所建立的POT模型中100 a一遇和200 a一遇估计结果的对比。图9和图10分别为5个不同站台的地磁场水平分量BH和地磁场水平分量变化率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）地磁场变化率的重现水平随着重现期增加呈非线性增大趋势，不能简单地采用统一放大倍数构建极端磁暴场景。

参考文献（References）：

[1] 刘连光，刘春明，张冰，等.中国广东电网几次强磁暴影响事件[J].地球物理学报，2008，51（4）：976-981.LIU Lianguang，LIU Chunming，ZHANG Bing，et al.Strong magnetic storms influence on Chinas Guangdong power grid[J].Chinese Journal of Geophysics，2008，51（4）：976-981.[2]刘连光，刘春明，张冰.磁暴对我国特高压电网的影响研究[J].电网技术，2009，33（11）：1-5.LIU Lianguang，LIU Chunming，ZHANG Bing.Effects of geomagnetic storm on UHV power grids in China[J].Power System Technology，2009，33（11）：1-5.

[3]刘连光.大规模电网应对空间灾害天气的问题[J].电网技术，2010，34（6）：1-5.LIU Lianguang.Scientific issues on how to cope with damage in large-scale power grid caused by disastrous space weather[J].Power System Technology，2010，34（6）：1-5.[4]李海明，陶勇，张俊双，等.基于1989年3月地磁暴的蒙东电网事故风险评估[J].电网技术，2020：1-8.LI Hai，TAO Yong，ZHANG Junshuang，et al.Risk assessment of Mengdong power grid accident based on geomagnetic storm in march 1989[J].Power System Technology，2020：1-8.[5]王泽忠，董博，刘春明，等.华北地区大地电性结构三维建模及磁暴感应地电场有限元计算[J].电工技术学报，2015，30（3）：61-66.WANG Zezhong，DONG Bo，LIU Chunming，et al.Three-dimensional earth conductivity structure modelling in North China and calculation of geoelectromagnetic fields during geomagnetic disturbances based on finite element method[J].Transactions of China Electrotechnical Society，2015，30（3）：61-66.[6]董博.磁暴情况下复杂地质结构三维地电流场有限元分析[D].北京：华北电力大学，2015.DONG Bo.Analysis on three-dimensional geoelectric field induced in complicated conductivity structure during geomagnetic storm by finite element method[D].Beijing：North China Electric Power University，2015.

[7]林晨翔.大地电导率横向变化对磁暴感应地电场的影响研究[D].北京：华北电力大学，2016.LIN Chenxiang.Effect of conductivity lateral change on induced geo-electric field during magnetic storms[D].Beijing：North China Electric Power University，2016[8]PULKKINEN A，BERNABEU E，EICHNER J，et al.Generation of 100 year geomagnetically induced current scenarios[J].Space Weather，2012，10（4）：1-19.[9]王建城，苏盛，盛小勇，等.输电线路多年一遇极值覆冰估计方法适用性分析[J].电网技术，2015，39（9）：2614-2620.WANG Jiancheng，SU Sheng，SHENG Xiaoyong，et al.Comparative study of applicability of methods for estimating transmission line icing return period based on various extreme value distributions[J].Power System Technology，2015，39（9）：2614-2620.[10]BEGGAN C D，BEAMISH D，RICHARDS A，et al.Prediction of extreme geomagnetically induced currents in the UK high-voltage network[J].Space Weather，2013，11（7）：407-419.[11]THOMSON A W P，DAWSON E B，REAY S J.Quantifying extreme behavior in geomagnetic activity[J].Space Weather，2011，9（10）：1-12.[12]KAPPERNMAN J G，ALBERTSON V D.Bracing for the geomagnetic storms[J].IEEE Spectrum，1990，27（3）：27-33.[13]徐文耀.地球电磁现象物理学[M].合肥：中国科学技术大学出版社，2009.[14]赵旭东，何宇飞，陈俊，等.基于地磁台站数据对磁暴期间环电流和场向电流的分布特征研究[J].地球物理学报，2019，62（9）：3209-3222.ZHAO Xudong，HE Yufei，CHEN Jun，et al.The distribution of ring current and field-aligned current during storms based on ground observatory data[J].Chinese Journal of Geophysics，2019，62（9）：3209-3222.[15]ESPINOSA K V，PADILHA A L，ALVES L R.Effects of ionospheric conductivity and ground conductance on geomagnetically induced currents during geomagnetic storms：case studies at low-latitude and equatorial regions[J].Space Weather，2019，17（2）：252-268.[16]吳伟丽，赵普志，蔡其东.基于Pareto分布的磁暴感应地电场极值测度模型[J].空间科学学报，2019，39（1）：28-35.WU Weili，ZHAO Puzhi，CAI Qidong.Extreme value model of mid-and low-latitude geoelectric field due to magnetic storms[J].China Journal Space Science，2019，39（1）：28-35.[17]史道济.实用极值统计方法[M].天津：天津科学技术出出版社，2006.

[18]COLES S.An introduction to statistical modelling of extreme values[J].Springer，London，2004.[19]钱小仕，王福昌，盛书中.基于广义帕累托分布的地震震级分布尾部特征分析[J].地震学报，2013，35（3）：341-350，450.QIAN Xiaoshi，WANG Fuchang，SHENG Shuzhong.Characterization of tail distribution of earthquake magnitudes via generalized Pareto distribution[J].Acta Seismologica Sinica，2013，35（3）：341-350，450.[20]李细顺，高登平，王利兵，等.1997—2016年中国大陆地区地磁台站观测环境变化定量分析[J].震灾防御技术，2019，14（3）：686-697.LI Xishun，GAO Dengping，WANG Libing.Quantitative analysis of the observed environmental changes of geomagnetic stations in mainland China from 1997 to 2016[J].Technology for Earthquake Disaster Prevention，2019，14（3）：686-697.[21]周锦屏，安郁秀.武汉地磁台环境磁污染的观测与初步分析[C]//中国科学院地球物理研究所40周年所庆论文集：中国岩石力学与工程学会，1990：172-174.ZHOU Jinping，AN Yuxiu.Observation and preliminary analysis of environmental magnetic pollution of Wuhan geomagnetic station[C]//Chinese Society of Rock Mechanics and Engineering，1990：172-174.[22]王晓祜，吕智，周江林，等.北京地磁台的磁环境监测[J].地震地磁观测与研究，2007（3）：52-61.WANG Xiaohu，LV Zhi，ZHOU Jianglin，et al.The geomagnetic gradients monitoring of the Beijing geomagnetic observatory[J].Seismological and Geomagnetic Observation and Reaserch，2007（3）：52-61.[23]韩英，杨莉，闫万生，等.兰州地磁台观测仪器工作状态分析[J].甘肃科技，2017，33（19）：45-47.HAN Ying，YANG Li，YAN Wansheng，et al.Analysis of working status of observation instruments of Lanzhou Geomagnetic Station[J].Gansu Science and Technology，2017，33（19）：45-47.[24]卓志，王伟哲.巨灾风险厚尾分布：POT模型及其应用[J].保险研究，2011（8）：13-19.ZHUO Zhi，Wang Weizhe.Thick-tailed distribution of catastrophe risk：POT model and application[J].Insurance Studies，2011（8）：13-19.[25]PULKKINEN A，PIRJOLA R，VILJANEN A.Statistics of extreme geomagnetically induced current events[J].Space Weather，2008，6（7）：1-10.