李中奎，张 波，杨东凯，张国栋

（北京航空航天大学 电子与信息工程学院，北京 100191）

0 引言

星载全球卫星导航系统反射计（global navigation satellite system-reflectometry, GNSS-R）技术是近年来发展起来的1 种新兴的海洋遥感技术。GNSS 反射信号的研究始于 20 世纪 90 年代前后，目前已在多个遥感领域得到验证，包括海风海浪[1-2]、土壤湿度[3]、雪深[4]、海冰[5]、海面溢油[6]等多个地表参数的反演。在星载GNSS-R 领域，海面风场的反演一直是研究热点。2000 年，文献[7]采用Kirchhoff 近似几何光学建立了 GNSS 海面散射信号时延多普勒2 维相关功率模型。同年文献[8]采用海浪谱模型，研究了全球定位系统（global positioning system, GPS）海面反射信号极化特性的变化情况，分析了海面测风的可行性。2002 年，文献[9]首次在星载平台上探测到了 GNSS 反射信号，随后英国国家空间中心在2003 年发射了首个搭载 GPS-R 接收机的英国灾难监测星座（United Kingdom-disaster monitor constellation, UK-DMC）卫星[10]。2013 年，文献[11]将 4 个延迟多普勒图（delay Doppler map, DDM）观测量（DDM 质心、DDM 几何中心、DDM Taxicab 质心以及 DDM 加权面积）与风速建立了线性关系，结果表明2 种DDM质心能够提供更好的反演精度。2014 年，英国Surrey卫星公司成功研制了技术演示卫星（technology demonstration satellite-1, TDS-1）的 GNSS-R 模块并采集到了大量DDM 数据[12]。2016 年，文献[13]利用 TDS-1 数据，通过 DDM 峰值反解双基散射系数，并与高级散射仪（advanced scatterometer,ASCAT）数据建立经验模型进行风速反演，得到反演风速的均方根误差为2.20 m/s。同年，美国航空航天局（ National Aeronautics and Space Administration, NASA）成功发射了旋风全球卫星导航系统（cyclone global navigation satellite system，CYGNSS）。针对CYGNSS 反演风速时，中高风速探测精度不足和空间分辨率低等问题，文献[14]提出了基于广义线性的 GNSS-R 观测量，结果表明，通过主成分分析定义的广义线性观测量，能够获得最佳反演性能。同年，文献[15]提出了基于最小方差的风速估计器，用于CYGNSS 风速反演。

本文针对现有的星载 GNSS-R 测风数据利用率低、空间覆盖率差等问题，提出了 1 种卡尔曼（Kalman）滤波模型用于星载GNSS-R 风速反演，并利用美国NASA 公开的CYGNSS 实测数据与欧洲中期天气预报中心（ European Centre for Medium-range Weather Forecasts，ECMWF）的再分析风速数据集，对该方法进行了验证。

1 星载GNSS-R 海面风速反演原理

GNSS 海面散射信号时延多普勒2 维相关功率模型[16]为

TDS-1 的DDM 数据是无量纲的计数值，且缺少相关材料将其转化为以瓦为单位的散射信号绝对功率值，因此文献[2]采用定义校正因子，对相关参数进行了校正。不同于 TDS-1，CYGNSS 的DDM 数据集提供了更多的校准参数[17]，因此，本文没有采用定义校正因子对0σ进行近似替代的方法。本文求解0σ采用的校准流程如图1 所示，首先由 DDM 原始数据得到校准后的双基雷达散射截面（bistatic radar cross section, BRCS），进而再结合有效散射面积查找表得到0σ。

图1 求解归一化双基雷达散射截面的校准流程

2 基于Kalman 滤波模型的风速反演算法

Kalman 滤波适用于线性、离散和有限维系统。每 1 个有外部变量的自回归移动平均系统或可用有理传递函数表示的系统，都可以转换成用状态空间表示的系统，从而能用 Kalman 滤波进行计算[18]。Kalman 滤波的关键是状态方程和观测方程的建立，本文运用差分自回归移动平均模型（autoregressive integrated moving average model, ARIMA）建立Kalman 滤波的状态方程和观测方程并得到风速的预测值，同时利用校准后的NBRCS 与ECMWF的海面风速，建立地球物理模型函数 （geophysical model function, GMF）经验模型，再将GMF 得到的风速值作为观测值。

2.1 GMF 经验模型测风算法

文献[15]首次提出了综合利用 DDM 均值、DDM 方差、前沿斜率等5 个DDM 衍生出的观测量进行风速反演。由于星载平台高度及卫星下行链路的传输速率限制，星载 GNSS-R 系统一般采用延迟多普勒平均图（delay-Doppler map average,DDMA）和前沿斜率反演海面风速；其次，受制于星载探测的空间分辨率的要求，星载平台风速反演中只能使用更小尺寸的DDM，这意味着相对于前沿斜率，DDMA 对风速变化具有更高的敏感性。因此，本文将围绕DDMA 对海面风速进行反演。

DDMA 定义为围绕镜面反射点的特定时延-多普勒窗内NBRCS 的均值。DDMA 是利用延迟多普勒图（delay-Doppler map, DDM）中对风速最敏感的1 块区域，并通过求均值降低了噪声的影响，因此，常采用DDMA 与风速建立经验函数关系。本文以ECMWF 的再分析数据集作为同比风速数据，建立了DDMA 与海面风速之间的GMF 经验模型。用于拟合的函数形式为

式中：A、B、C、D是待定的拟合参数；u10为海面风速；为DDMA；e为自然常数。

2.2 ARIMA 模型

ARIMA 模型是1 种时间序列建模方法[19]，它是由自回归滑动平均模型扩展而来。ARIMA 模型具有良好的短期预测能力，根据风速时间序列的自相关性，通过ARIMA 模型的短期预测，可得到下一时刻的风速预测值，并将其作为Kalman 滤波的先验值。同时，由 ARIMA 模型可方便得到Kalman 滤波所需要的状态方程。

ARIMA(p,d,q)模型总共有3 项，即自回归项（auto regressive, AR）、差分项（integrated, I）、移动平均项（moving average, MA），对应参数为p、d、q。通过对不同的（p、d、q）组合测试，可以找到最合适的模型参数。

建立ARIMA 模型的步骤如下：

步骤①：平稳化分析。对数据进行预处理使该时间序列满足建模的前提条件，同时能提高预测精度。

步骤②：模型定阶。在得到稳定的时间序列后，需要确定合适的p和q值来确定最终的预测模型。在自回归模型中，阶数p可由赤池信息量准则（Akaike information criterion, AIC）进行定阶，合理的阶数p，需要满足精度要求的同时，又能避免运算太过复杂。

步骤③：求解自回归系数和滑动平均系数。求解自回归系数及滑动平均系数，可用 Yule-Walker方程法或最小二乘估计法进行求解。

步骤④：风速预测。通过ARIMA 模型得到风速预测的表达式，对下一时刻的风速值进行预测。

为验证 ARIMA 模型的预测效果，本文将ECMWF 再分析风速数据通过双线性插值得到1 段连续的沿实际镜面反射点轨迹分布的风速时间序列，并用步骤①～步骤④对该时间序列上的风速进行逐步预测。图 2 显示了选取的部分时间段的ECMWF 风速和预测风速。

通过对选取的样本点进行计算，ARIMA 模型预测风速与ECMWF 风速的均方根误差为0.20 m/s。同时由图2 可以看出，ARIMA 模型在风速值波动大的时刻附近的预测效果较差，风速值较稳定的时刻附近，预测效果较好。该结果可以证明，利用ARIMA模型对遥感风速时间序列进行预测的可行性。

2.3 Kalman 滤波模型

本文由ARIMA 模型建立1 个线性的时间序列模型后，导出Kalman 滤波的状态方程和观测方程，然后将DDM 观测量做为GMF 经验模型的输入，求解得到初步的反演风速值即为实际观测值，从而建立起Kalman 滤波模型。

首先推导状态方程和观测方程，ARIMA 的预测方程可表示为

式中：v(t)为当前时刻的预测风速；φ1,φ2, …,为自回归系数，下际p为AR 模型的阶数；Δt为时间间隔；ε(t)为随机数，它的均值为 0。令则式（3）可表示为

式中εk为k时刻随机数。令则状态方程和观测方程可表示为：

用矩阵形式表达为：

式中：Xk、w k、yk、v k分别为k时刻的状态向量、噪声矩阵、观测值和观测误差；M为状态转移矩阵；H为观测算子矩阵。X、M、H的表达式分别为

Kalman 滤波方法对海面风速进行反演优化的迭代过程主要分为5 步：

1）对风速状态进行一步预测，其计算方法为

2）计算风速状态先验估计值的均方误差，其计算方法为

3）计算Kalman 增益矩阵，其计算方法为

4）风速优化更新，其计算方法为

5）对状态估计值的均方误差进行更新，其计算方法为

上述 Kalman 滤波在实现过程中的关键点在于：实际测量值yk、过程噪声的协方差Q及测量噪声的协方差R的取值问题。文献[18-19]中Q和R的数据取经验值，缺少取值的指导依据。本文实际测量值采用的是 GMF 经验模型得出的风速值，R值则是通过单独对GMF 经验模型反演方法的结果进行统计得到，一般取其均方根误差。而对于测量噪声协方差Q，这里则是沿用经验值来代替。文献[20]对此类不准确方差下，带随机系数矩阵的 Kalman 滤波器的稳定性进行了分析，在满足系数矩阵的有界性、条件能观测性及噪声和初始误差的有界性条件下，Kalman 滤波器可保持稳定。本文建立的Kalman 滤波器满足上述3 个条件。

3 实验与结果分析

3.1 数据集描述

2016 年12 月美国NASA 成功发射了由8 颗轨道高度为 510 km、轨道倾角为 35°的微小卫星组成的CYGNSS 星座。每个CYGNSS 气象台搭载有1 台延迟多普勒测绘仪器(delay Doppler mapping instrument,DDMI)，每个DDMI 拥有4 个反射信号接收通道，并同时以1 Hz 的频率对GPS L1 C/A 码信号同时进行采集，从而生成DDM 数据，因此CYGNSS 在1 s 时间内，可以在全球范围内进行32 次风速测量[17]。

本文采用2018-09-23—24 期间的CYGNSS 数据作为观测数据，并以 ECMWF 的再分析风速数据集作为同比风速。ECMWF 提供时空分辨率为1 h、0.5°的海面风速，是风速反演研究中的良好同比数据源。数据的匹配方法是对 ECMWF 风速数据集进行空间双线性插值和时间线性插值，得到对应时间和空间内的风速值。

3.2 风速反演与结果分析

3.2.1 GMF 经验模型反演结果

CYGNSS 遥感数据质量不及传统遥感卫星的级别，为建立更准确的GMF 经验模型，以便进行风速反演，需要对数据质量进行控制。为剔除部分信噪比较低的数据，引入距离校正增益（rangecorrected gain, RCG）阈值对数据进行筛选。RCG定义为镜面反射点处的接收机天线增益乘以对应的距离损耗，可表示为

图3 数据样本百分比随RCG 阈值变化曲线

将经过数据质量控制筛选后的2018-09-23 的DDMA 数据作为训练集，2018-09-24 的 DDMA 数据作为测试集。由训练集通过式（2）拟合得到DDMA 与 ECMWF 风速之间的GMF 经验模型，并使用测试集分析所建立的模型效果。拟合结果如图4 所示。将测试集按照该GMF 经验模型进行风速反演并与 ECMWF 风速进行对比，对比结果散点图如图5 所示，均方根误差为2.18 m/s。

图4 DDMA 与风速的关系散点图及拟合曲线

图5 GMF 经验模型风速反演值与ECMWF 同比风速对比散点图

3.2.2 Kalman 滤波模型反演结果

经验模型函数反演方法反演精度较低，同时在执行严格的质量控制后会损失大量数据样本。本文采用上述Kalman 滤波模型，对同样的测试集数据进行了处理，整体反演流程如图6 所示。

图6 基于Kalman 滤波模型的海面风速反演流程

在实际数据处理过程中，CYGNSS 的 DDMI的每个通道以1 Hz 的频率采集数据，因此24 h 内可得到86 400 个数据样本，但并非所有数据样本都可以用于风速反演。比如镜面反射点经过陆地时采集的数据，数据质量太差（RCG＜3）导致的不可用观测量，以及 DDMI 每隔 60 s 进行 4 s 的黑体负载校正期间导致的空白数据，此时通过 GMF经验模型无法得到风速的观测值。本文利用ARIMA 模型及 Kalman 滤波模型的预测机制可对这2 种无效数据进行补充，从而提升数据利用率。

基于上述方法，对2018-09-24 的CYGNSS 的01号星的数据进行了反演，为了证明模型的可用性，图 7（a）～图 7（d）展示了基于 Kalman 滤波的 1 个阶段的反演结果，图 7（a）～图 7（c）分别显示了Kalman 滤波的前150 s、中间的150 s，及最后的150 s的结果。而图7（d）显示的是1 个阶段内的总体效果。

图7 Kalman 滤波模型反演得到的局部时间序列

由图7（a）～图7（c）可得：相对于ECMWF同比风速的时间序列，GMF 经验模型反演得到的风速波动较大，同时存在空白数据及异常值；而通过 Kalman 滤波模型得到的反演风速值相对于GMF 更加收敛，减少了异常值，同时对短时间内的空白数据进行了有效补充。由于连续的长时间空白数据会导致 Kalman 滤波模型预测能力变差，因此本文对空白数据样本点的补充限制在连续 5 个样本点以内。

对2018-09-24 的24 h 的数据进行反演的结果如表1 所示。

表 1 2018-09-24 CYGNSS 01 号星 1 通道 24 h 数据反演结果

由表1 可知：利用Kalman 滤波模型进行反演的均方根误差为2.08 m/s，较经验模型法有所提升，同时反演得到的风速数据样本数量从经验模型的56 890 提升到了 63 175，提升了 11.04 %。与ECMWF 风速对比的散点图如图8 所示。

图8 Kalman 滤波模型反演风速值与ECMWF 同比风速对比散点图

通过图5 与图8 的对比可以发现，反演得到的整体数据质量得到了提升，散点图更加收敛，误差较大的异常值有所减少。为了进一步分析模型的风速反演性能，本文将RCG 以间隔值为3 进行分段，统计各个分段的风速反演精度，结果如图 9 所示。

图9 各RCG 段风速反演误差

由图9 可以看出，当RCG 较小时的风速反演效果较差，同时 Kalman 滤波反演模型在 RCG 较小时具有比经验模型更好的反演效果。这是因为RCG 较小时接收机天线增益较低，得到的数据信噪比较低，从而导致风速反演结果较差；而Kalman滤波模型对该情形导致的异常值进行了校正，从而得到了更好的结果。

4 结束语

本文对星载GNSS-R 海面风速反演方法进行了研究。首先推导了用于风速反演的归一化散射雷达截面及 DDMA 的计算流程；然后利用 ARIMA模型及 GMF 经验模型建立 Kalman 滤波模型，并进行风速反演。本文利用CYGNSS 星载数据对该Kalman 滤波模型进行了验证，反演得到的均方根误差为 2.08 m/s，较 GMF 经验模型反演方法的2.18 m/s 有所提升，与此同时数据利用率提高了11.04 %，进而提升了数据的空间覆盖率。

致谢：感谢美国航空航天局（NASA）提供的CYGNSSS 公开数据。