关于卡尔曼滤波器同化技术在曲靖气象雷达中获取风场信息的研究 孙星星， 蒋春华， 冯健， 刘毅， 吴瑜文， 徐彤， 乔建东， 邓中心， 周晨， 张玉强

《Remote Sensing》：Research on Kalman Filter Assimilation Technology for Wind Field Information from Qujing Meteor Radar Xingxin Sun, Chunhua Jiang, Jian Feng, Yi Liu, Yewen Wu, Tong Xu, Jiandong Qiao, Zhongxin Deng, Chen Zhou and Yuqiang Zhang

　　**亮点**
- 主要发现是什么？本研究关注80–100公里高度范围内的近空间区域的水平风场。利用曲靖站流星雷达的观测数据，基于卡尔曼滤波器的同化方法显著提高了HWM14背景模型的风场精度。

- 主要发现的意义是什么？通过结合卡尔曼滤波器同化技术，流星雷达观测数据被有效地整合到了HWM14背景场中。这种方法解决了由于流星在雷达测量中分布不均导致的数据缺失和误差问题。因此，同化后的风场分析误差大幅减少，特别是东西向和纬向风速的最大误差降低了50%以上。

**摘要**
全天候流星雷达被广泛用于观测近空间大气风场，这是近空间环境的一个关键参数。由于流星计数分布的时空不确定性，流星雷达在获取大气风场数据时可能会遇到误差和数据缺失。本研究利用卡尔曼滤波器同化技术与HWM14相结合，利用曲靖流星雷达（25.6°N, 103.7°E）的大气风场观测数据，研究了80–100公里高度范围内的水平风场。同化结果表明，经过卡尔曼滤波器同化后，HWM14的大气风场精度显著提高。同化风场分析与流星雷达风场分析值之间的差异显著减小：在90公里高度处，东西向风速的平均最大误差从同化前的14.7 m/s降低到12.0 m/s，提高了55.0%；在100公里高度处，纬向风速的平均最大误差从16.1 m/s降低到14.2 m/s，提高了53.4%。此外，同化风场分析与流星雷达风场分析值之间的偏差标准差也大幅减小。同化后的大气风场信息对于进一步研究近空间区域的大气扰动变化和动态过程具有重要意义。

**1. 引言**
近空间位于天空和外太空之间，大约在20–100公里的高度范围内，主要由平流层、中间层和下热层组成。重力波、潮汐波和行星波等常见于低层大气的波动可以通过耦合传播到这一区域。在此过程中，各种波动携带的动量和能量释放到近空间区域，导致该区域的大氣温度和背景风场发生剧烈变化[1,2,3]。温度和风场的变化影响该区域的大气结构，并进一步影响航空航天操作。因此，研究近空间环境变化以理解该区域的大气变化规律至关重要。

近年来，随着无线电技术和遥感技术的发展，基于地面的检测技术逐渐成为观测近空间大气环境的主要手段。地面检测技术主要包括激光雷达和各种无线电雷达（VHF雷达、MF雷达、流星雷达等）。其中，流星雷达由于其低成本、全天候运行和无人值守操作等优势，成为最常用的检测工具。流星雷达的发展和应用经历了多个阶段：早期工作主要集中在系统建设和区域观测验证上；全天空干涉技术的发展使检测精度和可靠性显著提升[4,5,6]。近年来，利用其高垂直分辨率的风场特性，这项技术在揭示大气波动动态[7,8]和建立关键大气参数（如密度和温度）的长期气候特征[9,10,11]方面取得了重要突破。这些进展确立了流星雷达在描述大气波动活动和传输过程中的不可或缺的观测作用，为深入理解中间层和下热层的大气动态提供了关键的数据基础。目前，全球有数十个单站全天空流星雷达在运行，部分雷达的数据已经公开发布。例如，美国麻省理工学院（MIT）发布了Buckland Park站的观测数据[6]，日本京都大学发布了赤道附近的Kototabang站的数据[12]，中国子午线项目数据中心提供了莫赫和武汉站的观测数据[13,14]。这些数据集对研究全球大气波动（如行星波和潮汐波）及其传播[15,16,17,18]以及下层和上层大气之间的耦合[19,20]做出了重要贡献，支持了近空间相关科学和应用研究。不过，这些流星雷达数据集提供的风场数据的时间分辨率为0.5–2.0小时，空间（主要是高度）分辨率为2–3公里。然而，单站全天空流星雷达获得的水平平均风场仅限于站点所在区域，缺乏水平分辨率[21]，无法获得近空间更大范围内的风场信息。因此，这一高度的实时风场数据无法满足大范围区域（国家或地区）的环境操作需求。

水平风模型（HWM）是一个全球参考大气风场模型。该模型利用各种观测和检测数据（如卫星、探空火箭和雷达），通过球谐拟合在地面以上0至500公里的高度范围内获取风场数据[22,23]。HWM（主要是HWM07和HWM14版本）在还原大气真实动态过程方面具有显著优势。江等人将气象火箭测量的风场与HWM07的结果进行了比较，发现两者趋势基本一致[24]。通过对比Fabry-Perot干涉仪数据和HWM14的数据，江等人发现HWM14能够很好地再现夜间热层风的变化特征[25]。然而，作为全球参考大气风场模型，HWM提供的风场结果是全局性的，与特定地区（如中国西南部）的实时风场结果存在显著差异。李等人分析了昆明站ST雷达的风场观测结果，发现该区域测量的风场与HWM07的结果存在差异[26]。周宝珠等人将十个全球站的流星雷达观测数据与HWM07模拟的数据进行了比较，发现HWM07的近空间风场结构与流星雷达观测结果一致，但在冬季两者之间存在显著差异[27]。

总之，基于地面的全天空流星雷达观测设备可以在站点上方80–100公里的高度范围内获取风场数据，具有良好的实时性能和可靠性[6]。相比之下，HWM模拟的全球风场数据表现出高度规律性。考虑到流星计数分布的时空不确定性（由于流星体的不确定性和难以准确检测微弱流星信号），本研究在80至100公里的高度范围内探讨了卡尔曼滤波器大气风场同化技术。它结合了单站流星雷达观测数据和HWM14数据，基于曲靖站的流星雷达大气风场观测数据和HWM14，分析了同化风场中的误差。该研究解决了全天空流星雷达观测数据在时间和空间上分布稀疏的风场同化问题，开发了定制的误差模型和局部协方差方案，有效校正了HWM14经验模型在中国地区的系统偏差。结果表明，单站同化得到的分析结果比背景场更可靠，为操作风场反演提供了一种经过验证的低成本方法。

**2. 数据来源**
2008年，中国无线电波传播研究所于云南省曲靖市（25.6°N, 103.7°E）建立了一个全天空流星雷达系统，主要用于近空间大气风场的观测和研究[28,29,30,31]。曲靖站的流星雷达系统是ATRAD Pty Ltd.（位于澳大利亚Thebarton）的成熟产品，在整体结构设计和内部配置上与Buckland Park全天空流星雷达非常相似[5,6]。该系统包括发射机、接收机、信号处理子系统、数据处理子系统、天线子系统等。详细技术参数见表1[28]供参考。天线子系统采用多天线空间布局，包括一个双极交叉偶极子发射天线和五个双极半波偶极子交叉天线组成的干涉接收天线阵列，接收天线呈“十字”排列。数据收集和分析由两台独立的计算机完成。本研究使用曲靖站2022年的流星雷达风场观测数据，覆盖80–100公里的高度范围。近空间风场的背景场数据来自HWM14。

**3. 基于卡尔曼滤波器的同化方法**
**3.1. 算法原理**
根据卡尔曼滤波器理论，观测数据的同化过程可以分为两部分：更新和预测[32]。更新部分负责根据当前时间的观测数据校正背景模型的模拟和预测结果，其表达式如下：
$$
\mathbf{X}_a = \mathbf{X}_b + \mathbf{K}(\mathbf{Y} - \mathbf{H}_2\mathbf{X}_b)
$$
其中，K是增益矩阵，用于利用观测数据调整模型生成的背景场；$\mathbf{X}_a$是同化后的分析值，HWM14可以获取风场背景场$\mathbf{X}_b$；$\mathbf{Y}$表示流星雷达设备观测的状态值。增益矩阵K可以通过背景场误差协方差矩阵$\mathbf{B}$和观测场误差协方差矩阵$\mathbf{R}$根据以下公式求解：
$$
\mathbf{K} = \mathbf{B}^\mathbf{H}^2(\mathbf{H}^2\mathbf{B}\mathbf{H}^2+\mathbf{R})^{-1}
$$
$\mathbf{B}_a$是同化后分析场的协方差矩阵，可以通过以下公式求解：
$$
\mathbf{B}_a = \mathbf{B} - \mathbf{B}\mathbf{H}^2(\mathbf{H}^2\mathbf{B}\mathbf{H}^2+\mathbf{R})^{-1}\mathbf{H}^2\mathbf{B}
$$
在近空间风场同化过程中，基于当前时间分析场值预测下一个时间间隔的近空间风场值定义为该过程的预测部分。预测值作为下一个时间步的背景场，预测误差作为下一个时间的背景误差，其表达式如下：
$$
\mathbf{B}_{a, t+1} = \mathbf{B}_a - \mathbf{B}\mathbf{H}^2(\mathbf{H}^2\mathbf{B}\mathbf{H}^2+\mathbf{R})^{-1}\mathbf{H}^2\mathbf{B}
$$
$$
\mathbf{X}_{a, t+1} = \mathbf{M}\mathbf{X}_a\mathbf{t} + \mathbf{\omega}_t
$$
其中，$\mathbf{M}$是状态转移矩阵，$\mathbf{\omega}_t$是状态转移误差；这两个函数用于将当前时间的分析场转换为下一个时间点的背景场。

**3.2. 数据预处理**
在对近空间风场进行同化之前，需要预处理观测数据以减少异常值和噪声，从而提高数据质量[33]。本研究采用双权估计方法进行数据质量控制[34]。为数据点分配不同的权重可以减轻异常值对整体估计的影响。图1展示了数据质量控制前后东西向风场结果的对比。左侧面板显示了2022年1月1日的原始水平风速剖面（蓝色散点），实线红色曲线表示使用上述双权估计方法计算出的动态阈值。可以看出，原始数据中包含许多超出此阈值的散点（异常值），尤其是在80公里高度附近。这些值可能来自检测噪声或间歇性干扰。右侧面板显示了通过基于阈值的过滤有效去除了这些异常值。由此产生的风场剖面显示出一个更加集中且符合物理规律的分布，与中间大气的典型风结构特征更为吻合。这一比较表明，双权重估计方法能够自动且有效地识别并剔除风数据中的异常观测值。因此，它为后续的同化分析提供了更高质量、更可靠的数据基础。图1. 数据质量控制前（左）和质量控制后（右）的纬向风场结果。

3.3 协方差模型构建
协方差模型由背景误差协方差矩阵、观测误差协方差矩阵和观测算子组成，它反映了数据之间的相关性。本研究中描述的近空间风场同化过程指的是将气象雷达观测数据与背景模型进行整合的过程。因此，误差协方差模型包括了观测数据的误差协方差和背景模型的误差协方差。观测算子表示观测误差协方差与背景模型误差协方差的相对大小和空间分布，这决定了观测场和背景场对同化结果的相对影响权重和空间结构。

3.3.1 观测场误差协方差矩阵的构建
由于风场观测是独立事件，本研究中使用对角矩阵来表示观测场误差协方差矩阵：
$$
R_{ij} = \begin{cases}
\delta_0 & (i = j) \\
0 & (i \neq j)
\end{cases}
$$
公式中，非对角线元素的值为0，$\delta_0$是权重参数。$\delta_0$的大小决定了观测值在同化系统中的权重。较大的$\delta_0$表示较大的观测场误差协方差，较小的$\delta_0$表示较小的误差协方差。在近空间风场同化系统中，观测结果比HWM14更加符合实际情况；因此需要给$\delta_0$设定一个较小的值。在本研究中，$\delta_0 = 0.01$。

3.3.2 背景场误差协方差矩阵的构建
由于观测数据彼此独立，风场的空间相关性信息只能从背景场误差协方差矩阵中得出。背景场误差协方差矩阵是使用基于高斯理论的协方差模型来建立的[35]。近空间风场的空间相关性信息可以细分为水平相关性和垂直相关性。背景场误差协方差的具体计算如下：
$$
B_{ij} = \sigma_{ij} \sigma_{ij} \lambda_{ij}
$$
公式中，$B_{ij}$是背景场之间的误差协方差值。
$$
\sigma_{ij} = \delta_{bj} x_{bj} x_{bj}
$$
$\delta_{bj}$是一个超参数，取值$\delta_{bj} = 20\delta_0$。
$$
\lambda_{ij} = \frac{1 - \rho^2 x_{ij}^2}{2\ell^2 x_i^2 + 2\ell^2 x_j^2 - \rho^2 \ell^2}
$$
$\lambda_{ij}$是相关系数，$\ell_i^2$是水平相关长度，$\ell_j^2$是垂直相关长度，$x_{ij}$表示网格点之间的水平距离，$\ell_j$表示网格点之间的垂直距离。
对于传统的高斯模型，相关长度是一个固定值，不随时间、经度或纬度变化。即认为近空间风场是各向同性的，或者在各向经度、纬度和垂直方向上都是各向同性的，这通常不符合近空间风场的实际特性。一般来说，近空间风场表现出明显的日变和季节性变化。结合统计协方差模型的思想，根据近空间风场的变化规律，可以利用大量历史数据来近似获得满足不同条件的近空间风场相关系数或相关长度。将其代入上述公式可以得到一个既具有统计协方差模型的各向异性特性，又具备高斯模型便利性和快速性的协方差模型。
本研究假设近空间风场在水平和垂直方向上是可分割的。近空间风场的背景场误差协方差由水平背景场误差协方差和垂直背景场误差协方差分别计算。水平背景场误差协方差可以进一步划分为子午背景场误差协方差和纬向背景场误差协方差。

3.3.3 水平背景场误差协方差
在本研究中，使用HWM14模拟的风场数据作为背景场数据的来源。由于近空间风场的周期性特征，通常表现为季节性[36]，因此本研究将模型在80–100公里高度获得的季节平均值视为真实值。在此基础上，对背景场的误差进行统计分析。背景场数据的划分条件见表2。

Ledoit–Wolf收缩协方差估计是一种基于协方差矩阵和单位矩阵的广义线性组合。当样本数量小于变量维度时，它的性能通常更好，并且可以确保协方差矩阵是正定的[37]。在本研究中，使用Ledoit–Wolf收缩协方差估计器来估计80–100公里高度的背景场误差协方差。得到的值设为$M$，结合以下公式可以获得风场相关系数：
$$
\lambda = M \sqrt{\frac{t \tau a c e}{M T}
$$
在计算水平相关系数时，忽略了垂直方向变化对水平相关系数的影响。以下图表显示了不同季节纬向风和子午风的水平相关长度分布。相关长度是指参考点与相关区域边界之间的距离（本研究将相关系数的阈值设为0.75）。图2、图3、图4和图5展示了基于HWM14模拟的背景场数据计算的曲靖地区春季、夏季、秋季和冬季的平均风场相关长度。图表显示，在不同季节，近空间风场的纬向相关长度不会随季节变化，纬向风和子午风的纬向相关长度是一致的。同时，子午风的子午相关长度也不随季节变化，而在夏季和秋季变小，其空间分布与春季和冬季有显著差异。

3.3.4 垂直背景场误差协方差
类似地，在计算垂直相关系数时，忽略了垂直方向变化对垂直相关系数的影响。根据曲靖气象雷达数据集的有效数据范围和垂直分辨率，本研究将HWM模拟的数据范围设置在垂直方向上为80公里到100公里，高度分辨率为2公里。根据垂直相关长度的定义，近空间纬向风和子午风的垂直相关长度分别为2公里或4公里，可以从垂直相关长度值得到垂直相关系数。
图6和图7展示了不同季节纬向风和子午风的垂直相关系数分布。图表显示，在不同季节，近空间风场的垂直相关系数基本不随季节变化，纬向风和子午风的垂直相关系数是一致的。

3.3.5 观测算子
观测算子H表示观测误差协方差$R$与背景模型误差协方差$B$的相对大小和空间分布，决定了观测场和背景场对同化结果的相对影响权重和空间结构。在本研究中，观测算子是一个由0和1组成的矩阵，其中矩阵的行数和列数分别与观测场和背景场的数据维度相对应。0和1的值取决于该位置的观测数据是否缺失；如果观测值缺失，则该位置的值为0。
H的功能是实现背景场与观测场之间的匹配。通过误差协方差的空间传播来实现背景场的校正。

3.3.6 时间相关性
一般认为，近空间风场的状态变化与时间有很强的相关性。在本研究中，使用高斯–马尔可夫算法来描述这种相关性。基于马尔可夫理论[38]，在背景场模型的计算结果中加入了一个时间衰减误差项作为卡尔曼滤波器的预测方程。误差项是当前时间步的状态分析值与模型计算值之间的差异。衰减系数是与时标相关的指数函数，因此上述公式可以表示为：
$$
X_{b}^{t+1} = (X_{a}^t - X_{b}^t) \exp(-\Delta_t \tau) + X_{b}^t
$$
$$
X_{b}^{t+1} = X_{a}^t \exp(-2\Delta_t \tau) - X_{b}^t \exp(-2\Delta_t \tau) + X_{b}^t
$$
上述公式中，$\tau$是时间相关尺度，随宇宙时间变化。下图显示了$\tau$与局部时间的分布。图8表明，HWM14风场的输出结果受到半日潮汐效应的支配，时间相关尺度随宇宙时间变化呈现出双峰结构。

4. 同化结果
根据上述同化过程，本研究基于2022年曲靖地区的气象雷达观测数据和HWM14进行了同化实验。图9（下图）展示了未经同化的独立气象雷达风场剖面与相应背景场和同化后的分析场的对比分析。实线红色代表独立的气象雷达观测值（真实值）。黑色虚线表示HWM14的背景风场。卡尔曼滤波器同化后的风场分析剖面为紫色实线。为了更全面地验证所提出方法的性能，我们还加入了全球大气社区气候模型（Waccm）和三维变分（3D-Var）同化的结果。图表显示，在80–100公里高度范围内，背景风场值与气象雷达风场的真实值之间存在显著差异。与背景风场值相比，卡尔曼滤波器同化后的风场分析值更接近气象雷达的真实值，同化效果明显。此外，如图9所示，卡尔曼滤波器同化得到的水平风场（紫色线）在垂直结构和幅度上都与3D-Var同化结果（蓝色线）非常吻合。这种一致性进一步验证了卡尔曼滤波器同化技术在可靠分析和校正近空间区域水平风场方面的能力。

3.3.4 垂直背景场误差协方差
在计算垂直相关系数时，同样忽略了垂直方向变化对垂直相关系数的影响。根据曲靖气象雷达数据集的有效数据范围和垂直分辨率，本研究将HWM模拟的数据范围设置在垂直方向上为80公里到100公里，高度分辨率为2公里。根据垂直相关长度的定义，近空间纬向风和子午风的垂直相关长度分别为2公里或4公里，可以从垂直相关长度值得到垂直相关系数。
图6和图7展示了不同季节纬向风和子午风的垂直相关系数分布。图表显示，在不同季节，近空间风场的垂直相关系数基本不随季节变化，纬向风和子午风的垂直相关系数是一致的。

3.3.5 观测算子
观测算子H表示观测误差协方差$R$与背景模型误差协方差$B$的相对大小和空间分布，决定了观测场和背景场对同化结果的相对影响权重和空间结构。在本研究中，观测算子是一个由0和1组成的矩阵，其中矩阵的行数和列数分别与观测场和背景场的数据维度相关。0和1的值取决于该位置的观测数据是否缺失；如果观测值缺失，则该位置的值为0。
H的功能是在背景场和观测场之间实现匹配。通过误差协方差的空间传播来实现背景场的校正。

3.3.6 时间相关性
普遍认为，近空间风场的状态变化与时间有很强的相关性。在本研究中，使用高斯–马尔可夫算法来描述这种相关性。基于马尔可夫理论[38]，在背景场模型的计算结果中加入了一个时间衰减误差项作为卡尔曼滤波器的预测方程。误差项是当前时间步的状态分析值与模型计算值之间的差异。衰减系数是与时间尺度相关的指数函数，因此上述公式可以表示为：
$$
X_{b}^{t+1} = (X_{a}^t - X_{b}^t) \exp(-\Delta_t \tau) + X_{b}^t
$$
$$
X_{b}^{t+1} = X_{a}^t \exp(-2\Delta_t \tau) - X_{b}^t \exp(-2\Delta_t \tau) + X_{b}^t
$$
上述公式中，$\tau$是时间相关尺度，随宇宙时间变化。下图显示了$\tau$与局部时间的分布。图8表明，HWM14风场的输出结果受到半日潮汐效应的支配，时间相关尺度随宇宙时间变化而呈现出双峰结构。

4. 同化结果
根据上述同化过程，本研究基于2022年曲靖地区的气象雷达观测数据和HWM14进行了同化实验。图9（下图）展示了未经同化的独立气象雷达风场剖面与相应背景场和同化后分析场的对比分析。实线红色代表独立的气象雷达观测值（真实值）。黑色短虚线表示HWM14的背景风场。卡尔曼滤波器同化后的风场分析剖面为紫色实线。为了更全面地验证所提出方法的效果，我们还加入了全球大气社区气候模型（Waccm）和三维变分（3D-Var）同化的结果。图表显示，在80–100公里高度范围内，背景风场值与气象雷达风场的真实值之间存在显著差异。与背景风场值相比，卡尔曼滤波器同化后的风场分析值更接近气象雷达的真实值，同化效果明显。此外，如图9所示，卡尔曼滤波器同化得到的水平风场（紫色线）在垂直结构和幅度上都与3D-Var同化结果（蓝色线）非常吻合。这种一致性进一步验证了卡尔曼滤波器同化技术在可靠分析和校正近空间区域水平风场方面的能力。比较结果显示，同化风场分析值与气象雷达风场值之间的差异显著减小。在90公里高度处，纬向风速的平均最大误差为12.0米/秒，比之前增加了55.0%；而在100公里高度处，经向风速的平均最大误差为14.2米/秒，比之前增加了53.4%。此外，同化风场分析值与气象雷达风场之间的偏差标准差也显著降低。图10展示了卡尔曼滤波器同化前的误差统计分析（纬向和经向）。图11展示了卡尔曼滤波器同化后的误差统计分析（纬向和经向）。

5. 讨论
本研究通过计算不同季节的相关长度，实现了背景误差协方差矩阵的季节性自适应设置。这种方法使同化系统能够根据大气动力结构的季节性变化，更适当地调节气象雷达观测信息的空间传播和权重分配。与使用固定相关长度的方案相比，这一改进显著提高了跨季节同化分析的准确性和稳定性。此外，虽然图2和图5中经向和纬向成分的模式非常相似，但在图3和图4中它们在水平方向上显示出显著差异。这些结果直接反映了低纬度和中纬度高层大气环流的季节性变化。在春季和冬季，环流结构相对稳定，以西风为主导，大尺度气流占主导地位[39,40]，这导致风场扰动在水平方向上具有相对均匀的统计特性，从而使经向和纬向成分的图表非常相似。相比之下，在夏季和秋季，环流发生显著调整，表现为西风急流的北移、东风系统的形成、对流活动的增强以及方向性明确的大气波动[41]，风场扰动的空间协方差结构表现出强烈的各向异性。这导致了经向和纬向成分在水平相关长度上的系统差异。这种对比生动地表明，背景误差协方差的结构并非固定不变，而是与主导环流和波动活动的季节性演变紧密相关。

6. 结论
基于卡尔曼滤波器同化技术，并结合HWM14模型，本研究评估了利用曲靖地区的气象雷达观测数据，在80-100公里高度范围内对二维大气风场进行数据同化的技术。同化结果显示，基于卡尔曼滤波器的同化技术对80-100公里高度范围内的大气风场具有良好的重建能力。与HWM14模型相比，同化后的大气风场数据的整体准确性提高了30%以上。因此，所提出的卡尔曼滤波器同化技术能够整合HWM14数据和雷达风场观测数据，非常适合用于建立大气风场数值预报的初始值。