基于SAR 多普勒質(zhì)心頻移的海面流場(chǎng)迭代反演算法

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

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

1 引言

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

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

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

2 數(shù)據(jù)源

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

采用的SAR 數(shù)據(jù)源于Radarsat-2 衛(wèi)星，為標(biāo)準(zhǔn)成像模式下的單視復(fù)數(shù)據(jù)。結(jié)合星地匹配實(shí)驗(yàn)中現(xiàn)場(chǎng)實(shí)測(cè)數(shù)據(jù)的位置及時(shí)間信息，本文選用了兩幅Radarsat-2 SAR 圖像進(jìn)行流場(chǎng)反演，其UTC 時(shí)間分別為2019 年6 月23 日21 時(shí)53 分和2019 年6 月25 日10 時(shí)11 分。數(shù)據(jù)信息如表1 所示，SAR 數(shù)據(jù)覆蓋區(qū)域如圖1 所示。SAR 覆蓋區(qū)域內(nèi)，海表面流即包括有一定規(guī)律性的地轉(zhuǎn)流和潮流，也包括無(wú)規(guī)律的風(fēng)生流，所以整體來(lái)說(shuō)該區(qū)域的海表面流場(chǎng)分布無(wú)明顯規(guī)律。

圖1 所用Radarsat-2 SAR 數(shù)據(jù)覆蓋區(qū)域Fig.1 Coverage area of Radarsat-2 SAR data

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

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

3 模型與方法

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

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

衛(wèi)星和地球表面的相對(duì)運(yùn)動(dòng)產(chǎn)生的多普勒中心頻率稱作預(yù)測(cè)多普勒中心頻率，其計(jì)算方法主要有兩種：第一種方法基于衛(wèi)星姿態(tài)和速度等參數(shù)，利用幾何關(guān)系估算該多普勒中心頻率；第二種方法從元數(shù)據(jù)中讀取多普勒系數(shù)、斜距時(shí)間、標(biāo)準(zhǔn)斜距時(shí)間等數(shù)據(jù)估算預(yù)測(cè)多普勒中心頻率。由于Envisat ASAR、Radarsat-2 等衛(wèi)星，SAR 原始數(shù)據(jù)處理時(shí)已經(jīng)估算了預(yù)測(cè)多普勒中心頻率，并給出了擬合多項(xiàng)式，所以本文基于第二種方法估算預(yù)測(cè)多普勒中心頻率，表達(dá)式為[29]

式 中，dopcoef1、dopcoef2、dopcoef3、dopcoef4、dopcoef5為多普勒系數(shù)；t為斜距時(shí)間；t0為標(biāo)準(zhǔn)斜距時(shí)間。

將衛(wèi)星和地球表面的相對(duì)運(yùn)動(dòng)產(chǎn)生的預(yù)測(cè)多普勒中心頻率從實(shí)測(cè)多普勒中心頻率中去除之后得到的剩余頻率稱作多普勒質(zhì)心頻率異常[30]

式中，fDc為實(shí)測(cè)多普勒中心頻率。經(jīng)過(guò)式（2）的運(yùn)算，衛(wèi)星相對(duì)地球表面的運(yùn)動(dòng)效應(yīng)被去除，因此多普勒質(zhì)心頻率異常對(duì)應(yīng)的只是海面水質(zhì)點(diǎn)的運(yùn)動(dòng)速度。將多普勒質(zhì)心頻率異常轉(zhuǎn)化為徑向多普勒速度，再將其投影到地表局部切平面坐標(biāo)系得到地距多普勒速度

式中，k為雷達(dá)電磁波波數(shù)；θ為雷達(dá)入射角。

上述計(jì)算得到的地距多普勒速度對(duì)應(yīng)海面水質(zhì)點(diǎn)運(yùn)動(dòng)，該運(yùn)動(dòng)主要包括海表面流場(chǎng)、海浪軌道速度等，其中海浪受到海表面風(fēng)場(chǎng)的調(diào)制作用，即風(fēng)場(chǎng)可間接影響多普勒速度。Chapron 等[23]提出的多普勒模型通過(guò)波浪譜表示了風(fēng)對(duì)多普勒速度的影響，并給出了計(jì)算風(fēng)場(chǎng)和流場(chǎng)產(chǎn)生的多普勒速度的經(jīng)驗(yàn)表達(dá)式。

3.2 海面流場(chǎng)反演的弦截下山法

3.1 節(jié)獲得的地距多普勒速度中既包含了流場(chǎng)的貢獻(xiàn)又包含了風(fēng)場(chǎng)的貢獻(xiàn)，只有準(zhǔn)確去除風(fēng)場(chǎng)的貢獻(xiàn)，才能獲得海面流場(chǎng)。利用M4S 模型仿真的多普勒譜信息可以計(jì)算多普勒質(zhì)心偏移，進(jìn)而獲得M4S模擬的多普勒速度場(chǎng)Vdop_m4s，即

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

上述方程是關(guān)于海面流場(chǎng)UC的 非線性方程，通過(guò)求解此方程可以獲得與SAR 圖像對(duì)應(yīng)的海面流場(chǎng)。

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

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

Uk+1、Uk為了滿足上述下山條件，在式（6）中引入下山因子A，其形式變?yōu)?/p>

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

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

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

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

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

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

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

（5）基于更新的Uk?1、Uk及Vk?1dop_m4s、Vkdop_m4s，轉(zhuǎn)入第（2）步，進(jìn)行下一步迭代運(yùn)算。

3.3 風(fēng)貢獻(xiàn)因子 γ的估算方法

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

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

4 結(jié)果與精度驗(yàn)證

4.1 反演結(jié)果

采用3.1 節(jié)描述的方法，對(duì)兩幅Radarsat-2 SAR 數(shù)據(jù)進(jìn)行處理，結(jié)果如圖3 和圖4 所示。從圖3 中可以發(fā)現(xiàn)，相比于預(yù)測(cè)多普勒中心頻率，多普勒質(zhì)心頻率異常要小的多，這是由于兩者產(chǎn)生的原因不同，預(yù)測(cè)多普勒中心頻率由衛(wèi)星和地球相對(duì)運(yùn)動(dòng)產(chǎn)生，而多普勒質(zhì)心頻率異常由海表層運(yùn)動(dòng)產(chǎn)生。

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

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

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

4.2 精度驗(yàn)證

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

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

表2 海流計(jì)實(shí)測(cè)海流數(shù)據(jù)Table 2 Current data measured by ocean current meters

圖5 兩景SAR 數(shù)據(jù)對(duì)應(yīng)的海表面徑向流速Fig.5 The radial velocity of the sea surface corresponding to the SAR data of the two scenes

5 總結(jié)

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