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

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

（北京航空航天大學 電子與信息工程學院，北京 100191）

0 引言

星載全球衛(wèi)星導航系統(tǒng)反射計（global navigation satellite system-reflectometry, GNSS-R）技術(shù)是近年來發(fā)展起來的1 種新興的海洋遙感技術(shù)。GNSS 反射信號的研究始于 20 世紀 90 年代前后，目前已在多個遙感領(lǐng)域得到驗證，包括海風海浪[1-2]、土壤濕度[3]、雪深[4]、海冰[5]、海面溢油[6]等多個地表參數(shù)的反演。在星載GNSS-R 領(lǐng)域，海面風場的反演一直是研究熱點。2000 年，文獻[7]采用Kirchhoff 近似幾何光學建立了 GNSS 海面散射信號時延多普勒2 維相關(guān)功率模型。同年文獻[8]采用海浪譜模型，研究了全球定位系統(tǒng)（global positioning system, GPS）海面反射信號極化特性的變化情況，分析了海面測風的可行性。2002 年，文獻[9]首次在星載平臺上探測到了 GNSS 反射信號，隨后英國國家空間中心在2003 年發(fā)射了首個搭載 GPS-R 接收機的英國災難監(jiān)測星座（United Kingdom-disaster monitor constellation, UK-DMC）衛(wèi)星[10]。2013 年，文獻[11]將 4 個延遲多普勒圖（delay Doppler map, DDM）觀測量（DDM 質(zhì)心、DDM 幾何中心、DDM Taxicab 質(zhì)心以及 DDM 加權(quán)面積）與風速建立了線性關(guān)系，結(jié)果表明2 種DDM質(zhì)心能夠提供更好的反演精度。2014 年，英國Surrey衛(wèi)星公司成功研制了技術(shù)演示衛(wèi)星（technology demonstration satellite-1, TDS-1）的 GNSS-R 模塊并采集到了大量DDM 數(shù)據(jù)[12]。2016 年，文獻[13]利用 TDS-1 數(shù)據(jù)，通過 DDM 峰值反解雙基散射系數(shù)，并與高級散射儀（advanced scatterometer,ASCAT）數(shù)據(jù)建立經(jīng)驗模型進行風速反演，得到反演風速的均方根誤差為2.20 m/s。同年，美國航空航天局（ National Aeronautics and Space Administration, NASA）成功發(fā)射了旋風全球衛(wèi)星導航系統(tǒng)（cyclone global navigation satellite system，CYGNSS）。針對CYGNSS 反演風速時，中高風速探測精度不足和空間分辨率低等問題，文獻[14]提出了基于廣義線性的 GNSS-R 觀測量，結(jié)果表明，通過主成分分析定義的廣義線性觀測量，能夠獲得最佳反演性能。同年，文獻[15]提出了基于最小方差的風速估計器，用于CYGNSS 風速反演。

本文針對現(xiàn)有的星載 GNSS-R 測風數(shù)據(jù)利用率低、空間覆蓋率差等問題，提出了 1 種卡爾曼（Kalman）濾波模型用于星載GNSS-R 風速反演，并利用美國NASA 公開的CYGNSS 實測數(shù)據(jù)與歐洲中期天氣預報中心（ European Centre for Medium-range Weather Forecasts，ECMWF）的再分析風速數(shù)據(jù)集，對該方法進行了驗證。

1 星載GNSS-R 海面風速反演原理

GNSS 海面散射信號時延多普勒2 維相關(guān)功率模型[16]為

TDS-1 的DDM 數(shù)據(jù)是無量綱的計數(shù)值，且缺少相關(guān)材料將其轉(zhuǎn)化為以瓦為單位的散射信號絕對功率值，因此文獻[2]采用定義校正因子，對相關(guān)參數(shù)進行了校正。不同于 TDS-1，CYGNSS 的DDM 數(shù)據(jù)集提供了更多的校準參數(shù)[17]，因此，本文沒有采用定義校正因子對0σ進行近似替代的方法。本文求解0σ采用的校準流程如圖1 所示，首先由 DDM 原始數(shù)據(jù)得到校準后的雙基雷達散射截面（bistatic radar cross section, BRCS），進而再結(jié)合有效散射面積查找表得到0σ。

圖1 求解歸一化雙基雷達散射截面的校準流程

2 基于Kalman 濾波模型的風速反演算法

Kalman 濾波適用于線性、離散和有限維系統(tǒng)。每 1 個有外部變量的自回歸移動平均系統(tǒng)或可用有理傳遞函數(shù)表示的系統(tǒng)，都可以轉(zhuǎn)換成用狀態(tài)空間表示的系統(tǒng)，從而能用 Kalman 濾波進行計算[18]。Kalman 濾波的關(guān)鍵是狀態(tài)方程和觀測方程的建立，本文運用差分自回歸移動平均模型（autoregressive integrated moving average model, ARIMA）建立Kalman 濾波的狀態(tài)方程和觀測方程并得到風速的預測值，同時利用校準后的NBRCS 與ECMWF的海面風速，建立地球物理模型函數(shù) （geophysical model function, GMF）經(jīng)驗模型，再將GMF 得到的風速值作為觀測值。

2.1 GMF 經(jīng)驗模型測風算法

文獻[15]首次提出了綜合利用 DDM 均值、DDM 方差、前沿斜率等5 個DDM 衍生出的觀測量進行風速反演。由于星載平臺高度及衛(wèi)星下行鏈路的傳輸速率限制，星載 GNSS-R 系統(tǒng)一般采用延遲多普勒平均圖（delay-Doppler map average,DDMA）和前沿斜率反演海面風速；其次，受制于星載探測的空間分辨率的要求，星載平臺風速反演中只能使用更小尺寸的DDM，這意味著相對于前沿斜率，DDMA 對風速變化具有更高的敏感性。因此，本文將圍繞DDMA 對海面風速進行反演。

DDMA 定義為圍繞鏡面反射點的特定時延-多普勒窗內(nèi)NBRCS 的均值。DDMA 是利用延遲多普勒圖（delay-Doppler map, DDM）中對風速最敏感的1 塊區(qū)域，并通過求均值降低了噪聲的影響，因此，常采用DDMA 與風速建立經(jīng)驗函數(shù)關(guān)系。本文以ECMWF 的再分析數(shù)據(jù)集作為同比風速數(shù)據(jù)，建立了DDMA 與海面風速之間的GMF 經(jīng)驗模型。用于擬合的函數(shù)形式為

式中：A、B、C、D是待定的擬合參數(shù)；u10為海面風速；為DDMA；e為自然常數(shù)。

2.2 ARIMA 模型

ARIMA 模型是1 種時間序列建模方法[19]，它是由自回歸滑動平均模型擴展而來。ARIMA 模型具有良好的短期預測能力，根據(jù)風速時間序列的自相關(guān)性，通過ARIMA 模型的短期預測，可得到下一時刻的風速預測值，并將其作為Kalman 濾波的先驗值。同時，由 ARIMA 模型可方便得到Kalman 濾波所需要的狀態(tài)方程。

ARIMA(p,d,q)模型總共有3 項，即自回歸項（auto regressive, AR）、差分項（integrated, I）、移動平均項（moving average, MA），對應參數(shù)為p、d、q。通過對不同的（p、d、q）組合測試，可以找到最合適的模型參數(shù)。

建立ARIMA 模型的步驟如下：

步驟①：平穩(wěn)化分析。對數(shù)據(jù)進行預處理使該時間序列滿足建模的前提條件，同時能提高預測精度。

步驟②：模型定階。在得到穩(wěn)定的時間序列后，需要確定合適的p和q值來確定最終的預測模型。在自回歸模型中，階數(shù)p可由赤池信息量準則（Akaike information criterion, AIC）進行定階，合理的階數(shù)p，需要滿足精度要求的同時，又能避免運算太過復雜。

步驟③：求解自回歸系數(shù)和滑動平均系數(shù)。求解自回歸系數(shù)及滑動平均系數(shù)，可用 Yule-Walker方程法或最小二乘估計法進行求解。

步驟④：風速預測。通過ARIMA 模型得到風速預測的表達式，對下一時刻的風速值進行預測。

為驗證 ARIMA 模型的預測效果，本文將ECMWF 再分析風速數(shù)據(jù)通過雙線性插值得到1 段連續(xù)的沿實際鏡面反射點軌跡分布的風速時間序列，并用步驟①～步驟④對該時間序列上的風速進行逐步預測。圖 2 顯示了選取的部分時間段的ECMWF 風速和預測風速。

通過對選取的樣本點進行計算，ARIMA 模型預測風速與ECMWF 風速的均方根誤差為0.20 m/s。同時由圖2 可以看出，ARIMA 模型在風速值波動大的時刻附近的預測效果較差，風速值較穩(wěn)定的時刻附近，預測效果較好。該結(jié)果可以證明，利用ARIMA模型對遙感風速時間序列進行預測的可行性。

2.3 Kalman 濾波模型

本文由ARIMA 模型建立1 個線性的時間序列模型后，導出Kalman 濾波的狀態(tài)方程和觀測方程，然后將DDM 觀測量做為GMF 經(jīng)驗模型的輸入，求解得到初步的反演風速值即為實際觀測值，從而建立起Kalman 濾波模型。

首先推導狀態(tài)方程和觀測方程，ARIMA 的預測方程可表示為

式中：v(t)為當前時刻的預測風速；φ1,φ2, …,為自回歸系數(shù)，下際p為AR 模型的階數(shù)；Δt為時間間隔；ε(t)為隨機數(shù)，它的均值為 0。令則式（3）可表示為

式中εk為k時刻隨機數(shù)。令則狀態(tài)方程和觀測方程可表示為：

用矩陣形式表達為：

式中：Xk、w k、yk、v k分別為k時刻的狀態(tài)向量、噪聲矩陣、觀測值和觀測誤差；M為狀態(tài)轉(zhuǎn)移矩陣；H為觀測算子矩陣。X、M、H的表達式分別為

Kalman 濾波方法對海面風速進行反演優(yōu)化的迭代過程主要分為5 步：

1）對風速狀態(tài)進行一步預測，其計算方法為

2）計算風速狀態(tài)先驗估計值的均方誤差，其計算方法為

3）計算Kalman 增益矩陣，其計算方法為

4）風速優(yōu)化更新，其計算方法為

5）對狀態(tài)估計值的均方誤差進行更新，其計算方法為

上述 Kalman 濾波在實現(xiàn)過程中的關(guān)鍵點在于：實際測量值yk、過程噪聲的協(xié)方差Q及測量噪聲的協(xié)方差R的取值問題。文獻[18-19]中Q和R的數(shù)據(jù)取經(jīng)驗值，缺少取值的指導依據(jù)。本文實際測量值采用的是 GMF 經(jīng)驗模型得出的風速值，R值則是通過單獨對GMF 經(jīng)驗模型反演方法的結(jié)果進行統(tǒng)計得到，一般取其均方根誤差。而對于測量噪聲協(xié)方差Q，這里則是沿用經(jīng)驗值來代替。文獻[20]對此類不準確方差下，帶隨機系數(shù)矩陣的 Kalman 濾波器的穩(wěn)定性進行了分析，在滿足系數(shù)矩陣的有界性、條件能觀測性及噪聲和初始誤差的有界性條件下，Kalman 濾波器可保持穩(wěn)定。本文建立的Kalman 濾波器滿足上述3 個條件。

3 實驗與結(jié)果分析

3.1 數(shù)據(jù)集描述

2016 年12 月美國NASA 成功發(fā)射了由8 顆軌道高度為 510 km、軌道傾角為 35°的微小衛(wèi)星組成的CYGNSS 星座。每個CYGNSS 氣象臺搭載有1 臺延遲多普勒測繪儀器(delay Doppler mapping instrument,DDMI)，每個DDMI 擁有4 個反射信號接收通道，并同時以1 Hz 的頻率對GPS L1 C/A 碼信號同時進行采集，從而生成DDM 數(shù)據(jù)，因此CYGNSS 在1 s 時間內(nèi)，可以在全球范圍內(nèi)進行32 次風速測量[17]。

本文采用2018-09-23—24 期間的CYGNSS 數(shù)據(jù)作為觀測數(shù)據(jù)，并以 ECMWF 的再分析風速數(shù)據(jù)集作為同比風速。ECMWF 提供時空分辨率為1 h、0.5°的海面風速，是風速反演研究中的良好同比數(shù)據(jù)源。數(shù)據(jù)的匹配方法是對 ECMWF 風速數(shù)據(jù)集進行空間雙線性插值和時間線性插值，得到對應時間和空間內(nèi)的風速值。

3.2 風速反演與結(jié)果分析

3.2.1 GMF 經(jīng)驗模型反演結(jié)果

CYGNSS 遙感數(shù)據(jù)質(zhì)量不及傳統(tǒng)遙感衛(wèi)星的級別，為建立更準確的GMF 經(jīng)驗模型，以便進行風速反演，需要對數(shù)據(jù)質(zhì)量進行控制。為剔除部分信噪比較低的數(shù)據(jù)，引入距離校正增益（rangecorrected gain, RCG）閾值對數(shù)據(jù)進行篩選。RCG定義為鏡面反射點處的接收機天線增益乘以對應的距離損耗，可表示為

圖3 數(shù)據(jù)樣本百分比隨RCG 閾值變化曲線

將經(jīng)過數(shù)據(jù)質(zhì)量控制篩選后的2018-09-23 的DDMA 數(shù)據(jù)作為訓練集，2018-09-24 的 DDMA 數(shù)據(jù)作為測試集。由訓練集通過式（2）擬合得到DDMA 與 ECMWF 風速之間的GMF 經(jīng)驗模型，并使用測試集分析所建立的模型效果。擬合結(jié)果如圖4 所示。將測試集按照該GMF 經(jīng)驗模型進行風速反演并與 ECMWF 風速進行對比，對比結(jié)果散點圖如圖5 所示，均方根誤差為2.18 m/s。

圖4 DDMA 與風速的關(guān)系散點圖及擬合曲線

圖5 GMF 經(jīng)驗模型風速反演值與ECMWF 同比風速對比散點圖

3.2.2 Kalman 濾波模型反演結(jié)果

經(jīng)驗模型函數(shù)反演方法反演精度較低，同時在執(zhí)行嚴格的質(zhì)量控制后會損失大量數(shù)據(jù)樣本。本文采用上述Kalman 濾波模型，對同樣的測試集數(shù)據(jù)進行了處理，整體反演流程如圖6 所示。

圖6 基于Kalman 濾波模型的海面風速反演流程

在實際數(shù)據(jù)處理過程中，CYGNSS 的 DDMI的每個通道以1 Hz 的頻率采集數(shù)據(jù)，因此24 h 內(nèi)可得到86 400 個數(shù)據(jù)樣本，但并非所有數(shù)據(jù)樣本都可以用于風速反演。比如鏡面反射點經(jīng)過陸地時采集的數(shù)據(jù)，數(shù)據(jù)質(zhì)量太差（RCG＜3）導致的不可用觀測量，以及 DDMI 每隔 60 s 進行 4 s 的黑體負載校正期間導致的空白數(shù)據(jù)，此時通過 GMF經(jīng)驗模型無法得到風速的觀測值。本文利用ARIMA 模型及 Kalman 濾波模型的預測機制可對這2 種無效數(shù)據(jù)進行補充，從而提升數(shù)據(jù)利用率。

基于上述方法，對2018-09-24 的CYGNSS 的01號星的數(shù)據(jù)進行了反演，為了證明模型的可用性，圖 7（a）～圖 7（d）展示了基于 Kalman 濾波的 1 個階段的反演結(jié)果，圖 7（a）～圖 7（c）分別顯示了Kalman 濾波的前150 s、中間的150 s，及最后的150 s的結(jié)果。而圖7（d）顯示的是1 個階段內(nèi)的總體效果。

圖7 Kalman 濾波模型反演得到的局部時間序列

由圖7（a）～圖7（c）可得：相對于ECMWF同比風速的時間序列，GMF 經(jīng)驗模型反演得到的風速波動較大，同時存在空白數(shù)據(jù)及異常值；而通過 Kalman 濾波模型得到的反演風速值相對于GMF 更加收斂，減少了異常值，同時對短時間內(nèi)的空白數(shù)據(jù)進行了有效補充。由于連續(xù)的長時間空白數(shù)據(jù)會導致 Kalman 濾波模型預測能力變差，因此本文對空白數(shù)據(jù)樣本點的補充限制在連續(xù) 5 個樣本點以內(nèi)。

對2018-09-24 的24 h 的數(shù)據(jù)進行反演的結(jié)果如表1 所示。

表 1 2018-09-24 CYGNSS 01 號星 1 通道 24 h 數(shù)據(jù)反演結(jié)果

由表1 可知：利用Kalman 濾波模型進行反演的均方根誤差為2.08 m/s，較經(jīng)驗模型法有所提升，同時反演得到的風速數(shù)據(jù)樣本數(shù)量從經(jīng)驗模型的56 890 提升到了 63 175，提升了 11.04 %。與ECMWF 風速對比的散點圖如圖8 所示。

圖8 Kalman 濾波模型反演風速值與ECMWF 同比風速對比散點圖

通過圖5 與圖8 的對比可以發(fā)現(xiàn)，反演得到的整體數(shù)據(jù)質(zhì)量得到了提升，散點圖更加收斂，誤差較大的異常值有所減少。為了進一步分析模型的風速反演性能，本文將RCG 以間隔值為3 進行分段，統(tǒng)計各個分段的風速反演精度，結(jié)果如圖 9 所示。

圖9 各RCG 段風速反演誤差

由圖9 可以看出，當RCG 較小時的風速反演效果較差，同時 Kalman 濾波反演模型在 RCG 較小時具有比經(jīng)驗模型更好的反演效果。這是因為RCG 較小時接收機天線增益較低，得到的數(shù)據(jù)信噪比較低，從而導致風速反演結(jié)果較差；而Kalman濾波模型對該情形導致的異常值進行了校正，從而得到了更好的結(jié)果。

4 結(jié)束語

本文對星載GNSS-R 海面風速反演方法進行了研究。首先推導了用于風速反演的歸一化散射雷達截面及 DDMA 的計算流程；然后利用 ARIMA模型及 GMF 經(jīng)驗模型建立 Kalman 濾波模型，并進行風速反演。本文利用CYGNSS 星載數(shù)據(jù)對該Kalman 濾波模型進行了驗證，反演得到的均方根誤差為 2.08 m/s，較 GMF 經(jīng)驗模型反演方法的2.18 m/s 有所提升，與此同時數(shù)據(jù)利用率提高了11.04 %，進而提升了數(shù)據(jù)的空間覆蓋率。

致謝：感謝美國航空航天局（NASA）提供的CYGNSSS 公開數(shù)據(jù)。