引言

数值天气预报采用数学物理方法求解数值模型，基于初值和边界场条件，能够有效地预测大气未来的演变过程。资料同化技术是一种能够有效提高模式初始场准确度的方法(矫梅燕等，2006)。随着观测技术的不断提高，观测资料的种类和时空分辨率不断增加，近年来如何充分利用多源观测资料改善初始场准确度，已成为数值天气预报研究领域的热点(崔春光等，2021；师春香等，2019；龚建东，2013)。

强对流天气往往伴随着闪电活动，闪电活动是由于过程中存在强烈的上升气流，且强对流云环境条件有利于冰相水凝物粒子的生成。这些粒子在上升气流的作用下发生碰撞和分离，并最终产生闪电活动(Liu et al，2012)。观测研究表明，闪电活动在强对流天气过程中具有重要的指示作用(Makowski et al，2013；Rudlosky and Fuelberg, 2013)。总闪电率(地闪和云闪)与云中霰粒子的含量以及上升气流的强度密切相关(Ribaud et al，2016；Guo et al，2016)。近年来闪电探测技术得到了快速发展，为了实时探测闪电活动，已建成了许多地基和星载的闪电探测系统(王宇等，2015；Yang et al，2017)。在对流允许尺度上通过各种方法进行闪电资料同化研究上取得了显著进展。其中，基于牛顿松弛逼近(Nudging)技术的闪电资料同化方法通过调整水汽含量、潜热以及水凝物粒子质量混合比等，从而在有闪电观测但模式中无对流出现的位置促进对流的产生(Wang et al，2022；Sun et al，2022)。另外，基于变分同化技术和集合卡尔曼滤波(EnKF)技术的闪电资料同化方法，通过建立闪电观测数据与模式变量间的观测算子，对模式背景场进行调整(Xiao et al，2021；Kong et al，2020；徐国强等，2020；何邓新等，2021)。

地面气象观测站通过对气象要素的直接测量获得常规大气背景场信息，包括温度、压力、湿度和风速等，能够反映大、中尺度天气特征，对于数值天气预报具有重要意义。与其他非常规观测资料相比，地面观测资料具有多项优势，例如，地面观测要素与模式变量一致，质量稳定且在一定区域内具有较高的精度(何光碧等，2013)。近年来，研究人员利用各种方法对常规观测资料进行同化应用，取得了显著的成果。其中，三维变分方法被广泛应用于常规观测资料的同化，通过该方法可以在一定程度上改善地面要素以及降水的模拟效果(王平等，2017；盛春岩等，2018；张鑫宇等，2021)。此外，Nudging同化方法也被证实可以提高对雨带位置、强度的模拟准确度，并改善地面温度和风场的模拟效果(Bullock et al，2014；周易等，2016)。另外，具有流依赖特征的EnKF方法也被广泛应用于常规观测资料的同化研究中，改善了模式的模拟效果(邵长亮和闵锦忠，2019；孟智勇等，2019)。

地面常规观测在直接观测气象要素方面通常具有较高的精度。然而由于观测空间不连续性、观测高度受限、无法反映云的信息等缺陷，对其的同化应用存在诸多问题(Lin and Pu, 2020)。而闪电观测则范围广泛，受地表地形影响较小，且在空间和时间上具有较好的连续性，能够反映强对流系统对流区域的大气背景场信息。然而，闪电资料仅能反映强对流云中对流区域的信息，无法有效反映无云区域和层状云区域的大气状态，因此在某种程度上限制了闪电资料同化对强对流天气预报的作用(Wang et al，2018)。由于闪电资料和地面常规观测资料提供了互补信息，联合同化这两种资料相比于单独同化其中一种资料有进一步提高对强对流天气预报准确度的潜力。基于以上背景，本研究开展了对上述两种资料联合同化的研究，分析联合同化两种资料相比于只同化其中一种资料对中尺度对流系统模拟的改进作用。

1 资料、方法及模拟设置 1.1 资料

本文使用的闪电资料来源于广东省气象局和Earth Network公司建立的闪电探测系统，由16个地面传感器组成(图 1a)，其定位精度约为600 m，总闪电探测效率为77%(Zhang et al，2016)。为适应WRF模拟域网格(3 km×3 km)，将总闪电资料进行了插值处理，从原始经纬度坐标转换到相应的模拟域网格上。

地面自动气象站分布密集，观测频次高，其空间分辨率可达到3 km，利用其观测资料，包括气温、湿度、气压、风速、风向等气象要素进行同化。同化区域内共有约5200个观测站点(图 1a)，约有4000个站点的温度、3500个站点的风速、900个站点的水汽、1000个站点的气压参与同化过程。

为了对模拟结果进行评估，使用了覆盖广东省范围的雷达组合反射率因子产品和地面常规观测资料，对雷达反射率因子资料进行了插值处理，使其与数值试验的模拟域网格保持一致。

1.2 同化方案 1.2.1 闪电资料

本研究使用了Wang et al(2017)在WRF-FDDA(Four Dimensional Data Assimilation)系统中开发的闪电资料同化方案。该方案利用观测到的闪电率反演对流云内霰粒子的质量混合比，并充分考虑了水凝物粒子形成过程中潜热释放过程的影响，将观测到的闪电率与霰粒子质量混合比建立联系。在霰粒子质量混合比反演方案中，利用闪电率和霰粒子垂直积分质量直接的观测算子，建立了闪电率分档区间与霰粒子垂直积分质量分档区间之间的对应关系(Carey et al，2014)：

$ F=2.43 \times 10^{-8} M_{\mathrm{g}} $

(1)

式中：F为闪电率(单位：min^-1)，M_g为-40~-10℃的霰粒子垂直积分质量(单位：kg)。

模式气柱中不同霰粒子垂直积分质量分档区间的霰粒子质量混合比特征廓线可以用来反演霰粒子质量混合比的三维结构。通过对强对流天气个例的回测模拟，建立了一个霰粒子质量混合比廓线的数据库。当模式网格中某一像素的闪电率观测值确定后，可通过该闪电率对应的霰粒子垂直积分质量分档区间确定该像素的霰粒子质量混合比垂直廓线。

为了更准确地反映霰粒子分布情况，还采用了一种基于距离权重函数的方法来考虑闪电观测相邻区域中的霰粒子含量, 如式(2)所示：

$ q_{\mathrm{g}}(x)=\frac{\sum_{i=1}^N q_{\mathrm{g}}\left(x_i\right) \cdot W_i\left(r_i, R\right)}{N} $

(2)

式中：q_g(x)为观测的闪电密度为0次·km^-2的格点x上的霰粒子质量混合比(单位：g·kg^-1)，q_g(x_i)和N分别为格点x影响半径R范围内非零闪电观测值的格点x_i上的霰粒子质量混合比和格点数量。其中W_i为Cressman(1959)提出的空间距离权重函数，由式(3)计算所得：

$ W_i\left(r_i, R\right)= \begin{cases}\frac{R^2-r_i^2}{R^2+r_i^2} & 0<r_i \leqslant R \\ 0 & r_i>R\end{cases} $

(3)

式中r_i是两个网格点x和x_i之间的水平距离。经过测试，当R=6 km时产生的结果更加可靠，因此R值设为6 km。

在霰粒子质量混合比的同化中，使用FDDA系统的“Grid Nudging”模块将反演的三维霰粒子质量混合比以四维时间连续的方式同化到WRF模式中。“Grid Nudging”适用于高空间分辨率观测的松弛逼近方法，如式(4)所示：

$ \frac{\partial X}{\partial t}=P(X, t)+G \cdot W \cdot T(Y-X) $

(4)

式中：X为模式变量的初猜值，P(X, t)是WRF模式预报方程中的原始动力、微物理项，G是防止模式变量瞬间变化引起虚假噪音波的松弛系数，W为空间权重，T为观测时间和模式时间迟滞性的实际权重，Y是模式格点上对应的观测值。

潜热释放是指在霰粒子形成过程中释放的热量，其温度增量与霰粒子质量混合比增量成正比关系，如式(5)所示：

$ \Delta T=\frac{\Delta q_{\mathrm{g}}\left[L_{\mathrm{c}}\left(T_0\right)+L_{\mathrm{f}}\right]}{c} $

(5)

式中：ΔT和Δq_g分别为温度增量(单位：℃)和霰粒子质量混合比增量(单位：g·kg^-1)，L_c(T₀)为温度为T₀时的凝结潜热(单位：J·kg^-1)，L_f为0℃时的冻结潜热(单位：J·kg^-1)，c为水的绝对热容(单位：J·kg^-1·℃^-1)。

潜热释放带来的温度增量依然使用FDDA系统的“Grid Nudging”模块以四维时间连续的方式同化到WRF模式中，其公式如下：

$ \frac{\partial X}{\partial t}=P(X, t)+G \cdot W \cdot T \cdot \Delta T $

(6)

采用的潜热调整过滤方法为：为了防止在对流区域的边缘处激发过强的上升气流，潜热对应的温度增量只在闪电反演的霰粒子质量混合比大于1 g·kg^-1时才添加。

1.2.2 地面常规观测资料

采用WRFDA-3DVAR同化系统对地面常规观测资料进行同化。在该同化系统中，通过求解代价函数得到最优解，来表示分析时刻下的大气真实状态，如式(7)所示：

$ \begin{gathered} \boldsymbol{J}(\boldsymbol{x})=\frac{1}{2}\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{b}}\right)^{\mathrm{T}} \boldsymbol{B}^{-1}\left(\boldsymbol{x}-\boldsymbol{x}_{\mathrm{b}}\right)+ \\ \frac{1}{2}(H \boldsymbol{x}-\boldsymbol{y})^{\mathrm{T}} \boldsymbol{R}^{-1}(H \boldsymbol{x}-\boldsymbol{y}) \end{gathered} $

(7)

式中：J(x)为代价函数，x表示大气状态变量，x_b表示大气状态背景场，B为背景场x_b的误差协方差矩阵，H为观测算子，y为已知的观测值，R为观测值y的误差协方差矩阵。

采用的地面资料同化方案为Ruggiero方案(Ruggiero et al，1996)，该方案考虑了模式与实际观测站点地形之间的差异：当模式最低层高度低于地面观测站的高度时，则将地面观测资料作为高空资料同化进入模式；当模式最低层高度比观测站高度高，且高度差小于100 m时，则利用背景场信息将观测资料反演至模式最低层；当模式最低层高度超出地面观测站高度100 m时，则剔除该站点资料。风的同化基于莫宁-奥布霍夫相似理论(张鑫宇等，2021)进行订正；温度的同化则基于位温垂直变化率做订正；当相对湿度>100%或＜10%时，对其进行局部调整。

本研究中使用的是通过2017年6月研究区域一个月的模拟数据得出的具有区域性特点的背景误差协方差，而非全球范围通用，能更好地反映不同地区的大气变化特征和误差分布情况。背景误差协方差模拟中的模式设置以及微物理参数化方案与本研究资料同化试验中所使用的相同。

1.3 个例选取和试验设计

选取2017年6月16日广东省一次强对流过程进行模拟对比。16日03时(世界时，下同)，广东省上空出现了多个离散的强雷暴单体。随着时间的推移，单体逐渐向东南方向移动并不断加强。10时左右，部分雷暴单体逐渐合并形成了一个中尺度对流系统，并持续向东南方向移动。17时前后，中尺度对流系统形成的线状对流带移动至海上，此次对流过程结束。

使用WRF4.1.1模式，采用双重嵌套的三重网格进行数值模拟(图 2)。网格水平分辨率分别为27、9、3 km。共设43个垂直层，模式顶部设置在50 hPa高度。模式的初始场和边界条件来源于美国国家大气环境中心6 h间隔的0.25°×0.25°全球预报系统(GFS)再分析资料。Thompson微物理参数化方案(Thompson et al，2008)、MYJ大气边界层方案(Janjic '，1994)、Noah陆面模式(Ek et al，2003)、GCMs快速辐射传输模型(Iacono et al，2008) 以及Grell-Freitas积云参数化方案(Grell and Freitas, 2014)被应用到本研究的模拟设置中。其中，Grell-Freitas积云参数化方案在模拟域D01和D02中被使用，但在云可辨尺度模拟域D03中未被使用。

本研究所采用的联合同化方案如下：在同化周期中，通过Wang et al(2017)开发的基于Nudging技术的闪电资料同化方法对闪电资料进行四维连续同化，以15 min作为一个闪电累积窗口(即以15 min累积量计算闪电密度，并将该数据在对应的15 min的模式积分过程中连续同化入模式背景场中)。此外，在同化周期中，通过WRFDA-3DVAR系统每间隔1 h同化一次地面常规观测资料，即采用逐1 h循环同化更新的同化方式。闪电资料和地面常规观测资料的同化在模拟域D03中开展。

为了评估联合同化闪电和地面常规观测资料相比于只同化其中一种资料对中尺度对流系统模拟的改进作用，开展了以下数值模拟试验：未同化任何资料的控制试验(CTRL), 地面常规观测资料同化试验(CDA), 闪电资料同化试验(LDA), 地面常规观测资料和闪电资料的联合同化试验(LCDA)。所有试验均采用冷启动的方式进行预报。表 1给出了各试验的名称缩写及简要描述，流程如图 3所示。

1.4 预报技巧评分

为了对各组试验的预报结果进行定量评估，采用基于领域的预报技巧评分FSS(fraction skill score)(Roberts and Lean, 2008)，计算了各组试验模拟的雷达组合反射率因子的FSS：

$ \mathrm{FSS}=1-\frac{\frac{1}{N} \sum_{i=1}^N\left(p_{\mathrm{f}}-p_{\mathrm{o}}\right)^2}{\frac{1}{N} \sum_{i=1}^N\left(p_{\mathrm{f}}\right)^2+\frac{1}{N} \sum_{i=1}^N\left(p_{\mathrm{o}}\right)^2} $

(8)

式中：N为验证区域中的网格点数，p_f和p_o分别为预报和观测中的雷达组合反射率因子超过给定阈值的分数覆盖率(单位：%)。FSS的值越接近于1，说明模拟结果与观测结果越接近。

2 结果与分析 2.1 模拟结果分析

通过主观比较和预报技巧评分的方式，对三组同化试验(CDA、LDA、LCDA)进行评估，以探究联合同化闪电和地面常规观测资料相比只同化其中一种观测资料的改进效果。

将各组试验模拟的雷达组合反射率因子与地基雷达观测结果进行对比(图 4)。在同化期中，相比于CTRL和CDA，LDA和LCDA更好地模拟出了对流中心的位置，并成功模拟出了23.3°N、113°E附近的强对流，但是该区域的强对流在CTRL和CDA中未能被成功模拟(图 4红色箭头)。在LCDA和LDA中，闪电资料的同化缩短了模式预热时间，对流系统在同化期的第一个小时内就已经产生，而同时期的CDA和CTRL中，对流还未生成(图略)。

对于广东省东北部的虚假对流(图 4绿色箭头)，CDA中虚假对流被抑制，而在LDA中，抑制效果并不明显。相比之下，LCDA在准确模拟出对流活动的同时，对上述虚假对流抑制更明显。同时，LCDA还成功抑制了LDA中广东省西部地区出现的虚假强对流(图 4黑色箭头)。地面常规观测资料在同化周期内的引入，减少了大尺度范围内地面温度、风速、水汽的误差，改善了要素场的分布情况，进而减小了模式对对流活动的模拟误差，抑制了部分虚假对流的产生，使得对流活动的空间分布准确度也得到了提高。

在预报期内，CTRL和CDA中模拟的对流中心位置明显偏离于雷达观测的对流中心。相比之下，LDA和LCDA更准确地预报出了对流中心的位置和空间范围。然而，LCDA和LDA对于广东省西部对流区域的模拟均不理想，这是由于在同化窗口中该区域观测到的闪电非常少，缺少对流云的信息。此外LCDA中的虚假对流也要明显弱于CTRL和LDA试验中对应区域的虚假对流。由于大多数闪电产生于对流系统的对流区域，不能很好地反映层状云区域的信息(Lund et al，2009)，导致LCDA对于层状云区域模拟范围与观测值的结果偏差较大，而层状云区的信息对对流活动的模拟结果有很大影响(Peters and Schumacher, 2014)。随着时间的推移，一些对流单体在预报期中并没有被维持，模拟结果与观测结果之间逐渐出现偏差。

为了研究联合资料同化对模式热动力场的影响，解释其对强对流天气预报准确度提升的原因，对模式模拟的2 m气温和10 m风场与观测进行对比。如图 5、图 6中方框所示，CTRL中强对流部分区域的模拟温度比观测结果高4~6 K，即其模拟的冷池范围小于实际观测。模式初始温度场的偏差和相应的冷锋误差是导致CTRL的对流过程模拟失败的主要原因。相比于CTRL，LDA中的冷池范围变大，但仅对对流发生区域的温度场进行了调整。而CDA对整个模拟区域的温度场进行了调整，降低了整个模拟区域2 m气温的误差，但冷池范围相比LDA更小。相比于其他三组试验，LCDA在对流区域2 m气温误差为1~2 K，冷池范围更大，更加接近观测值，且其阵风锋的位置更接近实况。LCDA通过同化大范围的地面观测资料调整了非对流区域的温度场，进一步提高了地面冷池的准确度。

图 7为各试验分析场中的10 m经向风风速和纬向风与地面观测结果的偏差。相较于CTRL，CDA降低了整个模拟区域10 m经向风和纬向风的误差，大部分地区减少了约2 m·s^-1。LDA对风场的改善主要体现在强对流区域，整体上对于水平风场的改善效果不及CDA。LCDA对10 m经向风场的误差减小更为显著。CTRL中强对流部分区域的10 m经向风风速的误差约4~5 m·s^-1，而这一误差在LCDA中被减小到了2~3 m·s^-1。因此，相比于LDA，LCDA进一步减小了水平风场的误差。

如图 8所示，在CTRL中，广东省西南部地区的大气偏干，水汽含量相对较少，不利于强对流的产生，而LCDA试验同化周期中广东西南部的水汽更充沛，为中尺度对流系统的发展提供了水汽条件。LDA试验中的水汽质量混合比仅在闪电发生区域得到相应调整，在没有闪电发生的区域对水汽误差改善并不明显。CDA试验通过同化地面常规观测资料，地面水汽的整体分布得到调整。在LCDA试验中，强对流区域的地面水汽质量混合比误差得到进一步减小，约为0.001~0.002 kg·kg^-1，水汽的整体分布也更加接近观测值。LCDA通过补充大范围地面常规观测资料的水汽信息，较LDA进一步缩小了大范围水汽场的误差。然而，在后续预报期内，中尺度对流系统北部的水汽含量较少，导致向东南方向推进的阵风锋未能持续激发出新的对流(图略)。

图 9为各组循环同化试验中2 m气温和2 m水汽质量混合比的均方根误差(RMSE)。在同化期中，对于气温和水汽的模拟，各组同化试验相较于CTRL都有一定程度的改善。在大部分时刻，联合同化试验LCDA中温度、水汽含量的均方根误差更小，例如，10时CTRL中的2 m气温RMSE为2.48 K，LDA、CDA和LCDA的RMSE分别为1.92、1.59和1.46 K。此外，LCDA和CDA中水汽质量混合比的RMSE较其他试验也更小。

2.2 预报技巧评分

图 10为以24 km容错半径(即8个模式网格点)和40 dBz阈值为基准的FSS技巧评分结果。40 dBz的雷达组合反射率因子被选择为表征对流区域的反射率因子阈值。同化期内联合同化试验LCDA的FSS评分较其他试验更高。CDA由于缺少闪电资料提供的对流区域内的信息，FSS评分低于包含闪电资料同化的试验(LDA和LCDA)。这是由于闪电资料同化在观测到闪电的区域引入了潜热增量，促进了上升运动，此外同化过程中引入的霰粒子质量混合比增量通过冷云降水过程更快地在背景场中形成冷池和阵风锋，促进对流系统的进一步发展。在预报期的前4 h内，CDA的FSS评分高于LDA，LCDA的FSS评分仍然略高于LDA和CDA(图 10)。这是由于联合资料同化窗口中模拟出的大气温、压、湿、风以及微物理信息，以及更准确的对流云、阵风锋的位置使得LCDA在预报周期中的模拟也更为准确。综上所述，相比于只同化一种观测资料，闪电资料和地面常规观测资料联合同化可以进一步提高模式对对流区域的预报准确度。

3 结论与讨论

本研究以2017年6月16日发生在广东省的一次中尺度对流过程为例，分析了联合同化上述两种资料相比于只同化其中一种资料对中尺度对流系统模拟和预报的改进作用，结果显示：

(1) 相比于无资料同化的控制试验，单独同化地面常规观测资料、闪电观测资料以及联合同化地面常规观测资料和闪电资料均提高了对该中尺度对流系统的模拟和预报准确度。

(2) 相比于只同化地面常规观测资料的试验，联合同化常规和闪电观测资料试验更好地模拟出了对流中心的位置和整体形态，并成功模拟出了漏报的对流。在同化窗口中，联合同化试验的FSS预报技巧评分显著高于地面常规观测资料同化试验，在之后的预报周期中其FSS预报技巧评分整体略高于地面常规观测资料同化试验。这是由于闪电资料的引入通过潜热、霰粒子质量混合比增量，提高了背景场中上升气流、冷池和阵风锋的准确度，进而改善了对该对流系统的模拟准确度。

(3) 相比于只同化闪电观测资料，联合同化试验中地面常规观测资料的引入减小了更大范围温度、水汽、风场的背景场误差，抑制了部分地区虚假对流的发展，整体上提高了对该对流系统的模拟准确度。

需要指出的是，本研究为个例研究，后续的研究中将对连续的、大量的、不同类型的强对流天气个例进行试验，对闪电和地面常规观测资料联合同化的效果进行更为系统性的评估。此外，使用具有流依赖背景场误差协方差特征的EnKF同化方法对闪电资料和地面常规观测资料进行联合同化也是下一步研究的重点内容。