基于SAR 多普勒質心頻移的海面流場迭代反演算法

秦艷萍，范陳清，張玉濱*

(1.中國海洋大學 信息科學與工程學院，山東 青島 266100；2.自然資源部第一海洋研究所，山東 青島 266061)

1 引言

海流是海洋重要的運動形式之一，與海上污染物的擴散、海上軍事活動、水產養殖和全球氣候變化等密切相關[1]，所以研究其特征和變化規律具有重要的科學意義。遙感技術的快速發展使得大范圍、高精度的探測海洋信息成為可能[2]。目前，利用遙感技術反演海流的方法主要有光學遙感和微波遙感。光學遙感主要是根據示蹤物（葉綠素等）的位移時間關系來反演流速，其易受天氣的影響[3–4]；微波遙感不受天氣的影響，其海流反演手段主要有：高度計[5]、地波雷達[6–7]、多普勒雷達散射計[8–11]、合成孔徑雷達（Synthetic Aperture Radar，SAR）[12–16]等。高度計反演海流覆蓋范圍廣，但其流場產品分辨率較低，不適合小尺度海流的測量[5]；地波雷達雖具有反演分辨率高、全天候等優勢[6]，但其主要是岸基測量，覆蓋區域有限；多普勒雷達散射計可實現二維測流[8]，但因為其是真實孔徑雷達，流場產品分辨率較低；而合成孔徑雷達則具有高分辨率、全天候、全天時等優勢[12]，使其成為反演海流的新手段。目前，利用SAR 技術反演海流主要有兩種方法：順軌干涉（Along-Track Interferometry，ATI）法和多普勒質心頻移（Doppler Centroid Anomaly，DCA）法。順軌干涉法基于兩幅沿軌天線分別獲得的SAR 復圖像的相位差反演流場[17]；多普勒質心頻移法則根據海表層運動產生的多普勒質心頻移與其徑向速度的關系反演海表面流速[18–21]。相比DCA 方法，ATI 技術反演精度、分辨率更高[22]，但其成像條件苛刻，數據不易獲取，所以不具有普適性。雖然DCA 方法流速產品分辨率較低，但其適用數據豐富且容易獲取，如Radarsat-2、Envisat ASAR 和Sentinel-1 等數據都可用于DCA 方法反演海面流場，具有業務化觀測全球海洋表面流的潛力。

DCA 方法中由多普勒質心頻率異常反演的地距多普勒速度不僅與海表面流場有關還與海表面風場有關。2004 年Chapron 等[21]利用Envisat ASAR 數據分析了全球海洋多普勒測量數據，發現多普勒頻移與海面風場有很高的相關性；后來Chapron 等[23]又利用Envisat ASAR 數據基于多普勒質心頻率異常提出了一種簡單的模型：UD≈γU10+UC（γ為風貢獻因子），以表示海面風場和流場對多普勒速度的影響，并發現當入射角為23°且在中風和海態完全發育時，風貢獻因子 γ=0.3。國內學者在基于DCA 法反演海面流場并去除風場影響方面也做了大量研究。楊小波[24]利用多普勒質心頻移法反演海流時，分析了海表面風場對流場反演的影響，發現風對流場的分布和結構有一定的影響，并給出了SAR 成像時刻研究區域風場產生多普勒速度的關系式：UDW=0.03U10；侯富城等[12]利用多普勒質心頻移法提取了內波產生的海表面流，并基于Chapron 等[23]提出的模型，指出在中等風速下，風貢獻因子 γ取值為0.2～0.25。雖然Chapron 等[23]提出的模型給出了風場和多普勒速度的關系，可以很好地修正風場產生的反演誤差，但只給出了風貢獻因子 γ在特定雷達頻率、入射角下的經驗值。風貢獻因子 γ與雷達頻率、入射角以及雷達天線的偏航角等多個參數有關[23]，一般不能由簡單的線性關系描述，并且對流場反演結果的精度影響很大，所以需找到正確估算風貢獻因子 γ的方法，以提高流場反演精度。

針對以上反演過程中去除風場影響存在的問題，本文提出了一種基于M4S 模型的迭代方法，用于流場反演。M4S 是由德國科學家Romeiser 等[25]開發的海面微波成像仿真模型，其核心程序模塊主要包括海表面微尺度波波高譜計算模塊和雷達海面成像仿真模塊。在給定風場和流場的情況下，可仿真海浪譜、海面SAR 復圖像（強度、相位）以及多普勒譜等信息。因此，基于M4S 模型并結合弦截下山法可估算并去除風場的影響，反演海面徑向流速。首先利用頻譜擬合法從SAR 原始復數據中獲得實測多普勒中心頻率；再利用SAR 數據頭文件信息估算衛星和地球表面的相對運動產生的預測多普勒中心頻率，繼而得到多普勒質心頻率異常值并將其轉換為地距多普勒速度；然后基于M4S 模型，利用本文提出的弦截下山法迭代反演SAR 局部區域的流場，估算風場對多普勒速度的影響，繼而反演整幅SAR 圖像的海面徑向流速；最后通過與實測海面流場數據比對，驗證反演算法的有效性。

2 數據源

本文采用的實測海表流場數據和SAR 數據來自作者課題組組織的星地匹配實驗，用安德拉海流計現場采集海表流場數據，并與Radarsat-2 SAR 衛星數據作時空匹配。采用安德拉海流計測量海表面流時，將海流計懸掛于浮體之下，使其處于海表面以下約0.4 m處，在船只停航時，用繩子牽住浮體，并在船尾將浮體和海流計放置于離船200 m 的海面，記錄停止放繩的時間和收繩的時間，此時間段內為有效數據。所用安德拉海流計型號為SEAGUARD，這是一種自記錄海流計，具有二維流場測量能力，測量流速的分辨率為0.1 mm/s，平均誤差為±0.15 cm/s，測量流向的分辨率為0.01°，誤差為±5°?，F場測量海表面流時，實驗海流計位于海表面以下約0.4 m 處，測量的是近表面流速，通過SAR 數據測量得到的流速是海表層流速。因為兩者測量的深度基本匹配，所以可用此實驗海流計所測的實測流場數據對SAR 數據的流場反演精度進行驗證。

采用的SAR 數據源于Radarsat-2 衛星，為標準成像模式下的單視復數據。結合星地匹配實驗中現場實測數據的位置及時間信息，本文選用了兩幅Radarsat-2 SAR 圖像進行流場反演，其UTC 時間分別為2019 年6 月23 日21 時53 分和2019 年6 月25 日10 時11 分。數據信息如表1 所示，SAR 數據覆蓋區域如圖1 所示。SAR 覆蓋區域內，海表面流即包括有一定規律性的地轉流和潮流，也包括無規律的風生流，所以整體來說該區域的海表面流場分布無明顯規律。

圖1 所用Radarsat-2 SAR 數據覆蓋區域Fig.1 Coverage area of Radarsat-2 SAR data

表1 本文所用Radarsat-2 SAR 數據信息Table 1 Radarsat-2 SAR data information used in this paper

風場數據采用歐洲中期天氣預報中心（European Centre for Medium-Range Weather Forecasts，ECMWF）ERA5 中的再分析風場數據。數據中包括海面上空10 m處風場的東向分量、北向分量和經緯度、時間信息，其空間分辨率為0.25°×0.25°，時間分辨率為1 h。在使用該風場數據時發現其時空分辨率比反演得到的SAR流速產品分辨率小，無法直接使用，所以本文通過時空插值使其與SAR 流速產品匹配，以方便之后的處理。

3 模型與方法

3.1 基于多普勒質心頻率異常的多普勒速度估算方法

從SAR 實測數據估算的多普勒中心頻率稱作實測多普勒中心頻率，主要包含衛星和地球表面的相對運動導致的多普勒頻移與海面水質點運動導致的多普勒頻移。實測多普勒中心頻率的估計方法有基于幅度的頻譜擬合法和基于相位的相位增量法[26–27]。本文利用頻譜擬合法從SAR 復數據中估計實測多普勒中心頻率。首先選取一定大小的窗口，對該窗口內的數據進行快速傅里葉變換（FFT）。然后在得到的方位向能量譜中找到多普勒中心頻率。由于直接從方位向能量譜中尋找多普勒中心頻率較為困難，本文結合能量均衡法，將方位向功率譜與參考函數相關，尋找相關后的過0 點，并把過0 點處的頻率作為該窗口的多普勒中心頻率[28]。最后以一定的步長移動窗口得到整幅SAR 圖像的多普勒中心頻率。

衛星和地球表面的相對運動產生的多普勒中心頻率稱作預測多普勒中心頻率，其計算方法主要有兩種：第一種方法基于衛星姿態和速度等參數，利用幾何關系估算該多普勒中心頻率；第二種方法從元數據中讀取多普勒系數、斜距時間、標準斜距時間等數據估算預測多普勒中心頻率。由于Envisat ASAR、Radarsat-2 等衛星，SAR 原始數據處理時已經估算了預測多普勒中心頻率，并給出了擬合多項式，所以本文基于第二種方法估算預測多普勒中心頻率，表達式為[29]

式 中，dopcoef1、dopcoef2、dopcoef3、dopcoef4、dopcoef5為多普勒系數；t為斜距時間；t0為標準斜距時間。

將衛星和地球表面的相對運動產生的預測多普勒中心頻率從實測多普勒中心頻率中去除之后得到的剩余頻率稱作多普勒質心頻率異常[30]

式中，fDc為實測多普勒中心頻率。經過式（2）的運算，衛星相對地球表面的運動效應被去除，因此多普勒質心頻率異常對應的只是海面水質點的運動速度。將多普勒質心頻率異常轉化為徑向多普勒速度，再將其投影到地表局部切平面坐標系得到地距多普勒速度

式中，k為雷達電磁波波數；θ為雷達入射角。

上述計算得到的地距多普勒速度對應海面水質點運動，該運動主要包括海表面流場、海浪軌道速度等，其中海浪受到海表面風場的調制作用，即風場可間接影響多普勒速度。Chapron 等[23]提出的多普勒模型通過波浪譜表示了風對多普勒速度的影響，并給出了計算風場和流場產生的多普勒速度的經驗表達式。

3.2 海面流場反演的弦截下山法

3.1 節獲得的地距多普勒速度中既包含了流場的貢獻又包含了風場的貢獻，只有準確去除風場的貢獻，才能獲得海面流場。利用M4S 模型仿真的多普勒譜信息可以計算多普勒質心偏移，進而獲得M4S模擬的多普勒速度場Vdop_m4s，即

當輸入M4S 模型的風場UW和流場UC信 息與真實SAR 成像時海面的風場和流場信息一致時，在忽略M4S 模型建模誤差的條件下，M4S 模擬的多普勒速度Vdop_m4s應與從真實SAR 圖像反演的多普勒速度Vdop_sar相等，即Vdop_m4s=Vdop_sar。在基于多普勒質心頻移法獲得了真實SAR 圖像反演的多普勒速度Vdop_sar的條件下，通過M4S 模型獲得海面真實流場，也就是求解以下關于海面流場UC的方程（假設海面真實風場UW已知）

上述方程是關于海面流場UC的 非線性方程，通過求解此方程可以獲得與SAR 圖像對應的海面流場。

這里采用迭代方法求解式（5），為了兼顧計算速度與迭代收斂性，選擇弦截法并結合下山法進行求解。弦截法不需要求解函數f(UC)的導數，容易實現；而下山法可以在一定程度上保證迭代的收斂性。弦截法求解非線性方程f(UC)=0的迭代公式為[31]

式中，Uk?1、Uk和Uk+1分 別表示第k?1、k和k+1次迭代的流速近似解。運算中需要首先把兩個初猜值U0和U1代入式（6），開啟迭代運算，直到相鄰兩次迭代的結果滿足 |Uk+1?Uk|

Uk+1、Uk為了滿足上述下山條件，在式（6）中引入下山因子A，其形式變為

下山因子A的取值為，···，直到使得下山條件式（7）成立為止。

基于M4S 模型，采用上述弦截下山法迭代反演流場的流程如圖2 所示，具體計算按如下步驟。

圖2 弦截下山法迭代反演流場流程Fig.2 Flow chart of iterative retrieval of current field by secant downhill method

（1）基于SAR 反演的海面多普勒速度Vdop_sar設置迭代初猜值，弦截法需要兩個初猜值，這里將Vdop_sar作為U0，U0+V′(V′=±0.1)作 為U1，并將其輸入M4S 模型，獲得V0dop_m4s、V1dop_m4s，在上述4 個值以及風場、初始化參數k、A的基礎上，開啟迭代運算。

（2）將由式（5）獲得的f(Uk)和f(Uk?1)的絕對值進行比較，如果滿足下山條件|f(Uk)|<|f(Uk?1)|則進入第（3）步處理，反之，進入第（4）步處理（當k=1時，通過選擇不同的V′的 值，使得 |f(U1)|<|f(U0)|成立）。

（3）利用式（6）計算非線性方程（5）新的近似解Uk+1，判斷 |Uk+1?Uk|是否滿足誤差閾值條件，滿足則停止迭代，以Uk+1作為最終的流場；若不滿足誤差閾值條件，將Uk+1輸 入M4S 模型獲得Vk+1dop_m4s，更新迭代參數k=k+1，A=1，轉入第（5）步操作。

（4）下山因子減半A=A/2，采用式（8）基于Uk?2和Uk?1重新計算Uk，判斷 |Uk?Uk?1|是否滿足誤差閾值條件，滿足則停止迭代，以Uk作為最終的流場；若不滿足誤差閾值條件，將新的Uk輸 入M4S 模型獲得新的Vkdop_m4s，轉入第（5）步操作。

（5）基于更新的Uk?1、Uk及Vk?1dop_m4s、Vkdop_m4s，轉入第（2）步，進行下一步迭代運算。

3.3 風貢獻因子 γ的估算方法

由于M4S 模擬仿真過程計算量較大，3.2 節的弦截下山法一般只適用于SAR 圖像局部區域的流場反演。對于整幅SAR 圖像，可以通過分塊迭代計算海面流場，估算每塊區域的風貢獻因子 γ，進而估算整幅SAR 圖像的風貢獻因子 γ。本文將從整幅SAR 圖像上均勻的選取大小相等的幾個局部區域（本文中是5 個），利用3.2 節提到的方法分別迭代反演其流場。將迭代反演獲得的流場、相應SAR 圖像反演的多普勒速度Vdop_sar以及外部風場數據（ECMWF 風場數據）代入式（9），獲得每個局部區域的風貢獻因子 γ。

由于風貢獻因子 γ與雷達頻率、入射角以及雷達天線的偏航角等因素有關[23]，而在一幅SAR 圖像中這些參數基本保持不變，所以本文中將各局部區域的風貢獻因子 γ的平均值作為整幅SAR 圖像的風貢獻因子γ。

4 結果與精度驗證

4.1 反演結果

采用3.1 節描述的方法，對兩幅Radarsat-2 SAR 數據進行處理，結果如圖3 和圖4 所示。從圖3 中可以發現，相比于預測多普勒中心頻率，多普勒質心頻率異常要小的多，這是由于兩者產生的原因不同，預測多普勒中心頻率由衛星和地球相對運動產生，而多普勒質心頻率異常由海表層運動產生。

經上述處理獲得的地距多普勒速度（圖3d、圖4）中即包含了流場貢獻又包含了風場貢獻，文獻[23]表明風場的貢獻甚至更大一些，因此只有準確的去除風場的貢獻，才能獲得最終的海面流場。這里首先采用3.2 節提出的弦截下山法對從SAR 圖像中截取的小區域進行迭代計算，獲得其海面流場。例如，對2019 年6 月23 日SAR 圖像反演的地距多普勒速度分布圖（圖3d）中紅框1 所示局部區域（大小為方位向5 km×距離向6 km）進行迭代計算，計算3 次的誤差|Uk+1?Uk|依次是0.34 m/s、0.07 m/s、0.02 m/s，已經滿足誤差閾值條件（e<0.05 m/s），因此停止迭代，輸出最終流場。然后，將該流場和相應SAR 圖像反演的多普勒速度Vdop_sar及風場數據代入式（9）估算該局部區域的風貢獻因子 γ，得到圖3d 中紅框1 區域的風貢獻因子 γ為0.15。用同樣的方法估算其他局部區域的風貢獻因子 γ，并將所有的風貢獻因子 γ進行平均，得到整幅SAR圖像的風貢獻因子γ。2019 年6 月23 日SAR 圖像的平均風貢獻因子 γ為0.15，2019 年6 月25 日SAR圖像的平均風貢獻因子 γ為0.22。將兩幅SAR 圖像的平均風貢獻因子 γ分別代入式（9）計算風場產生的多普勒速度，并將其從地距多普勒速度Vdop_sar中去除，得到整幅SAR 圖像的海面徑向流速，如圖5 所示。對SAR 地距多普勒速度圖中局部區域進行迭代計算時，通常只需迭代兩三次就可以滿足誤差閾值條件，獲得該區域的徑向流速。該結果表明本文提出的弦截下山法具有良好的收斂性和較高的收斂速度。

圖3 2019 年6 月23 日SAR 圖像反演結果Fig.3 SAR image retrieval results on June 23,2019

圖4 2019 年6 月25 日SAR 數據反演的地距多普勒速度Vdop_carFig.4 Ground range Doppler velocity retrieved from SAR data on June 25,2019

4.2 精度驗證

為了驗證本文方法，把SAR 反演的海表面徑向流速與星地匹配實驗中安德拉海流計采集的流場實測數據比對。安德拉海流計是常用的海流測量儀器，能夠準確測量海水瞬時速度。實驗中將海流計懸掛于浮體之下，使其處于海表面以下約0.4 m 處，保證其測量的是近表面海水的速度。獲取海表流場需要對海流計所測時間序列的速度取平均，以去除海浪的速度貢獻。將星地匹配實驗中測得的實測數據與SAR 圖像進行時空匹配，得到兩個實測數據分別與兩幅SAR 數據對應，具體數據如表2 所示，表中的流速、流向是衛星過境時刻前后10 min 內的平均值。

由于SAR DCA 方法反演的流速是一維徑向流速，在將海流計測得的流場與其比較時，需要將海流計所測二維流場向SAR 徑向作投影。2019 年6 月23 日雷達視向角為279.93°，2019 年6 月25 日雷達視向角為79.63°。將與之匹配的海流計所測流速分別投影到SAR 徑向，得到實測徑向流速分別為0.23 m/s、–0.14 m/s（如表2 中最后一列所示）。在SAR 反演的徑向流速分布圖中尋找實測數據的匹配位置點，如圖5a 和圖5b 中黑色五角星所示，以該點為中心，1 km×1 km 范圍內對SAR 反演的徑向流速取平均，得到兩幅SAR 圖像DCA 方法反演的流速分別為0.19 m/s，–0.29 m/s。將DCA 法反演的徑向流速與實測的徑向流速進行比較，得到兩幅SAR 圖像在實測點處的反演誤差分別為0.04 m/s 和0.15 m/s。在比對實驗中，SAR反演結果的精度與產品分辨率之間往往存在矛盾。平均范圍越大流速精度越高而分辨率會越低，反之，則流速精度越低而分辨率越高。所以需在兩者之間進行權衡，本文在1 km×1 km 的范圍內取平均，可在相對較高的分辨率（1 km×1 km）下獲得較高的精度。

表2 海流計實測海流數據Table 2 Current data measured by ocean current meters

圖5 兩景SAR 數據對應的海表面徑向流速Fig.5 The radial velocity of the sea surface corresponding to the SAR data of the two scenes

5 總結

為了克服SAR 復圖像DCA 方法反演海面流場時風場貢獻去除困難的難題，本文提出了基于M4S 模型的弦截下山法迭代反演海面流場。首先采用傳統的DCA 方法反演海面多普勒速度；然后采用本文提出的弦截下山法，迭代反演SAR 地距多普勒速度分布圖中局部區域的海面流場，并估算其風貢獻因子γ；最后對不同局部區域的風貢獻因子 γ取平均，獲得整幅SAR 地距多普勒速度分布圖的風貢獻因子 γ。進而去除風場對多普勒速度的貢獻，獲得整幅SAR 圖像對應的海面徑向流速。將兩幅SAR 圖像的反演結果與實測近海表面流速進行比對，得到的反演海面徑向流速偏差分別為0.04 m/s 和0.15 m/s。研究結果表明，本文提出的弦截下山法不僅具有良好的收斂性和較高的收斂速度，而且對于本文使用的兩景SAR 數據，反演的海面徑向流速偏差在0.2 m/s 內。