基于GMF方法的捕風一號風速反演分析

范東棟，盧敏健，陳趁新，高涵，尉昊赟,*

1. 航天東方紅衛星有限公司，北京 100094

2. 清華大學 精密儀器系，北京 100084

1 引言

海面風是海洋各種運動的主要動力來源，是形成海面波浪的直接動力，維持著區域與全球氣候；同時海面風也是氣象預報的必要參數，直接影響海上航行、海洋工程、海洋漁業等。

傳統的海面風探測手段，如海洋浮標、船舶、海洋站等，存在觀測點少、觀測區域受限、難以實現大范圍和惡劣天氣條件下的實時有效觀測等問題。為此，基于微波散射計[1]、激光測風雷達[2]、全球導航衛星系統-反射測量(global navigation satellite system-reflectometry，GNSS-R)[3]等技術的星載海面風探測技術得到了廣泛關注。其中，GNSS-R技術直接利用衛星導航系統全球覆蓋的微波信號源，具備全天時全天候觀測、星載設備體積小成本低等優勢，發展前景良好，應用潛力巨大[4]。2000年，Zavorotny和Voronovich基于雙基雷達方程[3]，提出了二維時延-多普勒模型。進一步，研究者結合不同海面反射條件下GNSS-R波形不同的特點，利用波形的不同特征，如平均功率、前沿斜率、后沿斜率等，構建波形特征參量與風速的回歸模型，實現風速反演。2003年，UK-DMC衛星進行了首次星載GNSS-R探測試驗，驗證了在軌海面風探測應用的可行性[5-6]。之后，國外相繼啟動了多個星載GNSS-R計劃，其中最為著名的是2016年 NASA發射的用于熱帶氣旋監測的低軌衛星星座——颶風衛星導航系統(cyclone global navigation satellite system，CYGNSS)[7]。CYGNSS的發射，極大地推動了GNSS-R風速反演技術的發展，C. Ruf 等人基于歸一化雙基雷達散射截面(normalized bistatic radar cross-section，NBRCS)和前沿斜率兩個特征量，分別構建了經驗性地球物理模型函數(geophysical model functions，GMF)，并通過高低風速不同建模及反演風速最小方差融合等方法，實現了風速反演性能的提升并向高風速段邁進。特別值得一提的是，對颶風核心處的風速反演結果與NOAA P-3機載測量數據具有較好的一致性[8-9]。近年來，隨著人工智能和深度學習技術的發展，新型反演算法也得到了快速發展。在各級算法發展的大力支持下，CYGNSS星座持續提供L1數據及風速反演L2、L3數據[8-9]，并將反演數據產品延拓到土壤濕度[10]等。在CYGNSS取得一系列成果和國內GNSS-R技術走向成熟[11-13]的共同推動下，2019年6月5日，中國成功發射捕風一號A/B 衛星，開展GNSS-R技術在軌試驗。

本文基于GMF方法，實現了中國首顆GNSS-R衛星捕風一號數據的風速反演及初步性能評估，可為后續衛星星座設計提供參考。

2 風速反演原理

捕風一號衛星接收機與CYGNSS接收機類似，在軌接收到的前向散射信號主要來自導航衛星、地表、探測衛星三者在滿足以地表為“鏡面”的鏡面反射幾何條件下的反射信號，如圖1所示。其中，T、R分別為GNSS衛星和GNSS-R衛星；XYZ、XTYTZT、XRYRZR分別為地心慣性坐標系、GNSS衛星本體坐標系和GNSS-R衛星本體坐標系。S為鏡面反射點，以S為原點，建立鏡面反射點本體坐標系，其中Z0軸指向鏡面反射點法線方向，Y0軸位于OTR平面內，X0軸垂直上述平面。此外，VT，VR分別為GNSS衛星和GNSS-R衛星的速度矢量；Re、HT、HR分別為地球半徑、GNSS衛星高度和GNSS-R衛星高度；RT、RR分別為GNSS衛星和GNSS-R衛星到鏡面點的距離；θ為鏡面點入射角。

圖1 GNSS-R觀測幾何示意Fig.1 Geometry of GNSS-R measurement

根據Zavorotny和Voronovich模型(ZV模型[3])，導航衛星、鏡面點、探測衛星三者位置與速度關系在鏡面點附近形成橢圓環帶的等延遲面(等傳輸路徑長度)和呈曲線狀分布的等多普勒條帶(等相對運動速度)。按相應時延-多普勒二維分割得相關功率信號，形成延遲-多普勒映射圖(delay-Doppler mapping，DDM)。每一DDM單元對應的相關功率可表示為：

(1)

進一步考查式(1)中NBRCS (即σ0)項，圖2給出了按照捕風衛星在軌參數及鏡面點入射角20°條件仿真得到的σ0分布情況。可以看出在風速1 m/s時，以鏡面點為中心20 km×20 km仿真區域內的σ0為曲面分布，此時洋面與海面風場之間的耦合也較弱；而在風速2 m/s時，σ0趨于一個恒定值。更高風速的仿真同樣表明，σ0在鏡面點區域附近保持恒定。因此，GNSS-R分析建模中，可將式(1)中σ0項在整個鏡面點附近區域視作單一取值置于積分項外，這樣，通過對公式中其他參數的定標及校準計算即可由觀測量導出與洋面風速直接相關的鏡面點區域σ0。圖3給出了仿真得到的不同鏡面點入射角時，NBRCS與風速的對應關系，可以看到，在給定鏡面點入射角的情況下，兩者呈一一對應關系。正是基于這一內在關聯，可通過建立NBRCS和風速的回歸GMF模型，實現風速反演。

圖2 不同風速下NBRCS仿真Fig.2 Simulation of NBRCS with different wind speed

圖3 不同入射角時NBRCS與風速關系仿真結果Fig.3 Simulated relationship between NBRCS and wind speed for different incident angles

3 風速反演流程與建模

3.1 反演分析流程

為了便于反演性能的驗證評估，本文采用CYGNSS風速反演的同一算法，該風速反演算法以文獻[14]為基礎，利用在軌觀測數據和約定風速參考構建經驗性參數化GMF回歸方程，實現風速反演。圖4給出了風速反演的基本流程。其中地理數據用于篩選出洋面數據，以全球陸地1 km級基準海拔數據(global land one-kilometer base elevation，GLOBE[15])為基礎，并將陸面區域延拓15 km以避免近岸效應的影響。參考風速集選用歐洲中期天氣預報中心(European center for medium-range weather forecasts，ECMWF)ERA-Interim[16]小時級再分析海面風速數據。通過將ERA-Interim再分析數據的空間位置二維線性差值和時間一維線性差值與觀測數據進行時空匹配，并將匹配后的風速作為建模的參考風速。

圖4 風速反演流程Fig.4 Flow chart for wind speed retrieval

3.2 反演輸入觀測數據集

捕風一號試驗衛星1級數據包含15個數據變量：觀測時間參數(week，second，time_utc)；觀測時刻捕風衛星位置速度參數(remote_ecef)；導航衛星類型、PRN和位置速度(gnss_type，gnss_prn，gnss_ecef)；鏡面點位置速度和反射信號角度(reflect_ecef，elevation_angle)以及觀測特征量降采樣DDM、功率校準DDM、信噪比和NBRCS(ddm_0，ddm，ddm_snr和NBRCS)。圖5給出了典型的功率校準DDM，呈典型的馬蹄形分布。NBRCS數據是用于反演分析的基礎觀測特征量，每個觀測點對應一個7 delay×5 Doppler的二維數據，其中delay為-0.5 ～ 1 chip，0.25 chip分辨間隔；Doppler為-1.0～ 1.0 kHz，500 Hz分辨間隔。

圖5 典型的功率校準DDMFig.5 Typical power calibrated DDM

1級數據產品中A星和B星日均觀測點數均約為6 000點，最高可到8 000點。反演建模中以2019年7～9月的數據作為建模數據集。圖6給出了1級數據觀測的洋面空間分布圖。可以看出，捕風一號可實現±45°范圍內的有效覆蓋。圖中顏色深淺代表了觀測頻次的大小，可見在當前的軌道設計下，在南北緯30°～45°之間觀測的頻次最高。相比較于CYGNSS，捕風一號可以獲得更高緯度的觀測數據。

圖6 2019年7月全月衛星數據觀測區域及頻次分布Fig.6 Satellite observation region and frequency distribution in the full month of July, 2019

3.3 基于GMF方法的經驗模型構建

按理論模型，逐像元的NBRCS表征的是圖2中的洋面特征量，在鏡面點附近近似為一常數，但實際探測信號由于定標及噪聲等因素，逐像元的NBRCS存在明顯的數據波動。為此，綜合考慮區域NBRCS均值及標準差，選取3 delay×5 Doppler區域的NBRCS均值作為反演建模分析觀測特征量。

分析捕風一號1級數據中影響GMF建模的敏感性因素發現：除鏡面點入射角因素外，不同衛星和不同增益天線明顯影響NBRCS和風速之間的關系。圖7給出了A星天線1和天線2在2個月時段內的數據對數密度圖，可以發現在同樣的參考風速段內，天線2得到的觀測數據更為分散，直觀可視。數據的分散性差異會直接影響GMF建模的模型參數值及模型參數誤差大小。

圖7 2019年8～9月A星不同接收天線NBRCS與參考風速關系的對數密度圖Fig.7 Log (density) scatterplot of NBRCS of satellite A antennas in August and September, 2019 vs. reference wind speed

為了直觀展示由此帶來的差異，圖8給出了A星和B星不同天線數據在20°±2.5°鏡面點入射角范圍內的提取曲線。其中曲線上的數值點為各參考風速對應的實測NBRCS數據點均值。需要特別指出的是，當高風速對應NBRCS均值大于低風速對應值時，以低風速值替代，以保持曲線的單調性。對比圖8可以看到，兩星的天線2數據明顯比天線1數據小，且B星兩天線間數據的分離更為顯著。據此，建模過程中，除了理論上需要按入射角分開建模外，還需要根據數據實際情況，對不同衛星、天線建立不同的GMF模型參數。

圖8 A星、B星不同接收天線NBRCS與參考風速關系Fig.8 Relationship between NBRCS and reference wind speed of satellite A and satellite B with different antennas

進一步，分析觀測數據SNR和參考風速關系，如圖9所示，發現高信噪比數據對應低風速數據，而高風速數據則整體上處于低信噪比狀態，SNR>4的數據表現出較好的綜合曲線特征，但易缺失較大比例高風速數據。因此，建模過程中，按SNR對數據進行分組構建不同的GMF模型參數，能更好地體現觀測數據特征。

圖9 DDM SNR與參考風速關系密度圖Fig.9 Log (density) scatterplot of measured DDM SNR vs. reference wind speed

綜上，考慮到現階段數據的上述敏感性特征，本論文在CYGNSS按鏡面點入射角建模的基礎上，引入更多的建模細分，按衛星(A星/B星)、天線(1號天線/2號天線)和DDM SNR等影響因素進行分類篩選。實際建模過程，還包括以下幾方面特色處理：1)為避免不同風速段建模數據量差異造成的模型對高風速段不敏感性，以參考風速值給定窗口(通常±0.5 m/s)內建模數據的NBRCS均值為該風速下特征值的表征，從而獲得離散的特征數據點。2)為了進一步提升高風速段數據在建模中的有效性，在離散數據構建GMF參數化模型時，中低風速(<12 m/s)和中高風速(>10 m/s)獨立構建模型， 并以10～12 m/s區域為過渡區，通過兩段擬合參數銜接獲得最終的GMF模型數據。邏輯上，建模依賴的地面參考風速具有更小的誤差，所以在模型建立過程，選用參考風速為自變量，觀測特征量為應變量。中低風速段和中高風速段選用的GMF擬合函數分別如式(2)及(3)所示：

(2)

(3)

式中：σ0為NBRCS；ai，bi為模型參數；uref為參考風速。

圖10給出了其中一個條件組合下的建模示例，其中對數密度圖反映的是觀測數據的實際分布，粉色曲線為提取出的離散的特征量統計表征值，黑色曲線為擬合得到GMF參數化模型曲線。圖11給出了以A星天線1數據為例提取的不同角度離散特征值和GMF參數化模型曲線。總體而言，這種差異性與理論分析基本一致，且與CYGNSS相關數據處理理論與算法分析中反映的特征基本相符，說明可以以此進行捕風一號數據的風速反演。

圖10 GMF參數化曲線構建Fig.10 Diagram of GMF parameterized modeling

圖11 A星天線1按鏡面點入射角劃分的GMF參量化曲線Fig.11 Discrete and continuous GMF curves with different specular incident angles for antenna 1 of satellite A

4 風速反演結果分析

利用模型構建導出模型參數，對捕風1級觀測數據進行風速反演分析，此處以GNSS-R同類型觀測模式CYGNSS的結果作比對來評估風速反演的性能。圖12給出了反演獲得捕風一號A、B雙星及CYGNSS星座8顆衛星在2019年8月6日同一天內的觀測結果。其中捕風一號數據采用NBRCS數據反演得到2級軌道數據，CYGNSS數據為融合NBRCS和LES觀測量反演結果的3級網格風速產品。比對數據可以很明顯看到兩者數據具有很好的全局一致性，如日本東側太平洋區域、馬達加斯加島附近有明顯的低風速區域，孟加拉灣和西風帶具有明顯的高風速區域。需要指出的是，因為比對數據級別不同，時空對應關系也不能完全匹配，因此圖12僅是一個定性的比對。對捕風一號衛星而言，可用于反演的反演特征量目前僅為NBRCS，反演數據質量分析也存在多方因素受限，數據產品中離散出現的高風速數據點較CYGNSS偏多，也制約了更為詳盡的定量比對。后續尚需要提取更多的數據特征量和反演算法的持續改進。下文將主要從反演數據的統計特性上與CYGNSS作比對分析，間接驗證本文的反演水平。

圖12 捕風一號和CYGNSS同日(2019-08-06)觀測反演數據比對Fig.12 Comparison of the retrieved wind speed of BF-1 and CYGNSS data in the same day (2019-08-06)

為了進一步展示本反演結果的統計特性，在此展示與CYGNSS ATBD[17]文檔相對應的兩組結果，并與之相對比。

圖13給出了本文和CYGNSS ATBD利用NBRCS反演得到的結果與參考風速之差隨鏡面點入射角的變化特性對數密度圖。其中顏色越暗紅，表明這類結果出現的頻次越高。圖13(a)為捕風一號反演結果，(b)為CYGNSS L2 ATBD的結果。可以看到，兩者反演結果的偏差趨勢基本一致。兩者反演風速偏差大部分落在±5 m/s的區域內，且兩者反演風速大于參考風速的情況明顯偏多。分析數據可以發現，這部分多出現于低信噪比和高風速情況。圖10中目前的1級數據中低參考風速仍存有不少偏小的NBRCS值即直接造成大負偏差數據的產生。比較圖13(a)和(b)還可以發現，捕風一號衛星反演數據在入射角10°～20°范圍偏差相對偏大，這可能和載荷與衛星太陽翼間的相對關系有關。

圖13 反演偏差隨鏡面點入射角分布特性對比Fig.13 Comparison of retrieval deviation vs incidence angles of BF-1 NBRCS and CYGNSS result in ATBD

圖14給出了不同參考風速對應的風速反演偏差特性，圖中顏色表示數據的對數頻次，其中(a)圖為捕風一號，(b)圖為CYGNSS L2 ATBD的結果。可以看到，與圖13結果類似，兩者反演結果的偏差趨勢基本一致，在風速較低的情況下，反演偏差對稱性較好，偏差比較集中；隨著風速的提高，風速偏差分布也逐漸增大。這與兩方面因素密切相關：1)NBRCS特征觀測量變化靈敏度降低，無論圖3理論分析還是圖10實際數據特征，都表明在高風速時NBRCS特征觀測量隨風速變化趨于陡直，變化不明顯；2)觀測數據信噪比降低，高風速信號散射明顯，NBRCS特征量值較小，從而造成探測信噪比下降，如圖9所示，普遍低于5 dB。上述兩個因素造成NBRCS特征量相對誤差明顯，進而體現到風速偏差分布顯著增大。對比(a)圖和(b)圖可以看出，本文針對捕風一號數據特征，引入比CYGNSS更多的建模分類因子，以增加模型復雜性和犧牲模型適用性為代價，得到反演結果的整體分散性略優。后續針對1級數據的高精度定標，特別是在低信噪比區域，減小定標誤差，將是提升反演性能重要切入點之所在。

圖14 反演偏差隨參考風速的分布特性對比Fig.14 Comparison of retrieval deviation vs reference wind speed of BF-1 NBRCS and CYGNSS result in ATBD

5 結論

本文以捕風衛星在軌觀測獲得的1級數據為基礎，開展了風速反演實踐，并與CYGNSS相關結果進行了比對分析。主要結論如下：

1)通過對觀測特征量敏感性分析，表明捕風一號當前1級數據對數據源衛星、觀測天線、鏡面點入射角、DDM SNR等參量具有顯著的敏感性差異。

2)結合捕風一號數據敏感性，對數據進行了分類和篩選，分別構建參考風速與特征量的GMF模型，實現風速數據的反演生成。

3)對反演結果的全球分布特征、隨鏡面點入射角和參考風速的統計特性與CYGNSS L1 ATBD數據進行了定性比對分析，表明兩者一致性良好。

本文工作及國內同行對捕風一號數據開展的研究均表明捕風衛星試驗的成功[18]。反演實踐也發現隨著1級數據定標性能的改進和更多特征量的引入，以及2級反演新方法的采納和多特征量反演結果的融合，捕風一號能獲得更有價值的數據產品。本實踐相關成果可為后續星座發展提供參考。