21、天气雷达极化测量中的变分反演与数据同化方法

天气雷达极化测量中的变分反演与数据同化方法

在天气雷达极化测量领域，为了更准确地获取降水的微观物理信息，变分反演和数据同化（DA）方法发挥着重要作用。本文将详细介绍这些方法的原理、操作步骤以及应用情况，并分析其面临的挑战和潜在的发展方向。

1. 背景误差协方差矩阵与新状态变量

在天气雷达数据处理中，背景误差协方差矩阵 (B) 是一个 (m\times m) 的矩阵，其中 (m) 是状态向量 (x) 的大小，等于二维分析区域内的网格数量乘以状态参数的数量。为了避免对 (B) 这样的大矩阵进行求逆运算，引入了新的状态变量 (v)，其表达式为：
[v = D^{-1}\Delta x]
其中 (\Delta x = x - x_b)，且 (DD^T = B)。这里的 (\Delta) 表示增量，(D) 是矩阵 (B) 的平方根。通过引入这个新的状态变量，可以简化后续的计算。

2. 代价函数及其最小化

代价函数 (J) 经过改写后，其最小化可以通过代价函数梯度 (\nabla_v J) 来实现。代价函数 (J) 的表达式较为复杂，这里不再详细列出。而代价函数梯度 (\nabla_v J) 的表达式为：
[
\nabla_v J =
\left(D^T H_{H}^T R_{H}^{-1} H_{H} Dv + d_{H}\right) +
\left(D^T H_{DR}^T R_{DR}^{-1} H_{DR} Dv + d_{DR}\right) +
\left(D^T H_{DP}^T R_{DP}^{-1} H_{DP} Dv + d_{DP}\right)
]
其中 (d) 是观测的创新向量，即 (d = y - H(x_b))。通过不断调整状态变量 (v)，使得代价函数 (J) 达到最小，从而得到最优的状态向量 (x)。

3. 背景误差协方差矩阵的空间影响

背景误差协方差矩阵 (B) 决定了观测的空间影响。矩阵 (B) 的元素 (b_{ij}) 可以用高斯相关模型来建模：
[b_{ij} = \sigma_b^2 \exp\left(-\frac{r_{ij}^2}{2r_L^2}\right)]
其中 (i) 和 (j) 表示分析空间中的两个网格点，(\sigma_b^2) 是背景误差协方差，(r_{ij}) 是第 (i) 个和第 (j) 个网格点之间的距离，(r_L) 是背景误差的空间去相关长度，在风暴尺度雷达数据分析中通常假设为常数（2 - 4 km）。

4. 前向观测算子与迭代过程

4.1 前向观测算子

给定每个网格点的两个雨滴谱分布（DSD）参数 ([N_0’, \Lambda])，可以确定雨滴的 DSD，并使用 T 矩阵方法从预先计算的散射振幅值计算雷达固有变量 (Z_H)、(Z_{DR}) 和 (K_{DP})，以及水平极化的特定衰减 (A_H) 和特定差分衰减 (A_{DP})。前向算子 (Z_H’) 和 (Z_{DR}’) 在每个距离门的表达式分别为：
[\sum_{i = 1}^{n - 1} Z_H’(n) = Z_H(n) - 2A_H(i)\Delta r]
[\sum_{i = 1}^{n - 1} Z_{DR}’(n) = Z_{DR}(n) - 2A_{DP}(i)\Delta r]
其中 (i) 和 (n) 分别表示从雷达位置起的第 (i) 个和第 (n) 个距离门，(\Delta r) 是距离分辨率。

4.2 偏导数计算与查找表

为了计算前向算子 (H) 的偏导数，需要计算每个极化测量变量（(Z_H’)、(Z_{DR}’) 和 (K_{DP})）相对于两个状态变量（(\Lambda) 或 (N_0’)）的偏导数，并将这些偏导数存储在 10 个查找表中。在每个查找表中，参数 (\Lambda) 的取值范围从 0 到 50，参数 (N_0’) 的取值范围从 0 到 15，为了确保足够的精度，每个参数的范围以 0.02 的间隔进行离散化。因此，每个查找表有 (2501\times751) 个元素。通过这些查找表，可以找到任意给定的 (\Lambda) 和 (N_0’) 值对应的前向算子 (H) 的偏导数。对于查找表值之间的 (\Lambda) 或 (N_0’) 值，可以进行插值以进一步提高精度。

4.3 迭代过程步骤

为了最小化代价函数 (J)，迭代过程如下：
1. 输入数据 ：输入雷达数据 (Z_H’)、(Z_{DR}’) 和 (K_{DP})、查找表以及背景状态参数。为了进行数据质量控制，仅使用分析区域中信噪比（SNR）> 1 dB 的雷达测量值。
2. 初始化状态向量 ：变分反演的初始状态向量等于背景状态向量（即 (x = x_b)），迭代从 (v = 0) 开始。使用前向算子和散射振幅的查找表计算每个网格点的雷达变量 (Z_H)、(Z_{DR})、(K_{DP})、(A_H) 和 (A_{DP})。将这些雷达变量插值到观测点，并加上衰减以得到 (Z_H’)、(Z_{DR}’) 和 (K_{DP}) 以及代价函数。
3. 更新状态向量 ：根据观测的创新向量修改状态向量，然后重复计算雷达变量以进行下一次迭代。在搜索代价函数的最小梯度过程中，不断更新状态向量，直到迭代收敛。
4. 获得分析结果 ：当最小化过程收敛后，获得 DSD 参数的分析场。需要注意的是，第一次收敛的分析结果可能并不令人满意。为了提高反演精度，可以将第一次收敛的分析结果作为新的背景，重复分析（即迭代）过程。这种将前一次分析结果作为新分析的背景的重复过程被视为“外循环”迭代。通常，经过几次外循环迭代后，可以得到代价函数相对较小的满意分析结果。

5. 变分反演算法的应用

变分反演算法首先在模拟雷达数据上进行了测试和评估，然后应用于实际雷达数据。在模拟数据的情况下，使用俄克拉荷马州诺曼市的 S 波段极化 KOUN 雷达的实际数据模拟了 X 波段 CASA 极化雷达的测量值。具体步骤如下：
1. 模拟数据生成 ：假设 S 波段雷达的 (Z_H) 和 (Z_{DR}) 测量值不受降水衰减影响，使用第 7.2 节讨论的方法反演雨滴的 DSD。然后使用反演得到的 DSD 计算 X 波段雷达变量，将其视为“真实值”。通过公式（7.26）和（7.27）获得衰减后的 X 波段 (Z_H’) 和 (Z_{DR}’)。最后，在衰减后的 (Z_H’) 和 (Z_{DR}’) 上添加偏差和随机误差，以模拟 X 波段雷达观测值，用于测试变分反演算法。
2. 实验设计 ：总共设计了 12 个实验，所有这些实验都包含 DSD 模型误差，即模拟的真实值假设为指数 DSD，而反演假设为 C - G DSD。在这些测试中应用了一个恒定的背景 ((N_0’ = 3, \Lambda = 5))。测试 1 到 4 假设模拟观测值无偏差但有不同的随机误差；测试 5 到 8 假设随机误差与测试 1 相同，但有不同的偏差；测试 9 到 12 假设偏差与测试 5 到 8 相同，但随机误差不同。

实验结果总结在表 1 中：
| 测试编号 | (Z_H) 反演偏差/模拟偏差 (dBZ) | (Z_{DR}) 反演偏差/模拟偏差 (dB) | (K_{DP}) 反演偏差/模拟偏差 (°/km) | (Z_H) 反演均方根误差/模拟误差 (dBZ) | (Z_{DR}) 反演均方根误差/模拟误差 (dB) | (K_{DP}) 反演均方根误差/模拟误差 (°/km) |
| — | — | — | — | — | — | — |
| 1 | 0.091/0 | 0.027/0 | 0.004/0 | 0.393/0.5 | 0.107/0.1 | 0.084/0.1 |
| 2 | 0.083/0 | 0.009/0 | 0.006/0 | 0.409/1.0 | 0.108/0.2 | 0.083/0.2 |
| 3 | 0.178/0 | 0.023/0 | 0.012/0 | 0.476/1.5 | 0.110/0.3 | 0.088/0.3 |
| 4 | 0.267/0 | 0.036/0 | 0.019/0 | 0.537/2.0 | 0.120/0.4 | 0.093/0.4 |
| 5 | 0.440/0.125 | 0.115/0.025 | 0.020/0.025 | 0.597/0.5 | 0.159/0.1 | 0.084/0.1 |
| 6 | 0.841/0.25 | 0.219/0.05 | 0.037/0.05 | 0.952/0.5 | 0.253/0.1 | 0.092/0.1 |
| 7 | 1.575/0.5 | 0.411/0.1 | 0.067/0.1 | 1.687/0.5 | 0.445/0.1 | 0.114/0.1 |
| 8 | 2.879/1.0 | 0.755/0.2 | 0.118/0.2 | 3.037/0.5 | 0.807/0.1 | 0.160/0.1 |
| 9 | 0.448/0.125 | 0.113/0.025 | 0.022/0.025 | 0.606/0.75 | 0.157/0.15 | 0.085/0.15 |
| 10 | 0.862/0.25 | 0.216/0.05 | 0.040/0.05 | 0.979/1.0 | 0.250/0.2 | 0.095/0.2 |
| 11 | 1.604/0.5 | 0.408/0.1 | 0.071/0.1 | 1.724/1.25 | 0.443/0.25 | 0.117/0.25 |
| 12 | 2.940/1.0 | 0.747/0.2 | 0.122/0.2 | 3.117/1.5 | 0.801/0.3 | 0.165/0.3 |

从实验结果可以看出：
- 在测试 1 到 4 中，反演的均方根误差值通常小于真实的均方根误差值，这表明该算法可以平滑观测误差，使得最终分析结果的误差小于所有使用信息源的误差，符合最优估计理论。同时偏差随着均方根误差的增加而增加，但总体偏差非常小，这与测量中没有系统偏差的事实相符。
- 测试 5 到 8 与测试 1 相比，除了包含不同的偏差外，随机误差相同。测试 5 到 8 的反演结果显示出明显的偏差和均方根误差，除了 (K_{DP}) 的一些值外，所有反演偏差或均方根误差都大于模拟测量偏差或误差。这表明变分算法对测量偏差比随机误差更敏感，测量偏差不仅会在反演中引入更大的偏差，还会增大反演的均方根误差。因此，在进行变分分析之前，应尽最大努力消除测量偏差。
- 测试 9 到 12 的测量偏差与测试 5 到 8 相同，但随机误差更大。然而，测试 12 的反演偏差和均方根误差与测试 8 几乎相同，再次证明了该算法对测量偏差比测量误差更敏感。

此外，在实际 X 波段雷达数据的应用中，对 40 km×40 km 的区域（401×401 个网格，网格间距为 100 m）进行了反演分析，假设 (Z_H’)、(Z_{DR}’) 和 (K_{DP}) 的观测误差分别为 2 dB、0.2 dB 和 0.2°/km，相关长度 (r_L) 假设为 2 km。结果表明，X 波段反演结果与 KOUN 雷达的观测结果吻合良好，同时能够捕捉到 (Z_H)、(Z_{DR}) 和 (K_{DP}) 中风暴结构的更多细节，验证了反演算法在实际 X 波段雷达数据应用中的有效性。

6. 数据同化方法的必要性与挑战

虽然变分方法考虑了测量误差、天气的空间相关性和衰减，但在考虑冰相、多种水凝物种类和/或多矩微物理时，用于反演云/降水物理的独立信息数量可能小于独立状态变量的数量。例如，WRF 双矩六物种（WDM6）方案有 12 个微物理状态变量，而 WSR - 88D 极化雷达数据（PRD）大多数时候只有四个测量值，最多提供三个独立变量。因此，仅依靠 PRD 无法确定所有 12 个状态变量，需要结合其他约束条件，如数值天气预报（NWP）模型的物理约束。

将 NWP 模型的预报结果纳入反演/分析的方法称为数据同化（DA）方法。然而，在实际应用中，实施 PRD DA 并不容易，其进展受到限制。主要原因包括：
1. 模型和参数差距 ：数值天气预报和天气雷达工程领域使用的理论、模型和参数存在差距。NWP 模型由一组动态、热力学和微物理方程组成，气象学家习惯处理雷达无法直接测量的状态变量，对极化雷达变量及其信息内容和测量误差结构可能不太熟悉。而雷达工程师则专注于开发天气测量和基于观测的反演技术。
2. 缺乏简单前向算子 ：极化雷达变量基于电磁波散射/传播和波统计定义，无法用建模者容易使用的状态变量的简单函数表示，缺乏简单的前向算子。
3. 相对误差大 ：PRD 的相对误差较大，这增加了数据处理和分析的难度。
4. 模型微物理参数化误差 ：模型微物理参数化存在较大的变化和误差，影响了反演的准确性。
5. 非线性关系复杂 ：模型状态变量和极化变量之间存在高度非线性关系，增加了数据同化的复杂性。

为了成功实施 PRD DA，需要有与 NWP 模型兼容的微物理参数化方案、准确快速的前向算子以及准确的误差表征。

7. 观测算子与误差分析

极化雷达变量表示为散射振幅的二阶矩，通过求解水凝物散射的麦克斯韦方程来计算。计算这些变量需要微物理信息（DSD/PSD、形状、取向和组成），但并非所有这些微物理属性都能在 NWP 模型模拟中预测。因此，为了便于从 NWP 模型输出计算雷达变量，需要做出许多假设，如 DSD 模型、形状大小关系、密度和取向等。

此外，在雷达分辨率体积中，按水混合比计算最占主导地位的物种不一定是极化特征的主要贡献者。例如，冰雹在雷达分辨率体积中的含水量可能仅占 <5%，但却主导了极化雷达的测量。同时，水 - 冰相变等复杂过程也增加了计算和分析的难度。

综上所述，变分反演和数据同化方法在天气雷达极化测量中具有重要的应用价值，但也面临着诸多挑战。未来的研究需要进一步解决这些挑战，提高反演和同化的准确性，为天气分析和预报提供更可靠的支持。例如，可以探索同时使用多频率 PRD 以获取更多独立信息，或者结合模型物理约束来进一步优化反演结果。此外，培养具有 NWP 和雷达极化测量背景和经验的科学家，开展深入研究以建立状态变量和雷达变量之间的准确关系，也是未来发展的重要方向。

天气雷达极化测量中的变分反演与数据同化方法

8. 数值天气预报与雷达工程的差异

数值天气预报（NWP）模型和雷达极化测量分别源于气象科学和雷达工程两个不同领域。以下是它们之间的详细差异：
- NWP 模型方面 ：NWP 模型由一系列动态、热力学和微物理方程构成。气象学家借助初始大气状态（如压力、温度、密度和云物理等）对这些方程进行数值求解，从而预测未来的天气状况。他们所处理的状态变量往往无法通过雷达直接测量，尽管有部分气象学文献对极化雷达变量、其信息内容以及测量误差结构有所涉及，但总体而言，NWP 气象学家对这些内容的熟悉程度有限。
- 雷达极化测量方面 ：极化雷达变量是基于电磁波散射/传播和波统计进行定义的，通过雷达脉冲信号进行估算。雷达工程师专注于研发用于气象测量和基于观测的反演技术，而气象学家则更关注观测、描述和解释极化特征，以开展微物理研究。

这种差异导致了在将雷达数据应用于 NWP 模型时存在一定的困难，需要既熟悉 NWP 又了解雷达极化测量的科学家进行深入研究，以建立状态变量和雷达变量之间的准确关系。

9. 数据同化方法的潜在解决方案

为了克服数据同化（DA）在实施 PRD DA 时面临的挑战，可以从以下几个方面寻找解决方案：
- 兼容的微物理参数化方案 ：开发与 NWP 模型相兼容的微物理参数化方案，确保模型能够准确描述云/降水的物理过程，同时能够与雷达数据进行有效的融合。
- 准确快速的前向算子 ：研究和开发准确且快速的前向算子，使得能够从 NWP 模型的输出中高效地计算出极化雷达变量。这可能需要结合先进的数值计算方法和物理模型。
- 准确的误差表征 ：对 PRD 的误差进行准确的表征，包括测量误差、模型误差以及它们的空间相关性。这有助于在数据同化过程中合理地利用雷达数据，减少误差的影响。
- 多源数据融合 ：除了 PRD 数据，还可以考虑融合其他数据源，如卫星数据、地面观测数据等，以增加独立信息的数量，提高数据同化的效果。
- 模型物理约束 ：充分利用 NWP 模型的物理约束，将其纳入数据同化过程中，以限制状态变量的取值范围，提高反演的稳定性和准确性。

10. 数据同化流程分析

数据同化（DA）的流程可以用以下 mermaid 流程图表示：

graph LR
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;

    A(初始条件):::process --> B(NWP 模型):::process
    B --> C(模型输出):::process
    D(雷达观测数据):::process --> E(前向算子):::process
    E --> F(模拟雷达变量):::process
    C --> G(数据同化系统):::process
    F --> G
    G --> H(分析场):::process
    H --> I(预报场):::process
    I --> B

该流程图展示了数据同化的主要步骤：
1. 初始条件 ：提供 NWP 模型所需的初始大气状态。
2. NWP 模型 ：根据初始条件进行数值模拟，得到模型输出。
3. 雷达观测数据 ：获取极化雷达的观测数据。
4. 前向算子 ：将模型输出转换为模拟雷达变量。
5. 数据同化系统 ：将模拟雷达变量与实际观测数据进行融合，得到分析场。
6. 分析场 ：作为下一次预报的初始条件，输入到 NWP 模型中。
7. 预报场 ：NWP 模型根据分析场进行预报，完成一个循环。

11. 总结与展望

天气雷达极化测量中的变分反演和数据同化方法在云/降水物理反演和天气预报中具有重要的应用价值。变分反演算法通过引入新的状态变量和迭代过程，能够有效地利用雷达数据进行雨滴谱分布参数的反演，并且在模拟数据和实际数据的测试中都取得了较好的效果。然而，该算法对测量偏差较为敏感，因此在应用前需要尽量消除测量偏差。

数据同化方法则为解决变分方法在处理复杂微物理问题时的局限性提供了一种途径。通过结合 NWP 模型的物理约束，可以更准确地反演云/降水的物理参数。但目前实施 PRD DA 面临着诸多挑战，包括模型和参数差距、缺乏简单前向算子、相对误差大、模型微物理参数化误差以及非线性关系复杂等。

未来的研究方向包括：
- 技术改进 ：进一步优化变分反演和数据同化算法，提高反演和同化的准确性和效率。例如，研究更有效的前向算子和误差表征方法，以更好地处理复杂的微物理过程。
- 多源数据融合 ：探索同时使用多频率 PRD 以及其他数据源（如卫星数据、地面观测数据等）的方法，以获取更多的独立信息，提高反演和预报的精度。
- 人才培养 ：培养具有 NWP 和雷达极化测量背景和经验的科学家，促进两个领域的交叉融合，推动天气雷达极化测量技术的发展。
- 应用拓展 ：将变分反演和数据同化方法应用到更广泛的气象领域，如强对流天气预警、定量降水估计等，为气象预报和灾害防御提供更有力的支持。

通过不断地研究和改进，天气雷达极化测量中的变分反演和数据同化方法有望在气象领域发挥更大的作用，为我们提供更准确的天气信息和更有效的灾害预警。