气象数据高精度融合技术研究

Research on High-Precision Fusion Technology for Precipitation Observation

摘要: 降水观测在气象领域至关重要，本文利用变分同化方法结合雷达与雨量计数据提升降水场估测精度，通过构建变分同化模型，并优化数据融合过程，实现了高精度降水估计。研究成果为气象预报、水文学及洪水预警等领域提供了可靠数据支持，并为未来研究提供了参考价值。

Abstract: Precipitation observation is crucial in meteorology. This study explores the strategy of enhancing precipitation field estimation accuracy through variational data assimilation (VDA) methods, combining radar and rain gauge data. By constructing a variational assimilation model and optimizing the data fusion process, high-precision precipitation estimation is achieved. The research findings provide reliable data support for meteorological forecasting, hydrology, and flood warning, and point to future research directions.

1. 引言

在气象领域，降水观测的精确性对于水文学、洪水预报等应用至关重要[1] [2]。雨量计，作为点观测的典范，以其高精度著称，能够精确捕捉局部地区的降水情况[3]。张欣[4]等使用对比实验验证了雨量计监测的合理性，雷梦佳[5]等制作了智能防堵磁致式雨量计原型样机，能达到与人工完全一致的雨量监测结果。但雨量计覆盖范围始终有限且数据不连续，难以满足大尺度降水估计的需求。相比之下，雷达观测具有大范围、连续性和快速响应的特点[6]。王文丽[7]将经多点校准法校准的雷达估测降雨数据、卡尔曼滤波校准后的雷达估测降雨数据输入到模型中进行径流模拟，评价其适用性。但雷达的估计精度存在明显不足，难以满足准确预测降雨量的需求。

为了克服单一观测手段的不足，学者们致力于发展数据同化技术，黄静[8]采用间接同化方式来同化雷达反射率因子，使得结果更接近真实值。张新[9]采用卡尔曼滤波算法对测雨雷达数据进行同化分析，降水误差得到有效降低，相比较而言，李梦迪[10]等验证了雷达融合雨量计降水估测方法，其准确性比单一手段同化技术高，且对极端降水有较好的捕捉和监测能力。

综上所述，本文综合分析雨量计及雷达的观测数据，结合L-BFGS-B迭代算法，为气象预报、水资源管理及防灾减灾提供更加准确、可靠的降水信息支持。

2. 模型基础

2.1. 模型背景与假设

背景1：现有两部具有不同探测精度雷达(文中表示为左雷达和右雷达)对同一区域的观测数据。左边雷达的降水测量误差为50%，结构估计误差为15%；右边雷达的降水测量误差为20%，结构估计误差为10%。

背景2：在背景1的基础上增加N点雨量计观测数据，建立变分同化模型。

设立背景的主要目的是应对一种特定情境，当已有观测数据受到限制，数据稀缺或质量不高的实际情况下，能够借助文中提出的方法来有效解决问题。这样的简化不仅有助于提升方法的普适性和易用性，还能促进其在不同领域和场景下的广泛应用。

假设：降水场结构仅考虑降水场的梯度，两部雷达资料区域均为正方形，来自同一种降水观测手段的权重是相等的。

2.2. 符号说明

本文所需的变量及其含义如表1所示。

Table 1. Variables and their meanings

表1. 变量及其含义

符号	含义说明
$x$	需要优化的降水场
$x_b$	雷达观测区域的背景场
$x_c$	雨量计的观测数据
$B$	背景误差协方差矩阵
$y_1$ , $y_2$	左雷达和右雷达的观测数据
$H_1$ , $H_2$	左雷达和右雷达的观测算子
$R_1$ , $R_2$	左雷达和右雷达的观测误差协方差矩阵
$\sigma$	雷达的相对误差的均方根
${\sigma }_{grad}$	雷达的相对梯度误差的均方根

3. 构建双雷达变分同化数学模型

3.1. 构建背景误差协方差矩阵B

为了精确估算背景误差协方差矩阵B，可以利用背景1中两个独立雷达观测数据间的差异作为依据，选取误差相对较小的雷达观测数据作为背景场，构建背景误差协方差矩阵，更精准地描述观测目标的背景状态不确定性。具体计算步骤如下公式所示：

误差：

$\text{error}={\text{left}}_{\text{rader}}-{\text{right}}_{\text{rader}}$ (1)

误差均值：

(2)

误差协方差：

$\text{ C(error)=}\left(\text{error}-{\text{error}}_{\text{mean}}\right){\left(\text{error}-{\text{error}}_{\text{mean}}\right)}^{\text{T}}$ (3)

添加正则化项：

$B=C\left(\text{error}\right)+\lambda I$ (4)

其中$I$ 表示单位矩阵，$\lambda$ 表示正则化因子，添加正则化项的目的是确保$B$ 是正定矩阵，避免在矩阵求逆时遇到问题[10]。

3.2. 目标函数的确定

通过将两部具有不同探测精度雷达对同一区域的观测数据分别作为观测数据和背景场数据，最小化一个目标函数，从而得到最优的降水场估计[6]，目标函数表达公式如下：

$\begin{array}{c}J\left(x\right)=\frac{1}{2}{\left(x-x_b\right)}^{\text{T}}{B}^{-1}\left(x-x_b\right)+\frac{1}{2c_1}{\left(y_1-H_{1}x\right)}^{\text{T}}\left(y_1-H_{1}x\right)+\frac{1}{2c_2}{\left(y_2-H_{2}x\right)}^{\text{T}}\left(y_2-H_{2}x\right)\\ +\frac{1}{2c_3}{\left(\nabla y_1-\nabla x\right)}^{\text{T}}\left(\nabla y_1-\nabla x\right)+\frac{1}{2c_4}{\left(\nabla y_2-\nabla x\right)}^{\text{T}}\left(\nabla y_2-\nabla x\right)\end{array}$ (5)

其中，$x$ 是需要优化的降水场，$x_b$ 是背景场，$B$ 是背景误差协方差矩阵，$y_1$ 和 $y_2$ 分别是观测场和背景场的观测数据，$H_1$ 和$H_2$ 是观测算子，假设为单位矩阵，$R_1$ 和 $R_2$ 是观测场和背景场的观测误差协方差矩阵，$c_1$ 、$c_2$ 、$c_3$ 、$c_4$ 分别是观测场的观测误差和梯度误差。

目标函数的梯度函数定义为：

$\begin{array}{l}\begin{array}{c}\nabla J\left(x\right)=B^{-1}\left(x-x_b\right)-\frac{1}{c_1}{H}_1^{\text{T}}\left(y_1-H_{1}x\right)-\frac{1}{c_2}{H}_2^{\text{T}}\left(y_2-H_{2}x\right)\end{array}\\ -\frac{1}{c_3}\left(\nabla x-\nabla y_1\right)-\frac{1}{c_4}\left(\nabla x-\nabla y_2\right)\end{array}$ (6)

3.3. 优化迭代算法L-BFGS-B (Limited-Memory Broyden-Fletcher-Goldfarb-Shanno with Box Constraints)

L-BFGS-B [11]用来最小化目标函数，从而找到最优的降水场估计。其基本思想是通过近似目标函数的二阶导数(即Hessian矩阵)来加速收敛过程。本文设置L-BFGS-B算法的最大迭代次数为100、梯度容忍度为1e-5。

3.4. 模型评价指标

1) 观测误差：

双雷达的相对误差的根方差$\sigma$ ：

$\sigma =\sqrt{\frac{1}{N}{\displaystyle \sum \frac{{\left(x_{\text{real}}-x_{\text{observe}}\right)}^2}{x_{\text{real}}^2}}}$ (7)

其中N为观测区域观测点个数，$x_{\text{real}}$ 是观测区域的实际降水量，$x_{\text{observe}}$ 是观测区域的观测降水量。由于观测场范围很大，降水量数据很多，所以观测数据准确时计算出的相对误差的根方差的应当贴近$\sigma$ ，因此可得：

$P_1=\frac{1}{2}{\left(\sqrt{\frac{1}{N}{\displaystyle \sum \frac{{\left(x_{\text{left}}-y_1\right)}^2}{x_{\text{left}}^2}}}-{\sigma }_{\text{left}}\right)}^2+\frac{1}{2}{\left(\sqrt{\frac{1}{N}{\displaystyle \sum \frac{{\left(x_{\text{right}}-y_2\right)}^2}{x_{\text{right}}^2}}}-{\sigma }_{\text{right}}\right)}^2$ (8)

其中$x_{\text{left}}$ ，$x_{\text{right}}$ 分别表示同化后降水场的左区域和右区域，${\sigma }_{\text{left}}$ ，${\sigma }_{\text{right}}$ 分别表示左雷达和右雷达的观测误差。

2) 梯度误差：

梯度误差反映的是降水场的结构特征，模型同化后的降水场的结构与观测场的结构的相对误差的根方差也应当尽可能地贴近雷达的梯度误差的根方差：

$P_2=\frac{1}{2}{\left(\sqrt{\frac{1}{N}{\displaystyle \sum \frac{{\left(\nabla x_{\text{left}}-\nabla y_1\right)}^2}{\nabla x_{\text{left}}^2}}}-{\sigma }_{\text{grad,left}}\right)}^2+\frac{1}{2}{\left(\sqrt{\frac{1}{N}{\displaystyle \sum \frac{{\left(\nabla x_{\text{right}}-\nabla y_2\right)}^2}{\nabla x_{\text{right}}^2}}}-{\sigma }_{\text{grad,right}}\right)}^2$ (9)

其中${\sigma }_{\text{grad,left}}$ ，${\sigma }_{\text{grad,right}}$ 分别为左雷达和右雷达的梯度误差。

3) 精度表达式：

$P=1/\left(1+\frac{1}{2}{P}_1+\frac{1}{2}{P}_2\right)$ (10)

综合观测误差和梯度误差得出最终的精度表达式，精度越高说明估测结果越准确。

4. 构建雷达–雨量计变分同化数学模型

4.1. 模型的建立

雨量计能够提供高精度的点雨量数据，但空间覆盖率有限，雷达能够实时探测一定区域内的降水情况，空间分辨率较高，但数据精度可能受到多种因素的影响。为了解决单一观测手段的缺陷，本文基于双雷达变分同化数学模型构建了一个能够高效整合雷达和雨量计多源观测数据的变分同化模型，构建过程如下：

1) 数据预处理，通过对需要优化的降水场x进行抽样，来同化雨量计和雷达的空间分辨率，提升降水估测的精度和可靠性；

2) 构建背景误差协方差矩阵B，见双雷达变分同化数学模型构建过程；

3) 确定目标函数，假设目标函数中除雨量计外的项权重相同，限制权重之和为1且最多为两位小数，雷达–雨量计变分同化模型目标函数如公式(10)所示：

(11)

4) 自适应交替优化算法

自适应交替优化算法[12]通过交替更新背景场和权重系数来迭代逼近最优解，初始化降水场和权重，固定权重，优化降水场，比较不同权重下目标函数优化后的最终值，选择最小值的权重为最佳权重${\omega }^{\ast }$ ，计算公式如下所示：

(12)

其中x为降水场。

$J_{min}=J\left(x,{\omega }^{\ast }\right)$ (13)

5) 反演

使用matlab工具对最终数据进行反演，分析比较反演结果。

4.2. 模型的检验与分析

1) 双雷达检验

为了提高模型的泛化性，减少单次观测误差带来的偏差，需要进行观测误差调整，多次模拟观测，每次引入一个随机误差项，取这些模拟的平均值作为背景场。

$D_{\text{sim},i}=D\times \left(1+{ϵ}_i\right)$ (14)

其中$D$ 为雷达观测数据，$\sigma$ 为观测误差头${ϵ}_i$ 是服从正态分布$\mathcal{N}\left(0,\sigma \right)$ 的随机误差项。

为了更好地反映降水结构，需要进行梯度误差调整，让背景场很好地反映出降水的梯度特征。

$G_m=\frac{\partial x_b}{\partial m},G_n=\frac{\partial x_b}{\partial n}$ (15)

其中$m,n$ 分别为纬度索引和经度索引，$G_m,G_n$ 分别为观测误差调整后的背景场对$m,n$ 的偏导数。

${G^{\prime }}_m=G_m\times \left(1+{ϵ}_{grad,i}\right),{G^{\prime }}_n=G_n\times \left(1+{ϵ}_{grad,i}\right)$ (16)

其中${\sigma }_{grad}$ 为降水场结构的估计误差，${ϵ}_{grad,i}$ 是服从正态分布$\mathcal{N}\left(0,{\sigma }_{grad}\right)$ 的随机误差项，使用积分的方式对$m,n$ 方向上做梯度的修正：

${x^{\prime }}_{sim,i}\left(m,n\right)=\frac{1}{2}\left[B^{\prime }\left(m-1,n\right)+{G^{\prime }}_m\left(m,n\right)+B^{\prime }\left(m,n-1\right)+{G^{\prime }}_n\left(m,n\right)\right]$ (17)

其中${x^{\prime }}_{sim,i}\left(m,n\right)$ 为背景场中各个点经过梯度误差调整的数值，最后进行N次模拟取平均：

$x_b=\frac{1}{N}{\displaystyle \sum _{i=1}^N{x^{\prime }}_{\text{sim},i}}$ (18)

使用matlab画图得到对误差调整后的雷达观测数据进行反演，结果如图1所示。

Figure 1. Error-adjusted radar background field

图1. 误差调整后的雷达背景场

将所得到的数据代入双雷达变分同化数学模型，经变分同化模型优化后的降水场如图2所示：

Figure 2. Precipitation field after variational assimilation

图2. 变分同化后的降水场

依据背景1的数据，将同化后的降水场矩阵代入到精度表达式中，得到该降水场的精度为0.654007，即降水场与观测场的相关误差的均方根极为贴近雷达的真实误差的均方根，可以得出本文提出的变分同化模型可以很好地同化双雷达的数据，反演出较为真实的降水场。由图1和图2对比可得，图1边界相对较为粗糙，存在着较为明显的跳跃与突变，这在一定程度上反映了数据在空间连续性上的不足。然而，在图2中，经过同化技术处理后的降水场，其边界展现出了平滑性，这种平滑并非简单的模糊处理，而是在保持降水场内部结构特征完整性的基础上，实现了数据间的自然过渡与融合。

2) 雷达–雨量计检验

在双雷达数据的基础上增加雨量计观测数据，通过联立方程：

$\{\begin{array}{l}{\omega }_1+{\omega }_2+{\omega }_3+{\omega }_4=1\hfill \\ {\omega }_1={\omega }_2={\omega }_3\hfill \\ {\omega }_4>{\omega }_1\hfill \\ \omega 最多有2位小数\hfill \end{array}$ (19)

由公式(19)可得ω的集合，使用L-BFGS-B优化算法对不同的目标函数进行优化求解，得到使得目标函数最小的权重，将其代入雷达–雨量计变分同化数学模型，最终得出最优权重下的反演场如图所示：

Figure 3. Precipitation field assimilated by increasing rain gauges

图3. 增加雨量计同化后的降水场

由图2和图3对比可得，在增加了雨量计数据之后，优化后的降水场通过平缓的过渡带展现了从高强度降水(红色区域)到低强度或无降水(蓝色区域)的自然过渡，有效减少了预测中的不合理突变，提升了预测结果的可信度和实用性。

在深入分析降水预测模型的性能提升过程中，将经过同化处理的降水场数据代入精度评估表达式(10)中，计算结果显示，引入雨量计数据后的降水场精度提升至0.759577，相较于未采用雨量计数据时的降水场精度，提升了0.105570，实现了16.1%的增长。这一数据直观地反映了雷达融合雨量计观测数据在提升降水预测精准度方面的重要作用。

综上所述，将雨量计数据纳入降水预测模型的同化过程中，是提升模型预测精度、优化降水场空间分布特征的有效途径。这一做法不仅增强了模型对复杂降水现象的捕捉能力，也为相关领域的研究和应用提供了更加精准的数据支持。

5. 总结

为提升降水场估测精度，降低灾害发生风险，本文构建了双雷达变分同化模型以及雷达–雨量计变分同化模型，设计不同来源数据的权重系数选择模型，并进行反演实验比较，进一步验证了雷达雨量计变分同化模型对于降水估计准确率的优势。该模型能有效提高降水场估测精度，满足提高气象预报的准确性和提升灾害预警能力的要求，并为未来研究提供了参考价值。