一種星載GNSS-R 海風反演的卡爾曼濾波模型

李中奎，張 波，楊東凱，張國棟

（北京航空航天大學 電子與信息工程學院，北京 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 公開數據。