星載SAR 圖像方位角計算及其對風速反演影響評估*

李慧敏 胡登輝 林文明 王 臣 何宜軍

1（南京信息工程大學海洋科學學院 南京 210044）

2（自然資源部空間海洋遙感與應用重點實驗室 北京 100081）

3（中國科學院微小衛星創新研究院 上海 201203）

0 引言

自1978 年全球首顆衛星SEASAT-A 發射成功以來[1]，星載合成孔徑雷達（SAR）經過40 多年的快速發展，已成為地球表面常規監測的重要手段之一。特別是SAR 對于海面的觀測及關鍵物理參數提取，為物理海洋、生態環境及海氣相互作用等科學研究提供了重要數據[2-6]。國際SAR 衛星計劃的引領與持續投入主要以中國和歐洲為主，例如歐洲空間局的哨兵一號衛星（Sentinel-1）和中國的高分三號衛星（GF-3）、1 m C-SAR 等。其他衛星計劃也日漸增多，涵蓋不同波段、不同極化、不同軌道和不同數據獲取模式，為實現多顆SAR 衛星的星座組網觀測奠定了堅實基礎。星載SAR 衛星得到廣泛重視的原因在于其全天候、全天時、高分辨率的數據優勢，這對于地球表面持續監測、定位關鍵目標信息、響應環境高頻突發事件、聚焦亞中小尺度圈層科學問題等有切實效益。

SAR 衛星數據的應用涵蓋海洋學多個研究領域，雷達后向散射信號主要由海面布拉格波貢獻，其他尺度海洋現象通過改變或調制海面粗糙度，可以在SAR 圖像上留下不同散射信號特征[7]。在過去幾十年里，對于SAR 后向散射信號的處理算法逐漸成熟，發展了針對不同SAR 衛星的絕對定標、輻射定標、數據分級、參數反演等業務化流程，為SAR 數據分發共享與科學應用提供了諸多便利[8]。然而，SAR 數據使用過程中最常用到、又容易被忽略的是方位角信息，即SAR 圖像方位向在對應地球坐標投影下相對于地理北極的角度大小[2,7]。該角度準確與否直接影響地球方向性物理參數的反演精度，也間接影響反演過程中需要使用方位角信息的物理參數精度。隨著SAR 數據業務化服務的不斷提高，各級數據中經常提供多個角度參數，對其仔細篩選和應用尤為重要，往往決定著后續研究中是否會引入系統性誤差。

在與SAR 圖像方位角有關的眾多地球物理反演參數中，海面風場是最常見的一個，其算法原理與散射計海面風場反演相似，都是通過建立海面10 m 高度處近中性風場矢量與后向散射系數及入射角、極化、波段等雷達參數的經驗關系，構建地球物理模式函數（Geophysical Model Function，GMF）[9]。其中針對C 波段散射計的CMOD 是應用最為廣泛的GMF[10,11]，該系列函數一直在不斷的改進與提高。在CMOD4 函數建立不久，Alpers 和Brümmer 就將其應用到ERS-1 的SAR 圖像海面風場反演研究中[12]。基于SAR 反演數據，Young 與Mourad 團隊對海面風場精細結構進行深入探究，為揭示海洋大氣邊界層千米尺度動力過程開拓了思路[13,14]。隨后，Mouche等[15]針對星載C 波段SAR 后向散射數據開展了系統研究，在綜合考慮海面電磁散射特性及不同風力作用下的波浪情況，進一步改進發展了C-SARMOD 函數，可以較好地應用于哨兵一號SAR 數據的風場反演。

1 星載SAR 成像與坐標系統

1.1 星載SAR 成像幾何

需要特別指出的是，SAR 作為單天線觀測設備，無法像散射計一樣獨立提取風向，SAR 風速反演依賴外部風向輸入，同時因為GMF 表征的是雷達后向散射系數與相對風向的關系，即順風方向與SAR 雷達視向的夾角，計算該角度大小同時需要確定圖像方位角信息。為形象簡明地展示星載SAR 成像幾何關系，很多現有圖文資料[2,7]以圖1（a）所示為主，這類概念圖容易產生圖像方位角?IMG等于衛星飛行方向?SAT的誤導信息。相比之下，圖1（b）給出了考慮地球表面為曲面的SAR 衛星成像關系，可以看出?IMG與?SAT存在差異，兩者之差與星載SAR 信號處理、空間投影和坐標系轉換有密切關系。本文以歐洲空間局哨兵一號波模式（Wave mode, WV）SAR 數據為例，首先概述星載SAR 成像原理與坐標轉換數學關系，明確?IMG與?SAT無法等同的原因；然后給出基于地面控制點（Ground Control Point，GCP）的?IMG計算方法，并統計分析其與?SAT的關系；最后對比?IMG和?SAT用于SAR 風速反演的精度差異。為了最大程度量化SAR 圖像方位角的影響，實驗數據選取不受其他海洋大氣現象影響的SAR 數據[6]，風場參考資料選用歐洲中期天氣預報中心的ERA5 再分析數據，時間分辨率為逐小時，空間分辨率為25 km。本文內容將對SAR 圖像方位角的正確計算及其在海面風場反演中的影響起示范作用。

圖1 不考慮（a）和考慮（b）地球曲面的星載SAR 成像幾何概念。?SAT 和?IMG 分別為衛星飛行方向和SAR 圖像方位向與地理北極的夾角Fig. 1 Illustration of SAR imaging geometries without （a） and with （b） consideration of the curved Earth’s surface.The satellite heading ?SAT is the angle of its flighting direction relative to the North. The image heading ?IMG is the angle of geolocated azimuth direction of each pixel relative to the North

SAR 傳感器大多以線性調頻信號的方式發射電磁波，即發射波形幅度在脈沖時間τ內保持恒定，而瞬時頻率fi=krt隨著時間t以線性方式進行變化，其中kr為線性調頻頻率，kr與τ決定了波形產生帶寬Br=krτ。每一個τ之后是雷達接收回波窗口時間，并將接收到的回波強度存儲下來。所以SAR 在衛星飛行過程中，持續不斷的向與飛行方向垂直的徑向發射和接收電磁波，其數字信號處理涉及很多方面，與空間幾何相關的包括：確定從衛星到地球表面任何點目標的斜距距離；確定每一時刻雷達信號在衛星參考坐標系中的位置；根據精確衛星軌道參數，確定衛星坐標系與地球地理坐標系的轉換。根據SAR 信號數據的不同視角，可以歸納為圖像信號坐標系、空間衛星坐標系和地球地理坐標系，三者之間的關系及轉換如圖2 所示。

圖2 星載SAR 圖像信號坐標系（紅色）、空間衛星坐標系（黑色）與地球地理坐標系（藍色）的關系Fig. 2 Spaceborne SAR coordinate systems from image signal （in red） to satellite-centered （in black）and the Earth-centered （in blue）

1.2 圖像信號坐標系

圖像信號坐標系是指SAR 數字信號對于軌道時間軸的記錄方式，在右手笛卡兒坐標系中，定義x軸為雷達徑向方向，y軸為雷達沿軌方向，z軸為雷達信號強度，如圖2 紅色箭頭部分所示。x軸記錄每個雷達脈沖的回波時間，根據脈沖返回時間將其劃分為斜距單元；y軸記錄所有雷達脈沖時間，根據脈沖重復頻率劃分為方位向單元；z軸記錄點目標（τ,t）的反射強度S。這實際是對SAR 信號傳播軌跡的解碼和轉換，即對于某一個脈沖時間τm和沿軌時間tk，利用多普勒頻率進行方位向壓縮。在此過程中，可以確定徑向的單元分辨率為δr=c/（2Br），其中c ≈3×108m·s-1為光速；而方位向分辨率δa=da/2由合成孔徑算法所得[2]，其中da為雷達天線長度。

1.3 空間衛星坐標系

空間衛星坐標系（xs,ys,zs）是指原點在衛星質心的右手笛卡兒坐標系，用于較好實現圖像信號坐標系和地球地理坐標系之間的轉換，如圖2 黑色箭頭部分所示，其中zs是衛星質心指向地球質心的方向，xsys面與zs垂直且ys與衛星飛行方向（即沿軌方向）一致。對于任何一個圖像信號坐標（τm,tk），其與空間衛星坐標系（xs,ys,zs）關系為

其中，H為衛星軌道半徑，和分別為緯度和經度對時間的一階導數。對于圓形軌道，H對時間的導數為0，所以式（1） 可簡化為

1.4 地球地理坐標系

如圖2 中藍色部分所示，地球地理坐標系（xe,ye,ze）是指原點在地球質心的經緯度坐標系，對于某一經度為θe和緯度為?e的點， 其坐標為（recosθecos?e,resinθecos?e,resin?e），其中re為地球半徑。同樣對于經緯度為的衛星位置B可表示為（Hcoscos,Hsincos,Hsin?︿），其中H為衛星軌道半徑。（xe,ye,ze） 與（xs,ys,zs）的坐標轉換可以通過移動衛星位置B以及旋轉矩陣M來實現，即

其中，

通過式（2）和式（3），可以建立SAR 數字信號坐標（τm,tk）與地球地理坐標系（xe,ye,ze）的數學關系，詳細推導和原理參見文獻[16-18]。一個星載SAR 傳感器的發射到業務化運行，首先需要進行周密的在軌測試，其中一個基本步驟就是對SAR 圖像進行地理編碼/解碼和精確定位，目前最主流的方式是利用GCPs 確定衛星姿態及坐標系轉換中所使用的參數，所以在每一景SAR 圖像數據中，都會對應提供一定數量的GCPs。

2 合成孔徑雷達圖像方位角計算方法

從SAR 獲取的原始信號推導圖像方位角需要精確的衛星軌道、衛星姿態和SAR 成像算法信息，過程復雜且不利于實際應用推廣[19,20]。相比之下，SAR數據中提供的GCPs 經過嚴格推算和驗證，SAR 圖像中的像素點精確對應地球地理坐標系下的經緯度，可以用于SAR 圖像方位角的快速高效計算。對于哨兵一號 WV 數據，每一景SAR 圖像寬度和長度分別為約為20 km，含有77 個GCPs，如圖3 所示。該兩景圖像案例均位于陸地邊緣，海岸帶信息清晰可見，通過對比圖像陸地邊緣與背景地圖邊緣，佐證了SAR 圖像像素與地理經緯度是一一對應的，大量相關研究工作也已經量化了這一準確性[21,22]。

圖3 哨兵一號波模式SAR 圖像升軌（a）與降軌（b）案例。藍色點為GCPs 位置，藍色箭頭為利用GCPs 計算的圖像方位角?IMG，紅色箭頭為衛星飛行方向?SATFig. 3 Two Sentinel-1WV cases of ascending （a） and descending （b） showing the GCPs in blue points,?IMGcalculated using GCPs in blue arrows and ?SAT in red arrows

從圖3 可看出，GCPs（藍色圓點）相對于SAR 圖像是均勻分布的，圖像方位角可看作是同一方位向上兩個GCPs 的方向角，如圖3 中的藍色箭頭指示相鄰GCPs 的方向角。這一計算方式可追溯到大航海時代，用于查找地球上兩點之間的航向角，目前仍廣泛應用于航空、航海和陸地導航等領域[23]。當假設地球為規則圓球形時，經緯度為（θ1,?1） 和（θ2,?2）的A，B兩點方向角?B計算公式為

然而地球是以赤道半徑a=6378.137 km、極半徑b=6356.752314245 km和扁率f=（a-b）/a的橢球體，需要使用更為精確的計算方法，例如Vincenty[23]提出的實現算法如下：

步驟1計算相關參數。

步驟2迭代計算以下參數，直到λ（初值為L）變化非常小。

步驟3計算方位角?B。

利用該方法，計算圖3 所示兩個案例的圖像方位角，發現在整景SAR 圖像刈幅寬度內，基于不同方位向上GCPs 數據點對得到的方位角變化非常小（均方根小于0.2°），這是因為地球表面在較小范圍內可以近似為平面。圖3 中的?IMG為各自圖像方位向第一個GCPs 點對的計算結果，與數據中提供的衛星飛行方向?SAT相比（紅色箭頭所指），存在明顯的差異。具體而言，?SAT明顯偏離了SAR 圖像方位向，而基于GCPs 計算的?IMG與理論SAR 圖像方位角吻合。

3 哨兵一號波模式圖像方位角與衛星飛行方向對比

為深入探究哨兵一號波模式 SAR 數據的圖像方位角與衛星軌道方向角之間的偏差，本研究基于2016 年獲取的約70 萬景開闊大洋SAR 數據，統一計算了?IMG與?SAT的數值，并將兩者偏差的全球空間分布以2.5°×2.5°經緯度網格顯示，如圖4 所示。圖4（a）～（d）分別對應23°入射角（WV1）升軌（ASC）、23°入射角降軌（DSC）、36°入射角（WV2）升軌、36°入射角降軌情況，顏色代表平均差值大小。由此可以發現，衛星升軌軌道的?IMG-?SAT值在南半球為負值，在北半球為正值，且數值大小隨著緯度的升高而逐漸變大。對于兩個入射角，?IMG與?SAT兩者角度差在赤道都為0°，而在高緯度60°S 和60°N，WV1 分別約為-3°和3°，WV2 分別約為-6°和6°。降軌?IMG-?SAT隨緯度變化趨勢與升軌恰好相反，但角度差值的數值大小在相同緯度保持一致。結合圖1（b）可以推斷，在同一緯度下?IMG-?SAT隨傳感器入射角增大而變大，當入射角為0 時（即星下點），圖像將被壓縮成一條線，方位角與衛星飛行方向一致。同樣結合圖3 可知，哨兵一號為右視SAR 雷達傳感器，升軌時衛星在北半球由赤道向極地飛行，雷達距離向掃描的經度覆蓋區間逐漸較小，衛星高度處（H ＞re）相較于地球表面（re）的方位向更加偏離地理北極，導致?IMG＞?SAT。在南半球時，降軌可以看作是北半球的升軌，因而與升軌的?IMG與?SAT變化規律相反。

圖4 2016 年全球哨兵一號波模式SAR 圖像方位角?IMG與衛星飛行方向?SAT差值分布。 （a）和（b）分別為WV1（23°入射角）升軌與降軌。（c）和（d）分別為WV2 （36°入射角）升軌與降軌Fig. 4 Global maps of the difference between image heading ?IMG and satellite heading ?SAT for S-1 A WV1 ASC（a）, WV1 DSC （b）, WV2 ASC （c） and WV2 DSC （d） SAR data acquired in 2016

圖5 分別給出了WV1 升軌與降軌以及WV2 升軌與降軌?IMG-?SAT隨緯度的散點分布，從圖5 中可以看出，除極個別異常點外（地面處理系統升級導致）?IMG-?SAT與緯度的非線性關系在特定入射角和升軌/降軌時非常固定。研究表明，五次多項式對該關系的擬合誤差小于0.001°，有

圖5 2016 年全球哨兵一號波模式SAR 圖像方位角?IMG與衛星飛行方向?SAT的差值散點圖。黑色虛線為利用五次多項式對WV1 升軌與降軌和WV2 升軌與降軌數據的擬合，擬合參數見表1Fig. 5 Scatters of the difference between image heading ?IMGand satellite heading ?SAT along latitude for S-1 A WV SAR data acquired in 2016. Curve fits in dashed lines are performed using the fifth order polynomial function with coefficients listed in Table 1

擬合系數如表1 所列。利用式（6）和表1 系數，根據一級和二級數據中給定的衛星方位角?SAT可直接計算每一景哨兵一號波模式圖像的圖像方位角?IMG，避免反復讀取GCPs 導致的計算量，這對于實際地球科學研究應用十分有效，同時對于20 km 幅寬的波模式SAR 圖像，方位角在有限空間范圍內的變化可忽略不計。值得注意的是，該圖像方位角矯正偏差擬合公式只適用于哨兵一號 WV 數據，對于哨兵一號其他模式SAR 數據，需要重新利用相應的GCPs進行計算，且寬刈幅的SAR 圖像因方位角變化較大，不能使用單一數值替代。

表1 圖5 中五次多項式對WV1 升軌與降軌和WV2 升軌與降軌數據的擬合系數Table 1 Coefficients of the fifth order polynomial fit of the curves shown in Fig. 5

4 圖像方位角對海面風場反演精度的影響

SAR 圖像方位角和衛星軌道角之間存在的系統偏差，會直接影響與方向相關的地球物理參數提取，例如海浪方向、海面風向、海洋內波傳播方向等。從SAR 圖像中得到的方向信息一般以SAR 圖像坐標系為參考，也就是以SAR 圖像方位向為0°角，當需要轉換為以北極為0°角的地球地理坐標系時，需要輸入SAR 圖像方位角大小。此時若默認?IMG=?SAT，就會引入上述如圖4 和圖5 所示的系統誤差，該誤差隨緯度增加而變大，最大可達10°以上；且對于升軌和降軌軌道，該誤差數值互為相反數，這對于分析某一地球物理現象的日變化有不可忽視的影響。同時由圖4和圖5 可知，當不區分升軌和降軌而進行空間平均時，?IMG與?SAT之間的偏差會相互抵消。由于目前絕大多數與方向相關的研究以統計為主，且物理參量的均方根誤差在30°左右，導致衛星軌道角引起的誤差尚無法得到報道，圖像方位角的準確計算方法也沒有明確指出和重視。

特別需要說明的是，SAR 圖像方位角引起的偏差會進一步累積影響相關物理參數的反演。以海面風速反演為例，一般已知SAR 后向散射系數及入射角、極化、波段等雷達參數，但相對風向需要外部絕對風向和?IMG共同確定，然后再結合經驗GMF 函數經過最小二乘法得到風速數值。毫無疑問，當?IMG角度值不準確時會使相對風向存在誤差，進而導致風速反演的誤差。圖6 給出了2016 年哨兵一號波模式 SAR數據利用?IMG和?SAT進行風速反演的偏差，顏色代表2.5° ×2.5°經緯度網格內的偏差均值。為了量化圖像方位角對SAR 風速反演的影響，使用同一個外部輸入風向，即對于每一景SAR 圖像，在時間和空間上將ERA5 再分析數據的10 m 高度風矢量進行插值，以獲取最佳匹配數據，同時使用同一個C-SARMOD GMF，這樣就保證了對于同一個SAR 數據，風速反演的唯一變量是圖像方位角。總體而言，不同入射角和軌道方向反演的兩者風速偏差大小有所差異，36°入射角的風速偏差大于23°入射角，這是因為大入射角下的?IMG與?SAT偏差更大，引起的相對風向偏差也更大；然而不同于圖4，風速偏差的正負對于升軌和降軌并不是相反的，這是因為風速反演數值與相對風向緊密相關，而非簡單的線性關系。

圖6 利用?SAT 和?IMG對2016 年全球哨兵一號 WV SAR 數據進行風速反演的偏差。（a）和（b）分別為WV1（23°入射角）升軌與降軌。（c）和（d）分別為WV2 （36°入射角）升軌與降軌。反演算法統一使用C-SARMOD GMF 和ERA5 再分析數據的10 m 高度處風向。正值和負值分別代表高估和低估Fig. 6 Global mean biases of SAR retrieved sea surface wind speed using between the satellite heading and image heading. （a） and （b） are WV1 ascending and descending, and （c） and （d） are WV2 ascending and descending, respectively. The C-SARMOD GMF is used with wind direction inputs of ERA5 U10.Positives and negatives here mean overestimates and underestimates, respectively

從圖6 可以看出，使用?IMG=?SAT所引起的風速反演誤差是不能忽視的。哨兵一號波模式的兩個入射角風速誤差在全球空間分布對于升軌或降軌是相似的，升軌時除赤道附近、北半球中高緯和30°S 附近外，總體為負值；降軌時除南半球中高緯外，總體為正值。然而對于同一入射角，升軌與降軌的風速偏差并不能抵消，且大入射角下的風速誤差更大，可達±0.5 m·s-1。因此，正確使用SAR 圖像方位角對于海面風速反演影響顯著，尤其是對于探究風速的早晚差異性。另外，哨兵一號波模式作為當下唯一在開闊海洋持續收集數據的SAR 衛星，可以提供超高分辨率的海面風場數據，對于高分辨率的數值模式同化有巨大應用潛力。為保障SAR 海面風場的反演精度，明確SAR 圖像方位角的定義、計算、使用及其對風速反演精度的影響，也是十分必要的[24,25]。

基于?IMG和ERA5 的外部風向，進一步評估哨兵一號波模式SAR 數據全球開闊海洋風速反演的精度，參考資料為ERA5 10 m 高度處的風速，如圖7 所示，其中顏色為歸一化到0～1 的數據密度。對于兩個入射角和升軌/降軌，SAR 反演風速與ERA5 風速的平均偏差在0.35 m·s-1以內，表現為略微高估，標準差在1.3 m·s-1左右。相比于全球散射計風速1.5 m·s-1的總體統計誤差[24]，開闊海洋的波模式SAR 反演風速具有很高的精度，這不僅僅因為可以通過機器學習智能分類[6]，最大程度上去除了SAR 數據中其他海洋大氣現象的影響，也得益于對SAR 數據處理的嚴格要求。鑒于SAR 高分辨率的優勢特點，其在提供高分辨率海面風場的潛力不可忽視，對當前海氣界面千米尺度動力過程的研究具有重要意義。

圖7 利用?IMG的2016 年全球哨兵一號波模式 SAR 數據風速反演精度。（a）和（b）分別為WV1 升軌與降軌。（c）和（d）分別為WV2 升軌與降軌。對比風速為ERA5 再分析數據Fig. 7 Accuracy of wind speed retrieved using ?IMG from the Sentinel-1 WV SAR data in 2016. （a）～（d） are for WV1 ASC, WV1 DSC, WV2 ASC and WV2 DSC, respectively. The reference wind speed is from ERA5 reanalysis product

5 結語

星載SAR 作為海洋觀測的重要手段之一，在過去幾十年的發展和應用中展現出巨大的潛力，SAR 獲取的高分辨率海面粗糙度圖像，可以為海洋與大氣千米尺度物理過程提供數據保障，有效服務地球科學國際研究。然而在涉及SAR 圖像方位角的物理參數提取和反演中，一直存在疑問和誤用的情況。對此，本文基于2016 年哨兵一號全球波模式數據，詳細闡述了SAR 圖像方位角的推算方法，及其與衛星軌道方位角的區別；在推導SAR 空間坐標轉換的過程中，厘清了SAR 成像中的幾何變化關系；提出利用GCPs 計算SAR 圖像方位角的方法，并利用個例和統計數據進行對比分析，發現?IMG與?SAT存在系統偏差，且該偏差與入射角和軌道方向有關；在此基礎之上，進一步分析了使用?IMG與和?SAT對海面風速反演的影響，通過量化兩者之間的偏差大小，指出?SAT會導致最大±0.5m·s-1的風速反演誤差。本文研究明確了SAR 圖像方位角的計算方法，指出了其在海面風場反演中的重要性，可以便利地遷移應用于其他SAR 衛星數據的業務化算法。同時，本文詳細闡述了SAR 圖像方位角與衛星軌道角的概念和兩者差別，并以實際觀測數據進行對比展示，這對于的SAR 的入門使用群體具有明確的示范作用。

SAR 圖像方位角不僅對海面風場反演有影響，其對任何方向性的SAR 物理參數提取都至關重要。當前哨兵一號全球波模式數據使用中，尚沒有文獻報道該如何計算圖像方位角，已有的波模式數據（例如海浪和海面風場）都是默認使用衛星軌道方位角進行替代，這一近似導致的偏差是不可忽略的。本文明確了圖像方位角的正確計算和使用方法，彌補了業界對于這一參數認識和重視的不足，對后續海洋SAR 圖像數據處理和科學應用具有重要的參考意義。

致謝本文使用的哨兵一號波模式數據由歐洲空間局提供，該數據可通過網站https://www.esa.int/獲取。歐洲中期天氣預報中心再分析風場數據由ECMWF 氣象存檔與檢索系統提供，該數據可通過https://www.ecmwf.int/公開獲取。