ADTD数据计算雷暴日与目测雷暴日的均一性分析

朱华亮，华连生，温华洋，庞文静，梁 丽，孔芹芹

(1.安徽省气象信息中心，安徽 合肥 230031；2.中国气象局气象探测中心，北京 100081)

引 言

雷电灾害是“联合国国际减灾十年”公布的影响人类活动的严重灾害之一[1]。在任何给定时刻，世界上都有1800多场雷雨正在发生，每秒大约有100次雷击[2]。据统计，中国每年因雷电造成的人员伤亡达上千人，财产损失约70～100亿元[3-5]。因此，加强雷电现象的监测与预警，开展有效的雷电防护工作尤为重要。我国气象站从建国初期就开始对雷电现象进行人工观测记录，根据《地面气象观测规范》规定[6]，在1日内(20:00至次日20:00，北京时，下同)只要闻雷声1次或多次，就记录为1个雷暴日，形成的雷暴日数是目前描述雷电活动的唯一长期记录。随着科学技术的发展，人工观测项目正在逐步被自动观测系统取代。2013年地面气象观测业务调整后，取消雷暴日的人工观测，采用ADTD(advanced TOA and direction system)闪电定位系统监测雷电现象的活动[7]。然而，雷暴日是反映当地雷电活动频繁程度及其强度的物理参数，在防雷工程设计、雷击风险评估和雷电科学研究中应用广泛[8-9]，短时间内还不可被取代，所以，利用闪电定位系统监测数据来延续雷暴日数据的研究更为迫切。

关于闪电定位资料能否替代人工观测的雷暴日，很多研究从两种资料的差异、原因及关系等方面进行了探讨，发现两种观测数据的分布趋势大体相同，具有较好的一致性，且呈现不规则的片状函数关系，认为闪电定位数据替代人工观测雷暴日具有较好的合理性和可行性[10-16]，并可通过直接替代法、地闪密度法或二元法等将闪电定位数据转换为年雷暴日数[17-19]。然而，大多数研究侧重于闪电定位数据与人工观测雷暴日之间的差异性分析及转换关系，而对于通过闪电定位数据计算的年雷暴日数与人工观测年雷暴日数之间的均一性却鲜有探讨。均一的长序列气象资料是气候变化研究的基础[20-21]，为此，本文对全国范围内的闪电定位数据与人工观测雷暴日之间的关系进行研究，形成ADTD监测数据与人工观测年雷暴日数之间的关系式，旨在探讨通过ADTD监测数据计算的年雷暴日数与人工观测年雷暴日数之间的均一性，以期为全国各地区的年雷暴日数资料序列延长、雷电防护工程设计、雷暴气候变化分析以及雷电灾害风险评估提供参考。

1 资料与方法

1.1 资 料

全国共有843个国家基准气候站和国家基本气象站，各气象站自1961年以来有完整的人工观测雷暴日记录，至今已经积累了约60 a的雷暴日数据，这些数据为各地的雷电灾害风险评估、防雷工程设计和雷电科学研究等发挥了不可或缺的作用。2013年全国地面气象观测业务改革后取消了雷暴日的人工观测。因此，选取全国843个国家基准气候站和国家基本气象站1961—2013年人工观测雷暴日资料作为研究对象，其中631个国家基本气象站的资料用于最优匹配半径的确定和年雷暴日数转换公式的构建，212个国家基准气候站用于转换公式的效果评估。人工观测雷暴日资料均经过“台站、省级、国家级”的三级质量控制，资料的完整性为100%。

ADTD闪电定位系统是由中国科学研究院空间科学与应用研究所研制的第2代地闪定位系统，该系统采用磁向和时差联合法进行闪电探测。2010年，全国已建成近300套ADTD闪电探测子站，初步形成全国闪电监测网，并投入业务应用，截止2018年底，全国已布设400余套ADTD闪电探测子站，对全国范围内的闪电探测能力进一步提升，其中中部和东部地区ADTD站点以150 km基线组网，西部地区ADTD站点以300 km基线组网，具体分布如图1。因此，选取ADTD系统监测的全国范围内2010—2018年闪电定位数据作为研究对象，其中2010—2013年ADTD闪电定位数据用于年雷暴日数转换公式的构建和验证，2014—2018年ADTD闪电定位数据用于年雷暴日数序列的延长。为保证ADTD监测数据的可靠性，需对其进行质量控制。有研究发现闪电定位系统会将一些电流值较小的云闪记录为地闪，称为“小幅值地闪”，本文采用2005年IEEE所推荐的界定值，剔除电流幅值小于2 kA的数据[22-24]。另据研究发现[25]，闪电定位数据中电流强度小于10 kA的正地闪也可能来源于云闪的干扰，对这一部分数据同样进行剔除。此外，ADTD数据中有一些闪电数据电流幅值大于200 kA，而大于200 kA的闪电一般认为误差较大[26]，因此对这一部分数据也进行剔除。

图1 全国ADTD闪电定位系统站点分布Fig.1 Distribution of stations for ADTD lightning location system in China

1.2 方 法

1.2.1 最优匹配半径确定方法

闪电是云与云之间、云与地之间或者云体内各部位之间的强烈放电现象，发生闪电时闪电通道的温度非常高，周围的空气被加热迅速膨胀，就会产生雷声，亦即闪电和雷声同时发生[27]。由于人的感官观测范围有限，通常只能听到其周围15 km范围内的雷声，较好的情况下能达到20 km，而ADTD监测闪电的范围远大于人耳听到雷声的范围，为得到ADTD数据计算雷暴日和人工观测雷暴日相匹配的范围，需对不同半径范围内统计得到的雷暴日与人工观测雷暴日进行匹配分析。在处理ADTD数据时，规定1日内(20:00至次日20:00)ADTD监测的闪电次数大于或等于1次，就记为1个雷电日。以区域内气象站所处位置为圆心，在1～40 km半径范围内每间隔1 km，统计区域内气象站不同监测半径R下的年雷电日数，计算区域内不同半径R下年雷电日数与年雷暴日数的均方误差(MSER)，并取MSER最小的监测半径R为与人工观测雷暴日最匹配的范围，亦即最优匹配半径。具体计算公式如下：

(1)

式中：Tdij(d)为第i个气象站第j年人工观测的雷暴日数；TRij(d)为监测半径R下第i个气象站第j年统计的雷电日数；m为区域内的站点数；j为年份(j=1,2,3,4分别对应2010，2011，2012，2013年)；N*为正整数集。

1.2.2 惩罚最大F检验(PMFT)法

为探究基于ADTD数据延长的年雷暴日数序列相对于人工观测序列是否发生突变，采用PMFT方法对2013年前后的雷暴日数序列进行均一性检验。PMFT方法[28-29]经验性地考虑了时间序列的一阶滞后自相关，并嵌入多元线性回归算法，能够用于检验、订正包含一阶自回归误差数据序列的多个间断点(平均突变)。该方法进行序列突变点检验，可以不使用参考序列，且有效地避免了断点位置对置信度的影响。具体检验方法如下：

设随机变量εt独立且服从均值为0，方差为σ2的高斯分布，对于存在线性趋势β的年雷暴日数序列{Xt}，要检验序列在2013年是否存在突变点，作如下假设：

原假设(H0)：

H0:Xt=μ+βt+εtt=1961,1962,…,2018

(2)

备择假设(Ha):

(3)

式中：t表示年份；μ、μ1、μ2为趋势线的截距，其中μ1≠μ2。当Ha为真时，则2013年就为突变点，用△=|μ1-μ2|来表示平均突变的大小，其检验统计量(PF)为：

PF=P(k)Fc(k)

(4)

式中：k=2013；P(k)是一个经验构造的惩罚因子，具体见参考文献[28]；统计量Fc(k)为服从自由度为1和55的F分布，其构造如下：

(5)

1.2.3 评估指标

(6)

2 基于ADTD数据的年雷暴日数拟合公式

2.1 雷暴日与雷电日的最优匹配半径

中国幅员辽阔，纬度跨度较广，距海远近差距较大，加之地形复杂，雷暴日的分布也不尽相同。为得到ADTD监测雷电日和人工观测雷暴日的最优匹配半径，以省为单位，将全国分为28个区域(不包含香港、澳门和台湾地区)，其中北京、天津和河北作为一个区(京津冀)，上海和江苏作为一个区(沪苏)。应用公式(1)得到各区域的最优匹配半径见表1，发现各地区对应的最优半径不尽相同，大部分省份的最优匹配半径为9～15 km，其中新疆、青海和西藏地区对应最优匹配半径较大，重庆地区对应的最优匹配半径最小。这是由于新疆、青海和西藏等地区人稀地广，多为高原和山脉，地形起伏较大，布设的ADTD闪电定位仪间距较大，形成的闪电定位网对该地区闪电的探测效率较低，导致较大范围内的雷电日数与人工观测雷暴日数一致，而重庆地区以山地为主，对雷声的传播具有较大的阻隔作用，导致观测员监听到的雷声范围较小，进而造成最优匹配半径较小。

表1 2010—2013年中国各地区ADTD监测雷电日和人工观测雷暴日的最优匹配半径Tab.1 Optimal matching radius between thunderstorm days calculated using data from ADTD and artificial thunderstorm days in different regions of China during 2010-2013

2.2 年雷暴日数拟合公式的建立

目前，关于闪电定位数据与人工观测雷暴日数据转换关系的研究中，有研究指出二元法对人工观测雷暴日的拟合效果最好[18-19]。为此，根据最优半径对应的年雷电日数T和年地闪密度Ng(最优半径内的年云地闪电次数与其对应面积之比)，建立年雷电日数T、年地闪密度Ng与年雷暴日数Td的二元线性回归方程，其中建模样本为2010—2013年的国家级基本气象站人工观测数据和ADTD监测数据。表2列出各地区回归方程的拟合结果，各地区的回归方程均通过显著性水平为0.05的F检验，回归系数均通过显著性水平为0.05的t检验。大部分地区的回归方程的拟合优度R2较高，而重庆、四川等地形较为复杂的地区，以及新疆、西藏等ADTD站点较为稀疏的地区，回归方程的拟合优度R2偏低。在均方根误差RMSE方面，绝大部分站点的RMSE小于10.0 d，而在ADTD站点稀疏和地形较为复杂的地区，RMSE偏大。总体表明，利用年雷电日数和年地闪密度建立的年雷暴日数拟合方程较为合理。

2.3 年雷暴日数拟合公式的效果评估

表2 2010—2013年中国各地区二元线性回归方程的拟合结果Tab.2 Fitting results of binary linear regression equations in different regions of China during 2010-2013

表3 2010—2013年中国各地区基于基准气候站观测资料的年雷暴日数拟合公式的效果检验Tab.3 The effect test of the fitting formula of annual thunderstorm days based on the observation data of the base climate stations in each region of China during 2010-2013

3 年雷暴日数延长序列的均一性检验

利用2014—2018年ADTD监测数据，通过各地区年雷暴日数拟合公式，对843个气象站2014—2018年的年雷暴日数进行计算，得到各气象站2014—2018年的年雷暴日数序列。采用PMFT法对843个气象站1961—2018年的年雷暴日数序列进行均一性检验，以2013年作为断点，判断2014—2018年的延长序列与1961—2013年的人工观测序列相比是否有明显的跳变。采用公式(4)计算了843个站点的检验统计量，结果显示843个站点中有776个站点的检验统计量小于显著性水平为0.05的临界值11.05，表明这些站点的年雷暴日数序列在2013年前后无明显跳变，占比达92.1%，仅有67个台站的延长序列相比人工观测序列有明显跳变，这些站点主要集中在新疆、青海、西藏等地区(图2)，其中新疆、青海和西藏地区分别有23个、7个和10个非均一站点(记为Ⅰ类非均一站点)，占该区域内站点数比分别达34.9%、20.0%、34.5%。这是由于新疆、青海和西藏等地区的ADTD站点布局较为稀疏，对闪电的探测效率较低，定位误差较大，加之复杂地形的影响，导致年雷暴日数的延长序列值相对人工观测序列值的偏差较大，进而造成了序列的非均一性。从图2还可以看出，剩余27个年雷暴日数序列未通过均一性检验的站点分布在云南、四川、东南沿海等地区，经分析发现主要有3种原因造成这部分站点的年雷暴日数序列非均一：一是人工观测序列在2003—2013年存在异常点所致(记为Ⅱ类非均一站点)，如广东省汕头气象站在2011年人工观测的年雷暴日数只有21.0 d，与前后年份相比，明显偏少，造成了序列的不均一，类似这种情况的站点共有19个；二是气象站周边的ADTD闪电定位仪的探测性能下降所致(记为Ⅲ类非均一站点)，如河北省蔚县气象站、重庆市酉阳气象站、四川省色达气象站和云南省澜沧气象站，其人工观测序列的年平均雷暴日数分别为45.1、50.1、77.9和101.6 d，而延长序列的年平均雷暴日数分别为25.1、25.8、44.6和66.0 d，延长序列相比人工观测序列的年雷暴日数明显减少，表明站点周边的ADTD闪电定位仪监测到的闪电回击次数较少，经分析发现这4个气象站观测场内或附近均布设有ADTD闪电定位仪且布设时间都较早，ADTD闪电定位仪的探测性能有所下降，导致在气象站周围监测到的闪电回击次数较少；三是延长序列中某年拟合的年雷暴日数异常偏高所致(记为Ⅳ类非均一站点)，如四川省绵阳气象站的人工观测年平均雷暴日数为30.4 d，2014—2018年计算的年雷暴日数分别为33.0、34.0、36.0、36.0和73.0 d，2018年的年雷暴日数相对前几年明显偏高，造成了序列的非均一性，类似这种情况的站点还有四川省温江气象站。此外，浙江省玉环气象站和云南省昭通气象站的年雷暴日数序列也未通过均一性检验，造成序列非均一的具体原因不明(记为Ⅴ类非均一站点)。检验结果总体表明，除新疆、青海和西藏等地区外，其余地区利用二元法延长的年雷暴日数序列与人工观测序列相比无明显跳变，均一性较好。

图2 1961—2018年中国年雷暴日数序列均一和非均一的站点分布Fig.2 Distribution of meteorological stations with homogeneity or heterogeneity of annual thunderstorm days sequence in China during 1961-2018

4 结论与讨论

(1)研究发现各地区ADTD监测雷电日和人工观测雷暴日的最优匹配半径不尽相同，大部分省份的最优匹配半径介于9～15 km之间，这与钟颖颖等[11]研究发现观测员通常只能听到8～15 km范围内的雷声的研究结论一致。新疆、青海和西藏等地区由于人稀地广，地形起伏较大，多为高原和山脉，布设的ADTD闪电定位仪间距较大，形成的闪电定位网对该地区闪电的探测效率较低，定位误差较大，导致较大范围内的雷电日数与人工观测雷暴日数一致，故新疆、青海和西藏等地区的雷电日和雷暴日最优匹配半径较大。

(2)利用2010—2013年的国家基本气象站人工观测数据和ADTD监测数据，建立的年雷暴日数拟合方程对大部分地区具有较高的拟合优度R2和较低的均方根误差RMSE，而新疆、青海和西藏等地区由于地形复杂、ADTD站点布局较为稀疏等原因，回归方程的拟合优度R2有所偏低，均方根误差RMSE有所偏大。经2010—2013年的独立样本检验发现，回归方程计算的年平均雷暴日数与人工观测年平均雷暴日数在空间上具有相同的分布规律，绝大部分地区的年雷暴日数计算值与人工观测值的平均绝对误差小于10.0 d，平均相对误差小于30.0%，但新疆、青海和西藏等地区的平均绝对误差和平均相对误差都较大。总体表明，除新疆、青海和西藏等地区外，基于年雷电日数、年地闪密度建立的二元回归方程对年雷暴日数的拟合效果较好，可用于雷暴日数序列的延长。

(3)利用PMFT方法对843个气象站1961—2018年的年雷暴日数序列进行均一性检验，结果显示843个站点中有776个站点的序列在2013年前后无明显跳变，占比达92.1%，有67个台站的延长序列相比人工观测序列有明显跳变，这些站点主要集中在新疆、青海、西藏、东南沿海等地，造成序列非均一的原因有高山阻隔、ADTD探测站网稀疏、ADTD探测仪探测性能下降等。总体表明，利用二元法延长的年雷暴日数序列与人工观测序列相比无明显跳变，均一性较好。

在国家级气象站1961—2013年的人工观测年雷暴日数序列中，仍可能存在因观测方式改变、台站迁移等原因造成的非均一点，后期可继续采用PMFT检验法对人工观测序列进行均一性检验，找出序列非均一的年份，分析其原因，并完成序列的均一性订正，形成均一的长序列年雷暴日数数据集，为雷暴气候变化分析、防雷工程设计、雷击风险评估和雷电科学研究提供科学依据。