20、气象雷达极化技术中的降水滴谱反演方法

气象雷达极化技术中的降水滴谱反演方法

1. 同时衰减与滴谱反演

1.1 反演原理与步骤

在气象雷达测量中，为了准确获取降水滴谱（DSD）信息，需要考虑信号的衰减。通过特定的积分方程来同时估计衰减和反演降水 DSD 及微物理参数。反演步骤如下：
1. 从最后一个距离门（r = rN）开始，该距离门内的衰减可忽略。利用公式 7.8 从 (Z_{DR}’(r_N)) 求解 (\Lambda(r_N))，再用公式 7.7 从 (Z_{H}’(r_N)) 得到 (N_0(r_N))。
2. 使用 (N_0(r_N)) 和 (\Lambda(r_N)) 计算最后一个距离门的衰减项（公式 7.7 和 7.8 中的最后一项），然后根据 (r_{N - 1}) 处的 (Z_{H}’(r_{N - 1})) 和 (Z_{DR}’(r_{N - 1})) 求解 (N_0(r_{N - 1})) 和 (\Lambda(r_{N - 1}))。
3. 重复步骤 2，直到初始距离门（r = r1）。

1.2 正反递归方法对比

原理上，正向和反向递归都可用于 DSD 反演。但实际应用中，正向递归由于误差和误差积累往往不稳定。而反向递归利用了 (PIA_H) 和 (PIA_{DP}) 约束，结果更稳定，因此在实际应用中通常采用反向递归。

1.3 模拟示例

假设一条射线的降水 DSD 恒定，(N_0 = 8000 m^{-3}mm^{-1})，(\Lambda = 2.0 mm^{-1})。计算得到真实的降水微物理参数和雷达本征变量：(R = 34.4 mm hr^{-1})，(D_m = 2.0 mm)，(Z_H = 48.8 dBZ)，(Z_{DR} = 2.60 dB)。然后计算 (PIA_H) 和 (PIA_{DP})，并在雷达变量中加入 (Z_H) 为 1 dB、(Z_{DR}) 为 0.2 dB 的误差，得到模拟的雷达测量值。模拟结果显示，反向递归能给出稳定的反演结果，而正向递归在距离 (r > 25 km)（(PIA_H > 30 dB)）时结果发散，无法得到有效反演。

这种同时衰减和 DSD 反演的方法，在考虑 DSD 变化的情况下，比单独的衰减校正更准确，尤其在存在 DSD 模型和测量误差时。双极化和双频率技术与该方法有相似的原理和方程。

2. 降水滴谱的统计反演

2.1 确定性反演的局限性

之前讨论的 DSD 确定性反演虽有改善定量降水估计（QPE）和微物理参数化的潜力，但存在诸多局限。它难以处理测量和模型误差，除了 DSD 模型约束（如 (\mu - \Lambda) 关系）外，无法利用 DSD 数据的先验统计信息，也不能量化反演性能。例如，当 (Z_{DR}) 测量值接近零或为负时，无法从相应方程求解 (\Lambda)，也就无法反演 DSD。此外，它不提供反演误差协方差信息，不利于优化利用反演结果。

2.2 贝叶斯统计反演方法

为克服这些局限，引入了基于贝叶斯理论的统计反演方法。该方法将降水状态（DSD）参数和极化雷达测量值视为由概率密度函数（PDF）表征的随机变量。设状态向量 (x) 表示 DSD 参数，测量向量 (y) 表示极化雷达数据（PRD）。反演的目标是找到在测量向量 (y) 后验条件下 (x) 的 PDF。根据贝叶斯定理，后验条件 PDF (p_{post}(x | y)) 为：
[
p_{post}(x | y) = \frac{p_f(y | x) \cdot p_{pr}(x)}{\int p_f(y | x) \cdot p_{pr}(x) dx}
]
其中，(p_{pr}(x)) 是状态 (x) 的先验 PDF，(p_f(x | y)) 是给定状态 (x) 时观测 (y) 的正向条件 PDF。找到状态向量的后验条件 PDF 后，通过对状态 (x) 的整个范围积分计算期望值 (\bar{x}) 和标准差 (\sigma_x)：
[
\bar{x} = \frac{\int x \cdot p_f(y | x) \cdot p_{pr}(x) dx}{\int p_f(y | x) \cdot p_{pr}(x) dx}
]
[
\sigma_x^2 = \frac{\int (x - \bar{x})^2 \cdot p_f(y | x) \cdot p_{pr}(x) dx}{\int p_f(y | x) \cdot p_{pr}(x) dx}
]

2.3 DSD 参数的先验分布

为了获得后验条件 PDF，需要先确定状态向量（DSD 参数）的先验 PDF 和给定状态向量时观测向量（PRD）的正向条件 PDF。通过将 2DVD 测量的 DSD 拟合到伽马分布模型（使用截断矩拟合方法，利用 DSD 的二阶、四阶和六阶矩），得到 DSD 参数 (N_0) 和 (\Lambda) 的分布。这些分布有很大的偏态、大动态范围和与物理参数的非线性关系。为减少动态范围和缓解非线性效应，将 DSD 参数转换为 (N_0’ = \log_{10} N_0) 和 (\Lambda’ = \Lambda^{0.25}) 用于贝叶斯反演。转换后，轻雨时 (N_0’) 和 (\Lambda’) 的动态范围显著减小，接近高斯分布，能更好地表示中/大雨情况。多数 DSD 的 (N_0’) 值在 3 到 5 之间（即 (N_0) 约为 (10^3 - 10^5 m^{-3} mm^{-1})），(\Lambda’) 值在 1.1 到 1.6 之间（即 (\Lambda) 约为 1.5 - 6）。将离散化后的 (N_0’) 和 (\Lambda’) 的联合 PDF 保存为查找表，便于后续计算。

2.4 正向条件分布

为便于贝叶斯反演，还需要正向条件 PDF (p_f(x | y)) 来连接状态变量和测量值。考虑到雷达估计的采样误差通常呈高斯分布，假设正向条件 PDF 遵循双变量高斯分布：
[
p_f(Z_H, Z_{DR} | N_0’, \Lambda’) = \frac{1}{2\pi \cdot \sigma_{Z_H} \cdot \sigma_{Z_{DR}} \sqrt{1 - \rho^2}} \times \exp\left{ -\frac{1}{2(1 - \rho^2)} \left[ \left(\frac{Z_H - \bar{Z} H}{\sigma {Z_H}}\right)^2 - 2\rho \left(\frac{Z_H - \bar{Z} H}{\sigma {Z_H}}\right) \left(\frac{Z_{DR} - \bar{Z} {DR}}{\sigma {Z_{DR}}}\right) + \left(\frac{Z_{DR} - \bar{Z} {DR}}{\sigma {Z_{DR}}}\right)^2 \right] \right}
]
其中，(\rho) 是 (Z_H) 和 (Z_{DR}) 之间的相关系数，通常很低，假设为 0。(\bar{Z} H) 和 (\bar{Z} {DR}) 可通过正向模型/算子从状态 ((N_0’, \Lambda’)) 计算得到。(\sigma_{Z_H}) 假设为常数 2 dB，(\sigma_{Z_{DR}}) 是 (Z_H) 和 (Z_{DR}) 的函数。在上下界内，(\sigma_{Z_{DR}}) 假设为常数 0.3 dB（主要是估计误差）；若观测的 PRD 超出该区域，认为 (\sigma_{Z_{DR}}) 受异常降水影响，误差大于 0.3 dB，具体表达式为：
[
\sigma_{Z_{DR}} =
\begin{cases}
0.3[1 + (Z_{DR} - Z_{DR}^{up})] & \text{above} \
0.3 & \text{within} \
0.3[1 + (Z_{DR}^{low} - Z_{DR})] & \text{below}
\end{cases}
]

2.5 结果与评估

通过 DSD 数据的模拟雷达测量对贝叶斯方法进行测试。根据测量的 DSD 模拟 (Z_H) 和 (Z_{DR})，应用贝叶斯反演公式得到 DSD 参数的反演值，进而计算降水积分变量的均值。结果表明，多数反演的 (R) 值与观测值非常接近，反演的 (D_m) 比 (R) 更分散，可能是由于 DSD 模型不准确，部分 DSD 不能很好地用伽马分布表示，或者在反演中使用的 (\mu - \Lambda) 关系的 C - G 模型不能准确描述。贝叶斯反演结合 C - G 模型表现良好，降水估计偏差小于 8%，标准差低于 18%。

将反演算法应用于实际的 PRD 数据（KOUN 雷达），与经验单极化和双极化 QPE 结果对比。贝叶斯反演得到的雨率与双极化经验 QPE 结果相似，但在强对流区域，经验单极化估计的雨率高于贝叶斯反演。同时，贝叶斯反演给出的 (\Lambda’) 和 (N_0’) 估计的标准差图像趋势相似，可作为反演置信度的指标。大标准差通常出现在杂波、融化的雪/冰雹或生物散射体污染的区域，导致降水反演的不确定性增大。

通过将反演的雨率和 1 小时降水累积量与原位测量值对比，验证贝叶斯反演结果。结果显示，贝叶斯反演和经验估计器都能给出合理结果，但贝叶斯方法在对流核心区域表现更好，尤其是雷达回波受冰雹污染、(Z_H) 和 (Z_{DR}) 测量值极大的情况。具体的偏差和均方根误差（RMSE）见表 1：

反演类型	CHIC	MINC	NINN	SHAW	SPEN	WASH	KFFL
贝叶斯 - 偏差（%）	27.3	2.5	8.3	7.8	15.5	10.7	10.0
贝叶斯 - RMSE（%）	21.0	11.5	21.7	14.7	19.6	16.4	15.5
经验 (R(Z_H, Z_{DR})) - 偏差（%）	1.1	1.6	12.3	22.4	28.0	19.8	20.2
经验 (R(Z_H, Z_{DR})) - RMSE（%）	17.3	26.4	33.5	29.9	33.9	34.4	33.3
经验 (R(Z_H)) - 偏差（%）	35.4	48.1	48.2	49.5	44.9	48.9	50.4
经验 (R(Z_H)) - RMSE（%）	80.4	103.5	106.8	72.3	56.0	95.3	92.5

贝叶斯方法具有以下优点：
1. 同时获得估计值和估计误差，可评估反演的可靠性。
2. 考虑并包含测量和模型误差。
3. 利用降水微物理的先验（历史）信息。
4. 无需固定的经验估计器。

2.6 贝叶斯反演流程

graph LR
    A[确定DSD参数先验分布] --> B[模拟雷达测量值]
    B --> C[计算正向条件PDF]
    C --> D[应用贝叶斯公式反演DSD参数]
    D --> E[计算降水积分变量均值]
    E --> F[与原位测量值对比评估]

3. 变分反演

3.1 现有方法的不足

前面讨论的同时衰减校正和 DSD 反演是确定性方法，未考虑误差影响，可能无法最优利用多参数 PRD，尽管在衰减校正中包含了 DSD 信息。贝叶斯 DSD 反演方法是统计方法，考虑了误差影响，但未考虑衰减和空间天气信息，且不够灵活，难以添加更多观测/约束。

3.2 变分反演的原理与分类

变分方法结合了贝叶斯和同时反演方法的优点，能最优利用所有可用的信息/测量值。仅使用观测值的变分反演称为基于观测的反演，而在分析中包含数值天气预报（NWP）模型预报结果的称为数据同化（DA）方法。

3.3 变分反演的一般公式

变分反演也是一种统计反演，且考虑了状态变量的空间协方差。基于高斯分布误差模型的假设，从贝叶斯理论推导出变分方法的公式。假设状态向量估计 (x) 服从高斯分布，协方差矩阵为 (B)，观测向量 (y) 也服从高斯分布，误差协方差为 (R)，且分母为常数，则有：
[
-2\ln[p(x | y)] = [y - H(x)]^t R^{-1} [y - H(x)] + [x - x_b]^t B^{-1} [x - x_b] + C
]
其中，(x_b) 是状态向量的平均背景或初始猜测，(H(\cdots)) 是观测算子。最优反演是找到使概率 (p(x | y)) 最大的状态向量，等价于最小化代价函数：
[
J = [x - x_b] B^{-1} [x - x_b] + [y - H(x)]^t R^{-1} [y - H(x)]
]
得到最优解：
[
\hat{x} = x_b + [B^{-1} + H’^t R^{-1} H’]^{-1} H’^t R^{-1} [y - H(x_b)]
]
其中，(H’) 是雅可比算子，是包含正向观测算子 (H) 对状态向量各元素偏导数的矩阵。该变分公式可用于基于观测的反演或基于模型的分析（数据同化）。在数据同化中，(x_b) 表示模型背景；在基于观测的反演中，(x_b) 表示早期的分析迭代，(\hat{x}) 通过迭代求解。一旦定义了误差协方差和正向算子，就可以使用该公式找到状态向量的最优估计。

3.4 从 PRD 进行 DSD 反演的公式

在从极化雷达测量中进行 DSD 反演时，使用之前提到的 C - G 模型，选择 (N_0’ = \log_{10} N_0) 和 (\Lambda) 作为两个状态变量，即 (x = [N_0’, \Lambda]^t)。观测向量为 PRD：(y = [Z_H’, Z_{DR}’, K_{DP}]^t)。定义代价函数 (J(x)) 为：
[
J(x) = J_b(x) + J_{Z_H’}(x) + J_{Z_{DR}’}(x) + J_{K_{DP}}(x)
]
其中：
[
J_b(x) = \frac{1}{2} (x - x_b)^t B^{-1} (x - x_b)
]
[
J_{Z_H’}(x) = \frac{1}{2} [Z_H’ - H_{Z_H’}(x)]^t R_{Z_H’}^{-1} [Z_H’ - H_{Z_H’}(x)]
]
[
J_{Z_{DR}’}(x) = \frac{1}{2} [Z_{DR}’ - H_{Z_{DR}’}(x)]^t R_{Z_{DR}’}^{-1} [Z_{DR}’ - H_{Z_{DR}’}(x)]
]
[
J_{K_{DP}}(x) = \frac{1}{2} [K_{DP} - H_{K_{DP}}(x)]^t R_{K_{DP}}^{-1} [K_{DP} - H_{K_{DP}}(x)]
]
(J_b) 表示背景的贡献，(J_{Z_H’})、(J_{Z_{DR}’}) 和 (J_{K_{DP}}) 分别对应 (Z_H’)、(Z_{DR}’) 和 (K_{DP}) 的测量。

3.5 变分反演流程

graph LR
    A[定义状态向量和观测向量] --> B[确定误差协方差和正向算子]
    B --> C[初始化状态向量背景值]
    C --> D[计算代价函数]
    D --> E[迭代求解最优状态向量]
    E --> F[得到DSD反演结果]

变分反演方法为从 S - 、C - 和 X - 波段的 PRD（(Z_H)、(Z_{DR}) 和 (\Phi_{DP}/K_{DP})）进行二维雨 DSD 反演提供了一种有效的途径，展示了其在气象雷达极化技术中的适用性。通过结合贝叶斯和同时反演方法的优点，变分反演能够更好地处理误差、利用空间信息和添加更多约束，有望在降水微物理参数反演和定量降水估计中发挥重要作用。未来可以进一步研究如何优化变分反演的参数设置，提高反演的精度和稳定性，以及探索在更复杂气象条件下的应用。

4. 数据同化方法在降水滴谱反演中的应用

4.1 数据同化的概念与作用

数据同化（DA）方法是在变分反演的基础上，将数值天气预报（NWP）模型预报结果纳入分析过程。它的主要作用是融合观测数据和模型预报信息，以提高对降水滴谱（DSD）等气象参数的反演精度。与仅基于观测的变分反演相比，数据同化能够利用模型的动态信息，更好地描述大气的演变过程，从而更准确地反演降水微物理参数。

4.2 数据同化的流程

数据同化的基本流程如下：
1. 模型预报 ：使用 NWP 模型进行气象要素的预报，得到预报场。
2. 观测数据收集 ：收集气象雷达等设备的观测数据，如 (Z_H)、(Z_{DR})、(\Phi_{DP}/K_{DP}) 等。
3. 分析步骤 ：将观测数据与模型预报场进行融合，通过变分反演等方法调整模型状态，使其更符合观测数据。这一步骤通常涉及最小化代价函数，以找到最优的状态向量。
4. 更新模型状态 ：根据分析结果更新模型的初始条件，然后进行下一轮的预报和分析。

4.3 数据同化在 DSD 反演中的应用示例

以 S - 、C - 和 X - 波段的极化雷达数据为例，数据同化方法可以将这些观测数据与 NWP 模型的预报结果相结合，反演降水 DSD 参数。在实际应用中，可以按照以下步骤进行：
1. 初始化 ：设置 NWP 模型的初始条件和参数，同时确定观测数据的误差协方差和正向算子。
2. 模型预报 ：运行 NWP 模型，得到预报的气象要素场。
3. 观测数据处理 ：对雷达观测数据进行质量控制和预处理，去除噪声和异常值。
4. 分析步骤 ：使用变分反演公式（如第 3 节中提到的公式），将观测数据与模型预报场进行融合，调整模型状态。
5. 更新模型状态 ：根据分析结果更新 NWP 模型的初始条件，然后进行下一轮的预报和分析。

4.4 数据同化的优势与挑战

数据同化方法的优势在于能够综合利用观测数据和模型信息，提高反演的精度和稳定性。它可以考虑到大气的动态变化和空间相关性，从而更准确地描述降水微物理过程。然而，数据同化也面临一些挑战，例如：
1. 模型误差 ：NWP 模型本身存在一定的误差，这些误差可能会影响数据同化的结果。
2. 观测误差 ：雷达观测数据也存在误差，如何准确地估计和处理这些误差是数据同化的关键问题之一。
3. 计算成本 ：数据同化通常需要进行大量的计算，尤其是在处理高分辨率的观测数据和复杂的模型时，计算成本会显著增加。

优势	挑战
综合利用观测和模型信息	模型误差
提高反演精度和稳定性	观测误差
考虑大气动态变化和空间相关性	计算成本高

4.5 数据同化流程示意图

graph LR
    A[模型预报] --> B[观测数据收集]
    B --> C[分析步骤]
    C --> D[更新模型状态]
    D --> A

5. 不同反演方法的比较与总结

5.1 不同反演方法的特点

反演方法	原理	优点	缺点
同时衰减与滴谱反演	通过积分方程同时估计衰减和反演 DSD，利用反向递归提高稳定性	考虑 DSD 变化，比单独衰减校正更准确	未考虑误差影响，正向递归不稳定
贝叶斯统计反演	基于贝叶斯理论，将 DSD 参数和雷达测量视为随机变量	考虑测量和模型误差，利用先验信息，可评估反演可靠性	未考虑衰减和空间天气信息，不够灵活
变分反演	结合贝叶斯和同时反演方法，考虑状态变量的空间协方差	能最优利用所有可用信息，可添加更多约束	计算相对复杂
数据同化	将 NWP 模型预报结果与观测数据融合	综合利用观测和模型信息，考虑大气动态变化	存在模型误差和观测误差，计算成本高

5.2 不同方法的适用场景

• 同时衰减与滴谱反演 ：适用于对 DSD 变化较为关注，且衰减影响较大的情况，但对误差要求不高的场景。
• 贝叶斯统计反演 ：当测量和模型误差较大，且有一定的先验信息可用时，该方法能提供较为可靠的反演结果。
• 变分反演 ：在需要综合利用多参数 PRD，且希望考虑空间相关性和添加更多约束的情况下，变分反演是一个较好的选择。
• 数据同化 ：对于需要考虑大气动态变化，且有 NWP 模型支持的场景，数据同化方法能够更准确地反演降水微物理参数。

5.3 总结与展望

气象雷达极化技术中的降水滴谱反演方法各有优缺点，在实际应用中需要根据具体情况选择合适的方法。同时，不同方法之间也可以相互结合，以充分发挥各自的优势。例如，在数据同化中可以采用变分反演的方法进行分析步骤，或者在贝叶斯反演中考虑衰减的影响。

未来，随着气象雷达技术的不断发展和观测数据质量的提高，降水滴谱反演方法有望进一步优化。例如，可以研究更精确的误差模型，提高反演的精度和稳定性；探索在更复杂气象条件下的应用，如强对流天气、台风等；以及开发更高效的计算算法，降低计算成本。通过不断的研究和改进，降水滴谱反演方法将在气象预报、灾害预警等领域发挥更加重要的作用。

5.4 不同反演方法选择流程

graph LR
    A[关注DSD变化和衰减影响?] -->|是| B[同时衰减与滴谱反演]
    A -->|否| C[测量和模型误差大且有先验信息?]
    C -->|是| D[贝叶斯统计反演]
    C -->|否| E[需要综合利用多参数和空间信息?]
    E -->|是| F[变分反演]
    E -->|否| G[需要考虑大气动态变化和有NWP模型支持?]
    G -->|是| H[数据同化]
    G -->|否| I[重新评估需求]