一种通过同化闪电资料改进强对流预报的方法

(19)中华人民共和国国家知识产权局(12)发明专利申请(10)申请公布号CN112558188A(43)申请公布日2021.03.26(21)申请号202110089247.5(22)申请日2021.01.22(71)申请人兰州大学地址730000甘肃省兰州市天水南路222号(72)发明人甘茹蕙　杨毅　郭树昌　(74)专利代理机构兰州中科华西专利代理有限公司62002代理人曹向东(51)Int.Cl.G01W1/10(2006.01)G06F30/20(2020.01)权利要求书2页说明书8页附图5页(54)发明名称一种通过同化闪电资料改进强对流预报的方法(57)摘要本发明涉及一种通过同化闪电资料改进强对流预报的方法，该方法包括以下步骤：⑴建立闪电频数和最大垂直速度之间的观测算子；⑵提取同化时刻前20分钟的闪电资料，将闪电资料转化为1 km的网格上的最大垂直速度，然后将该1 km的最大垂直速度分配到数值模式网格上；⑶根据GGF方法，利用背景场计算得到各个变量的垂直局地化函数；⑷利用集合均方根滤波方法，结合垂直局地化函数同化代理最大垂直速度，循环同化四次后得到分析场；⑸利用分析场继续预报得到预报场；⑹进行预报检验。本发明通过同化总闪资料有效地改进了背景场中的动力场，并且对湿度场、温度场也有一定的改进，方法合理有效，且显著改进了数值模式的降水及回波的预报效果。CN112558188ACN112558188A权　利　要　求　书1/2页1.一种通过同化闪电资料改进强对流预报的方法，包括以下步骤：⑴根据对流活动中闪电频数和云顶高度之间的关系为：，而最大垂直速度与云顶高度之间的关系为，建立闪电频数和最大垂直速度之间的观测算子：；式中：wmax 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示最大垂直速度，单位：m/s；F为每公里每分钟的闪电频数； Z为云顶高度，单位：km；⑵提取同化时刻前20分钟的闪电资料，根据所述闪电频数和最大垂直速度之间的观测算子将所述闪电资料转化为1 km的网格上的最大垂直速度，然后将该1 km的最大垂直速度分配到数值模式网格上；⑶根据GGF方法，利用背景场计算得到各个变量的垂直局地化函数；⑷利用集合均方根滤波方法，结合所述垂直局地化函数同化代理最大垂直速度，循环同化四次后得到分析场；⑸利用所述分析场继续预报得到预报场；⑹进行预报检验。2.如权利要求1所述的一种通过同化闪电资料改进强对流预报的方法，其特征在于：所述步骤⑶中垂直局地化函数按下述方法计算：①令N表示样本大小，L表示观测数，M表示模式的垂直层，表示来自第n个样本的第m层的第l个观测处的模式状态变量；其中，，；是的集合平均；②按下述方法根据背景场构造一组观测扰动： ，；式中：代表来自第n个样本的第l个观测；是由闪电转换的第l个最大垂直速度；是计算得到的第n个样本的观测扰动；是的平均，其中是的平均，，是观测算子；③观测和模式变量之间的相关系数可以通过以下公式计算：；式中：为相关系数；代表的集合平均；④相关系数随机地分为G个组，每个组有Q个样本，则相关系数表示为，以及，所有的G个组中Q个样本的估计相关的均方根差表示为：2CN112558188A权　利　要　求　书2/2页；均方根最小化得到：，是计算得到的自适应的垂直局地化函数，即GGF函数； G取16。3.如权利要求1所述的一种通过同化闪电资料改进强对流预报的方法，其特征在于：所述步骤⑷集合均方根滤波方法中的观测变量是最大垂直速度wmax；分析变量包括风场、扰动位温、扰动位势、水汽混合比、雨水混合比、冰混合比、霰混合比、云混合比、雪混合比。3CN112558188A说　明　书1/8页一种通过同化闪电资料改进强对流预报的方法技术领域[0001]本发明涉及气象数据处理方法，尤其涉及一种通过同化闪电资料改进强对流预报的方法。背景技术[0002]尽管天气预报不断改进，但由于这些事件涉及的尺度和物理过程的多样性，实现对流活动预警的准确性仍然是一个很大的挑战。强对流天气通常伴有频繁的闪电活动，闪电可用于监测中尺度对流活动的发生和发展，闪电资料同化也是一种非常有用的方法，可以提高初始条件的准确性，从而改善天气预报。[0003]Browning（1989）提出调整中尺度数值天气预报模式初始场的一种可能方法是改进水汽分析，目前许多闪电资料同化研究都是通过改善初始条件下的水汽场进而改进预报。Alexander等（1999）基于降雨率和闪电之间的关系实现了闪电资料的同化。Fierro等（2012）提出可以通过改变混合区（0°C和‑20°C）的水汽混合比例来同化对流分辨尺度下的闪电资料。自此之后，基于Fierro的方法开展了很多闪电资料同化的工作，结果表明：无论是采用不同的同化方法同化代理水汽混合比，或者是基于不同的模式同化代理水汽混合比都可以改善预报。[0004]此外，Qie等（2014）建立了总闪电率与冰相粒子（霰、冰、雪）混合比之间的经验公式，用以调整混合区冰相粒子的混合比例。Chen 等（2018）利用Fierro和Qie的方法，利用闪电资料同时调整了混合区低层水汽量和霰质量，并且考虑了模型网格内的热力学和动力学条件，维持模型平衡。此外，Allen等（2016）和Kong等（2020）也利用闪电和霰质量、霰体积之间的观测算子同化了闪电资料。[0005]调整数值模式初始场中的垂直速度分布是改善强对流预报的另一种可行的方法。目前的研究多集中在水汽和水物质的同化上，考虑同化风场特别是垂直速度的研究较少。已有研究表明，上升气流强度与观测到的闪电频数有定性的相关性(Deierling和Petersen, 2008)。根据闪电观测结果，Williams和Earl（1985）表明，陆地对流过程中的总闪频率与对流云高度有大约五次方关系。Price和Rind（1992）指出，根据许多对流事件，对流云中的最大上升速度与云高相关，他们分别建立了大陆云和海洋云中最大垂直速度和总闪频率的经验公式。Wang等（2020b）根据Price和Rind的公式，采用Nudging方法同化闪电数据，调整动态场。此外，很少有研究旨在基于闪电和上升气流的关系来优化动态场。[0006]集合均方根滤波（EnSRF）使用预测的集合来估计流依赖的背景误差协方差，它已被用于同化常规观测、雷达、卫星和闪电数据（Gao和Min，2018；Wang等，2015）。局地化（Anderson，2007b；Hamill等人，2001）和协方差膨胀（Anderson，2007a，2008）对于EnSRF来说是必不可少的，因为它们会受到小集合数导致的采样误差的影响（Lei和Anderson，2014）。常用的局地化函数是一个近似高斯的五阶多项式函数（Gaspari和Cohn，1999，简称GC函数），GC的宽度参数需要针对给定的情况进行调整，这对于各种观测和状态变量来说是不实用的。此外，还有一些关于自适应协方差的理论研究（Anderson和Lei，2013；Lei等，4CN112558188A说　明　书2/8页2014；Lei和Anderson，2014），但关于二维变量最大垂直速度、垂直局地化的概念并没有很好的定义。Lei等（2016）提出了一种GGF（global group filter）方法，可以提供垂直局地化函数的自适应估计，随后，GGF定位函数被应用于ANSU‑A辐射数据同化中，获得了较好的预报效果（Lei等，2016；Wang等，2020a）。[0007]本发明基于2019年7月6日发生的强对流事件，根据总闪和最大垂直速度之间的“瞬时”关系，使用EnSRF方法以及GGF局地化方法同化由总闪数据转换的最大垂直速度（EnSRF_Wmax），以期改善强对流过程的预报。发明内容[0008]本发明所要解决的技术问题是提供一种合理有效的通过同化闪电资料改进强对流预报的方法。[0009]为解决上述问题，本发明所述的一种通过同化闪电资料改进强对流预报的方法，包括以下步骤：⑴根据对流活动中闪电频数和云顶高度之间的关系为：，而最大垂直速度与云顶高度之间的关系为，建立闪电频数和最大垂直速度之间的观测算子：；式中：wmax表示最大垂直速度，单位：m/s；F为每公里每分钟的闪电频数； Z为云顶高度，单位：km；⑵提取同化时刻前20分钟的闪电资料，根据所述闪电频数和最大垂直速度之间的观测算子将所述闪电资料转化为1 km的网格上的最大垂直速度，然后将该1 km的最大垂直速度分配到数值模式网格上；⑶根据GGF方法，利用背景场计算得到各个变量的垂直局地化函数；⑷利用集合均方根滤波方法，结合所述垂直局地化函数同化代理最大垂直速度，循环同化四次后得到分析场；⑸利用所述分析场继续预报得到预报场；⑹进行预报检验。[0010]所述步骤⑶中垂直局地化函数按下述方法计算：①令N表示样本大小，L表示观测数，M表示模式的垂直层，表示来自第n个样本的第m层的第l个观测处的模式状态变量；其中，，；是的集合平均；②按下述方法根据背景场构造一组观测扰动： ，；式中：代表来自第n个样本的第l个观测；是由闪电转换的第l个最大垂直速度；是计算得到的第n个样本的观测扰动；是的平均，其中是的平均，，是观测算子；5CN112558188A说　明　书3/8页③观测和模式变量之间的相关系数可以通过以下公式计算：；式中：为相关系数；代表的集合平均；④相关系数随机地分为G个组，每个组有Q个样本，则相关系数表示为，以及，所有的G个组中Q个样本的估计相关的均方根（RMS）差表示为：；均方根最小化得到：，是计算得到的自适应的垂直局地化函数，即GGF函数； G取16。[0011]所述步骤⑷集合均方根滤波（EnSRF）方法中的观测变量是最大垂直速度wmax；分析变量包括风场（u,v,w）、扰动位温（prt）、扰动位势（ph）、水汽混合比（qv）、雨水混合比（qr）、冰混合比（qi）、霰混合比（qg）、云混合比（qc）、雪混合比（qs）。[0012]本发明与现有技术相比具有以下优点：1、本发明基于中尺度数值模式WRF，通过同化闪电资料来改进背景场动力信息的方法（简称EnSRF_Wmax），建立了总闪和最大垂直速度之间的观测算子，利用EnSRF方法同化由闪电资料转换的代理最大垂直速度，并采用GGF垂直局地化方案，有效改进了背景场中的动力信息，不但方法合理有效，而且同时对于背景场中的水汽场、温度场也有一定的调整，进而有效地改进了数值模式对于对流过程中的降水和回波的预报。[0013]2、相比于同化三维资料，本发明同化二维的最大垂直速度减小了计算量。附图说明[0014]下面结合附图对本发明的具体实施方式作进一步详细的说明。[0015]图1为本发明中依据GGF垂直局地化方案计算的模式基本状态量的垂直局地化函数图。其中：（u）、（v）、（w）、（prt）、（ph）、（qv）、（qr）、（qc）、（qg）、（qi）、（qs）分别表示u分量风、v分量风、垂直速度、扰动位温、扰动位势、水汽混合比、雨水混合比、云水混合比、霰混合比、冰混合比、雪混合比。[0016]图2为本发明中闪电转换为模式网格上的最大垂直速度的示意图。其中：（a）是观6CN112558188A说　明　书4/8页测到的20 min内地总闪；（b）是将闪电分配到1 km的网格上的闪电频数；（c）是在1 km网格上将闪电频数转换为最大垂直速度；（d）是将1 km网格上的最大垂直速度分配到模式模拟网格（3 km）。[0017]图3为本发明中同化时刻前20 min的闪电分布，以及转换得到的最大垂直速度图。其中：（a1）是2019年7月6日03 UTC前20 min的闪电分布；（b1）是2019年7月6日04 UTC前20 min的闪电分布；（c1）是2019年7月6日05 UTC前20 min的闪电分布；（d1）是2019年7月6日06 UTC前20 min的闪电分布；（a2）是根据闪电频数转换的2019年7月6日03 UTC的最大垂直速度；（b2）是根据闪电频数转换的2019年7月6日04 UTC的最大垂直速度；（c2）是根据闪电频数转换的2019年7月6日05 UTC的最大垂直速度；（d2）是根据闪电频数转换的2019年7月6日06 UTC的最大垂直速度。[0018]图4为本发明中观测模拟的累积降水分布图。其中：（a1）是观测的2019年7月6日06 UTC至09 UTC的累积降水；（a2）是控制实验（CTL）模拟的2019年7月6日06 UTC至09 UTC的累积降水；（a3）是同化试验（LDA_GGF）模拟的2019年7月6日06 UTC至09 UTC的累积降水；（b1）是观测的2019年7月6日09 UTC至12 UTC的累积降水；（b2）是控制实验（CTL）模拟的2019年7月6日09 UTC至12 UTC的累积降水；（b3）是同化试验（LDA_GGF）模拟的2019年7月6日09 UTC至12 UTC的累积降水。[0019]图5为本发明中降水预报性能图。其中：（a）是阈值取1 mm时降水预报的性能图；（b）是阈值取5 mm时降水预报的性能图。具体实施方式[0020]一种通过同化闪电资料改进强对流预报的方法，包括以下步骤：⑴由于闪电的观测资料（频数）不是数值模式的模式变量，对其进行同化时首先要建立闪电和其他模式变量的观测算子。[0021]已有的研究表明对流活动中闪电频数和云顶高度之间存在五次幂的关系，Ushio等(2001)利用TRMM卫星上降水雷达和闪电成像仪的数据，获得了“瞬时”闪电频数和对流云云顶高度之间的关系为，而最大垂直速度与云顶高度之间的关系为，因此，建立闪电频数和最大垂直速度之间的观测算子：；式中：wmax表示最大垂直速度，单位：m/s；F为每公里每分钟的闪电频数；Z为云顶高度，单位：km。[0022]此观测算子仅适用于1 km的网格。[0023]⑵提取同化时刻前20分钟的闪电资料，根据闪电频数和最大垂直速度之间的观测算子将闪电资料转化为1 km的网格上的最大垂直速度，然后将该1 km的最大垂直速度分配到数值模式网格上。[0024]具体过程：第一步将闪电资料分配到1 km的网格，第二步根据得到的观测算子将1 km的闪电频数转换为最大垂直速度，最后再将最大垂直速度转换为模式网格（3 km），具体过程如图2所示。[0025]⑶根据GGF（global group filter）方法，利用背景场计算得到各个变量的垂直局7CN112558188A说　明　书5/8页地化函数。利用EnSRF方法同化最大垂直速度难点在于垂直局地化方案的设置，本发明中垂直局地化方案采用GGF方案，GGF局地化方案观测和状态变量的扰动集合，利用样本相关计算得到。[0026]垂直局地化函数按下述方法计算：①令N表示样本大小，L表示观测数，M表示模式的垂直层，表示来自第n个样本的第m层的第l个观测处的模式状态变量；其中，，；是的集合平均；②在EnSRF中没有观测扰动，为了计算GGF垂直局地化函数，需要按下述方法根据背景场构造一组观测扰动： ，；式中：代表来自第n个样本的第l个观测；是由闪电转换的第l个最大垂直速度；是计算得到的第n个样本的观测扰动；是的平均，其中是的平均，，H'是观测算子；在本研究中H'也就是利用背景场计算最大垂直速度的方法。[0027]③观测和模式变量之间的相关系数可以通过以下公式计算：；式中：为相关系数；代表的集合平均；④相关系数随机地分为G个组，每个组有Q个样本，则相关系数表示为，以及，所有的G个组中Q个样本的估计相关的均方根（RMS）差表示为：；这里的q,i,j没有实际的意义，就是数组算法的表示，求和算法，q从1取到Q，i从1取到G，j从1取到G，其中j不等于i。[0028]均方根最小化得到：，是计算得到的自适应的垂直局地化函数，即GGF函数； G取16。[0029]计算得到的模式基本量的局地化函数的归一结果如图1所示。[0030]⑷利用集合均方根滤波（EnSRF）方法，结合垂直局地化函数同化代理最大垂直速8CN112558188A说　明　书6/8页度，循环同化四次后得到分析场。[0031]由于最大垂直速度是一个特殊的二维的变量，集合均方根滤波方法（EnSRF）易于建立最大垂直速度和其他变量之间建立平衡约束，此外EnSRF方法避免了产生扰动观测时生成的采样误差，所以本发明中选用EnSRF同化闪电转换的最大垂直速度。[0032]在EnSRF中，上标a，b，o分别表示分析场、背景场以及观测，xb是模式预报的背景场，yo为一系列观测，分析场xa的最小方差估计根据卡尔曼滤波方程可以表示为：，，其中：H是观测算子H'的线性化形式，观测算子H'是状态量x到观测yo之间投影关系，K是卡尔曼增益矩阵，Pb是背景误差协方差矩阵，R是观测误差协方差矩阵，在EnSRF中没有观测扰动（），集合平均（）和第i个集合的偏差（）可以表示为：，；其中 称为“简化”卡尔曼增益矩阵，定义如下：。[0033]在本发明中，观测变量是最大垂直速度wmax；分析变量包括风场（u,v,w）、扰动位温（prt）、扰动位势（ph）、水汽混合比（qv）、雨水混合比（qr）、冰混合比（qi）、霰混合比（qg）、云混合比（qc）、雪混合比（qs）。[0034]⑸利用分析场继续预报得到预报场。[0035]⑹进行预报检验。[0036]为了验证本发明对强对流预报的改善效果，选取了2019年7月6日发生在我国山东和江苏地区的一次强对流过程进行个例分析，本次试验总共设计一组控制实验以及一组循环同化试验，详细的试验设计如表1所示：表1： 试验设计表9CN112558188A说　明　书7/8页表1中spin‑up的理解为：因为中小尺度模式的都是由无云初始场开始的，所以开始需要一定的时间产生相应的云水信息，冷启动之后会在一段时间内产生实际的云水信息，但是这几个小时内的预报是不准确的，将这几个小时称为“spin‑up”时间。也就是说这段时间的模拟结果是不可靠的，所以分析的时候不对spin‑up这期间的结果进行分析。[0037]在上述试验中，试验Ⅰ（CTL）是控制试验，试验Ⅱ（LDA_GGF）是利用EnSRF方法同化闪电的试验，同化时间是2019年7月6日03、04、05、06 UTC。本研究选用WRF 3.6.1版本，21个全国集合成员预报系统的集合预报被用为WRF模式提供初始场和边界场，为了使样本离散度增强，选用RRTMG和Goddard辐射方案使21个样本变成42个样本，除此之外，选用的参数化方案有Thompson微物理方案、Monin‑Obukhov近地层方案、Mellor‑Yamada‑Janjic TKE行星边界层方案、RUC陆面方案，Tiedtke积云参数化方案（只在D01中使用）。本研究采用两重嵌套网格，D01的水平分辨率为9 km，模式网格为171×151，D02水平分辨率为3 km，模式网格为391×271，垂直方向有50层。[0038]本研究选用EnSRF方案同化最大垂直速度，垂直局地化方案选用的是GGF局地化方案，图1所示是根据GGF局地化方案计算的模式状态量的局地化函数。本研究还建立了最大垂直速度和闪电频数之间的观测算子，图2 是将闪电观测转换为模式网格上最大垂直速度的过程，图3所示是本次试验中用到的闪电资料的分布，图3（a1‑d1）是同化时刻整点前20 min的累积闪电，图3（a2‑d2）是依据建立的观测算子转换得到的最大垂直速度，可以看出这次对流过程发生了很多次闪电过程，03 UTC闪电多集中在山东地区，之后对流过程以较快的速度向南移动，至06 UTC闪电主要分布在江苏的北部。[0039]图4是观测和模拟的累积降水分布，第一行是2019年7月6如06‑09 UTC的累积降水，第二行是09‑12 UTC的累积降水，第一列是观测降水，第二列是控制实验模拟的降水，第三列是同化试验模拟的降水，06‑09 UTC观测降水主要集中在江苏北部，控制实验没有模拟出江苏地区主要的雨带，而同化试验较好地改进了江苏地区的降水，09‑12 UTC观测的降水主要集中在江苏的南部，而模拟的降水位置偏北，同化试验模拟的雨带和观测的雨带位置比较接近。[0040]图5是降水预报性能图，图（a）和图（b）分别表示选取1 mm、5 mm阈值的结果，预报性能图可以反映很多信息，包括预报偏差（FR，图中虚线所示，越接近1表明预报偏差越小，大于1表示高估，小于1表示低估）、命中率（POD，越高表明预报越好）、临界成功指数（CSI，图中曲线所示，范围0‑1，约接近于1越好）、成功率（SR，1减去空报率，越接近于1越好），综合来说在预报性能图中约接近于右上角，表明预报效果越准确。图中浅色圆圈表示控制试验结果，深色表示同化试验结果，可以看出对于较小的阈值，同化试验在6小时预报都显著优于10CN112558188A说　明　书8/8页控制实验，当阈值取5 mm时，同化试验在前4个小时的预报效果较好。总体来看同化试验预报效果显著优于控制试验。11CN112558188A说　明　书　附　图1/5页图112CN112558188A说　明　书　附　图2/5页图213CN112558188A说　明　书　附　图3/5页图314CN112558188A说　明　书　附　图4/5页图415CN112558188A说　明　书　附　图5/5页图516