姜祝輝,陳 建,沈曉晶
(1.西安測繪研究所,陜西 西安 710054;2.中國科學院大氣物理研究所 中層大氣和全球環境探測重點實驗室,北京 100029)
衛星遙感是獲取大范圍海面風場資料的主要手段[1]。Weissman等(1979)經分析發現合成孔徑雷達(Synthetic Aperture Radar, SAR)強度圖像和風速風向有相關性[2]。Alpers等(2011)指出SAR圖像中可以清晰的看到臺風的條紋,而該條紋與風向有關[3]。從此掀開了SAR反演海面風場研究的熱潮。散射計和SAR的雷達回波信號海面散射機理相同,所以可借鑒散射計反演海面風場的方法,利用SAR資料反演海面風場[4]。
基于散射計地球物理模型函數反演SAR海面風場方法(本研究稱之為“直接反演方法”)被絕大多數學者采用,該方法通過已知的風向、后向散射系數及雷達方位角和入射角等信息,利用散射計地球物理模型函數反演風速。散射計C波段垂直極化的地球物理模型函數已經發展得相當迅速且比較成熟,如CMOD4、CMODIFR2、CMOD5、CMOD7[5]等。楊勁松(2005)、張毅等(2007)對CMOD4地球物理模型函數進行了仿真,結果表明,后向散射系數在不同程度上依賴于風速、風向、入射角和極化方式[6-7]。Susanne等(2000)分別利用方位角交叉相關法和CMOD4地球物理模型函數法進行ERS-2 合成孔徑雷達波譜模式數據的風速反演,反演結果與ERS-2散射計進行比較后得出CMOD4地球物理模型函數法優于方位角交叉相關法的結論[8]。散射計利用不同視角后向散射系數通過最大似然法反演海面風場,He等(2005)提出了新的SAR反演海面風場的方法,風向模糊問題通過引入浮標等外部數據來解決[9]。艾未華等(2019)根據L/C雙頻共面SAR能夠獲取同一海域雙波段后向散射系數的特點,結合地球物理模型函數研究了一種新的SAR海面風場反演方法[10]。Fang等(2018)全面分析了不同極化方式的海面風場反演誤差[11]。李曉峰等(2020)系統總結了專用于SAR風場反演的地球物理模式函數[12]。方賀等(2018)評估了C-2PO和CMOD5.N地球物理模式函數SAR風速反演性能[13]。
直接反演方法需要將風向作為已知量反演風速[14]。該風向可以通過外部數據的引入來指定,比如用目標區域的浮標風向,或者散射計、模式預報結果中的風向來指定SAR反演風場中的風向,亦可以通過SAR圖像本身的條紋信息來指定風向[15-16]。Choisnard等(2008)指出直接反演方法會將風向誤差帶入后續的風速求解過程,30%的風向誤差將導致風速多達40%的不確定性[17]。Horstmann等(2000)指出,由于地球物理模型函數本身的對稱性,同樣的風向誤差導致的風速誤差在45°、135°、225°、315°附近最大[18]。
Portabella等(2002)[19]針對Horstmann等[18]指出的不確定性,首次提出了基于貝葉斯理論,將后向散射系數與通過數值預報結果得到的風向先驗估計引入到變分方程中,通過求取代價函數極小來確定最優風矢量的方法(本研究稱之為“變分方法”),代價函數求解方法是枚舉法,相當耗時。進而Choisnard等對該方法在不同背景風向情況下的反演結果進行了評估,代價函數求解方法是蒙特卡洛法,運算速度明顯提升[17]。Mouche等(2012)將多普勒頻移項加入上述代價函數,試驗證明在氣旋和鋒面天氣情況下反演結果更符合實際情況[20]。Jiang等(2017)提出利用阻尼牛頓法求取SAR反演海面風場變分方程代價函數極小的新解法,模擬試驗表明運算速度有極大提升[21]。
由于SAR資料刈幅窄,海上浮標資料稀疏,難以尋找足夠大樣本量的匹配數據,Jiang等僅找到了17組匹配數據進行了檢驗分析,初步得出了變分方法優于直接反演方法的結論[21]。考慮到這么少的數據量不能全面細致的刻畫在不同背景場條件下直接反演方法和變分方法的精度,本研究利用Jiang等提出的新解法[21],開展背景場對SAR反演海面風場變分方法風速反演結果影響分析。

(1)
式(1)中:v為待求海面風場,又稱分析場。σo為將v和SAR其他參數一同輸入到地球物理模型函數中計算得到的后向散射系數。Δσo為SAR后向散射系數標準差。Δv為數值預報結果中的海面風場標準差。J(v)的二階導數為?2J(v)。假定 ?2J(v) 連續,泰勒展開J(v)并保留前三項
J(v)=J(vB)+(?J(v)|B)T(v-vB)

(2)
因此
?J(v)=?J(v)|B+?2J(v)|B(v-vB)=0
(3)
假定 ?2J(v)|B是非奇異矩陣,其牛頓迭代公式為

(4)

假定β∈(0,1),σ∈(0,0.5),mk是滿足如下不等式非負實數
J(vk+βmdk)≤J(vk)+σβm(?J(v)|k)Tdk
(5)
所以利用阻尼牛頓法[22]迭代步長為
αk=βmk
(6)
式(4)可修改為
v=vB+αkdk
(7)
假定迭代初值為vB,可通過迭代得到最佳分析風場,詳細求解步驟見文獻[21]。
本研究首先簡單分析地球物理模型函數,研究風向、風速隨后向散射系數變化,然后著重分析背景場對變分方法風速反演結果影響。

直接反演方法是將已知的風向代入地球物理模型以求解風速,所以地球物理模型函數直接影響到風速反演精度。
當風速范圍為0~40 m/s,風向范圍為0°~360°,后向散射系數隨風速和風向的變化見圖1,其中實線代表風速,虛線代表風向。

圖1 后向散射系數隨風速和風向的變化Fig. 1 Variation of backscattering coefficient with wind speed and direction
當風向分別為0°、45°、90°時,隨著風速的逐漸增大,后向散射系數逐漸增大。當風速較小時,較小的風速變化將導致較大的后向散射系數變化;當風速大于20 m/s時,較大的風速變化將導致較小的后向散射系數變化。這就意味著如果采用直接反演方法,整體上講后向散射系數對較大的風速不敏感,風速反演誤差將會很大,這是變分方法要解決的問題之一。相比之下,風速較大時,風向為45°時較0°和90°時后向散射系數變化稍大些,也就是說如果采用直接反演方法,風向在45°附近時,風速反演結果比風向在0°和90°附近時要稍好一些。
當風速為3、8、25 m/s時,隨著風向從0°向360°變化,后向散射系數呈周期性變化。風速為3 m/s時,后向散射系數周期性變化的振幅小,且后向散射系數相對較小;風速為 25 m/s時,后向散射系數周期性變化的振幅大,且后向散射系數較大。風向為0°、180°、360°時,后向散射系數最大;風向為90°和270°時,后向散射系數最小;風向為45°、135°、225°、315°時,后向散射系數相對于風向的斜率最大,變化最劇烈。由于一切觀測均存在誤差,假定觀測風向有相同大小的誤差,真實風向為45°、135°、225°、315°時對應了最大的后向散射系數變化。直接反演方法是將已知觀測風向代入地球物理模型函數中計算風速,由于后向散射系數隨著風速的增大而增大,當真實風向為45°、135°、225°、315°時,不可避免的觀測誤差將導致后向散射系數更大的誤差,進而導致風速反演結果更大的誤差,這是變分方法所要解決的另一個問題。
考慮到圖1中風向為0°、45°、90°時后向散射系數和風速有不同的變化規律且風向為0°和90°時變化規律相似,假定背景場風向分別為45°和90°,背景場風速誤差分別為1 m/s和2 m/s,開展不同背景場條件下直接反演方法(虛線)風速誤差和變分方法(實線)風速誤差分析(圖2)。

圖2 不同背景場條件下直接反演方法(虛線)風速誤差和變分方法(實線)風速誤差分析Fig. 2 Wind speed error analysis by direct retrieval method (dotted line) and variational method (bold line) under different background wind fields(a)為背景場風速誤差1 m/s,背景場風向45°; (b)為背景場風速誤差1 m/s,背景場風向90°; (c)為背景場風速誤差2 m/s,背景場風向45°; (d)為背景場風速誤差2 m/s,背景場風向90°。
當背景場風速誤差為1 m/s,背景場風向為45°時,直接反演方法風速誤差呈近似對稱分布。背景場風向誤差的絕對值越大,直接反演方法風速誤差越大;背景場風速越大,直接反演方法風速誤差越大。當背景場風速約小于7 m/s時,直接反演方法風速誤差小于背景場風速誤差,隨著背景場風速的逐漸增大,直接反演方法風速誤差急劇增大,當背景場風速大于15 m/s時,部分直接反演方法風速誤差甚至超過了5 m/s。只要風速反演誤差大于背景場風速誤差,反演方法就不可用。背景場風向誤差為正值時,直接反演方法風速誤差為5 m/s的等值線投射到背景場風向誤差為負值時的等值線,對應的等值線約為-3.5 m/s的等值線,其他等值線亦有此規律,所以背景場風向誤差為正值時的直接反演方法風速誤差要小于背景場風向誤差為負值時的直接反演方法風速誤差。當背景場風向誤差在-3°~3°之間時,無論背景場風速為多大,直接反演方法風速誤差均小于背景場風速誤差。變分方法風速誤差隨著背景場風向誤差的增大而增大。當背景場風向誤差約小于15°時,隨著背景場風速的增大變分方法風速誤差略有增大;當背景場風向誤差大于15°時,隨著背景場風速的增大變分方法風速誤差略有減小。變分方法風速誤差最大的那條等值線出現在背景場風向誤差為20°附近,為1.2 m/s,略大于背景場風速誤差。整體上看,當背景場風速誤差為1 m/s,背景場風向為45°時,變分方法風速反演誤差要小于直接反演方法風速反演誤差,而在背景場風向誤差在-3°~3°之間時,直接反演方法風速誤差要小于變分方法風速反演誤差。
當背景場風速誤差為2 m/s且背景場風向為45°時,直接反演方法風速誤差與當背景場風速誤差為1 m/s且背景場風向為45°時的完全相同,變分方法風速誤差較背景場風速誤差為1 m/s時略大,趨勢相同,在大多數情況下小于背景場風速誤差。單純就變分方法風速誤差小于背景場風速誤差的情況看,背景場風向同為45°,背景場風速誤差為2 m/s時圖中的區域范圍要比背景場風速誤差為1 m/s時的要大。
當背景場風速誤差為1 m/s,背景場風向為90°時,直接反演方法和變分方法風速誤差均呈近似對稱分布。背景場風向誤差的絕對值越大,直接反演方法和變分方法誤差均越大;背景場風速越大,直接反演方法和變分方法誤差均越大。當背景場風向誤差在-12°~14°之間時,直接反演方法風速誤差小于背景場風速誤差。背景場風向誤差的絕對值較大時,直接反演方法風速誤差超過了3 m/s。整體上看,當背景場風速誤差為1 m/s,背景場風向為90°時,變分方法風速反演誤差均小于背景場風速誤差,體現出了優越性。與背景場風向為45°時類似的是,背景場風向誤差較小時,直接反演方法風速誤差比變分方法風速誤差要小些。
當背景場風速誤差為2 m/s且背景場風向為90°時,直接反演方法風速誤差與當背景場風速誤差為1 m/s且背景場風向為90°時的完全相同,變分方法風速反演誤差較背景場風速誤差為1 m/s且背景場風向為90°時略大,但均小于背景場風速誤差,體現出了方法的可用性。
當背景場風速誤差為1 m/s時,背景場風向為45°時直接反演方法誤差明顯比背景場風向為90°時要大得多,變分方法有效克服了直接反演方法的缺陷,整體誤差控制在了1.2 m/s以內。當背景場風速誤差為2 m/s時,也有類似的現象,不再贅述。
直接反演方法和變分反演方法的不同點有三個方面:一是輸入量不同。直接反演方法只需要輸入背景場風向,變分方法需要同時將背景場風速、風向作為輸入量;二是核心算法不同。直接反演方法中的地球物理模型函數算法簡單,而變分方法計算流程中的變分模型包含了地球物理模型函數及其反復迭代尋優,算法較為復雜;三是輸出量不同。直接反演方法分析場中的風向就是背景場風向,而變分方法對背景場風速風向都進行調整并遴選最優海面風場,分析場風速風向與背景場風速風向均不同。
本研究開展背景場對SAR反演海面風場變分方法風速反演結果影響模擬試驗分析。經分析,直接利用地球物理模型函數反演海面風速的直接反演方法存在兩大問題:一是當背景場風速較大的時候,SAR后向散射系數對風速不敏感,直接導致較大的風速反演誤差;二是當背景場風向存在誤差時,將導致后向散射系數更大的誤差,進而導致風速反演結果的更大的誤差。模擬試驗結果表明,當背景場風速較大時,變分方法風速誤差低于直接反演方法風速誤差。變分方法有效克服了直接反演方法不同背景場風向導致的風速較大的反演誤差的問題。而當背景場風向誤差較小時,直接反演方法風速誤差比變分方法風速誤差要小些。
后續將搜集SAR數據及同步觀測的浮標數據,開展變分方法核心參數的調整優化工作。