順軌干涉SAR海洋表面流場迭代反演算法

于祥禎 種勁松 洪 文

①(中國科學院電子學研究所微波成像技術國家級重點實驗室 北京 100190)

②(中國科學院研究生院 北京 100049)

1 引言

順軌干涉 SAR(Along-Track Interferometric SAR, ATI-SAR)是 1987年由美國 Goldstein和Zebker[1]首先提出來的一種技術，其干涉相位與表面流場的徑向速度分量成正比，利用順軌干涉SAR相位圖像可以直接獲得大面積、高分辨率的海洋流場和重要速度變化場信息。

由于海洋環境的復雜性以及SAR特殊的Bragg散射成像機制使得利用順軌干涉 SAR獲得的Doppler速度實際是各種波流速度(如表面流場速度、凈Bragg波相速度、大尺度波軌道速度等)的矢量和。順軌干涉SAR海洋表面流場反演方法研究，主要是不同速度分量分離方法的研究。

1993年，Shemer等人[2]在Monterey Bay的順軌干涉 SAR飛行試驗中考慮了沿不同方向傳播的Bragg波的貢獻。他們的試驗區域內包含一個蓄水池，其水體基本不流動，可以用該區域的順軌干涉相位來估計Bragg波的貢獻，但是其場景設置不具有普遍性。另外，他們通過空間平均的方法來消除大尺度波軌道速度的貢獻，但是他們的研究并沒有考慮大尺度波波前和波后對雷達截面調制的不同所帶來的附加速度分量的影響。1996年，Graber等人[3]給出了兩種方法來消除海面波浪運動的影響：一種是根據先驗風場和波浪信息并結合微波散射模型來計算波浪運動的貢獻；另一種是利用離散點的真實流場數據對順軌干涉相位獲得的表面流場進行定標。兩種方法獲得的表面流場基本都與岸基HF雷達獲得的表面流場相一致。傳統的凈Bragg相速度計算方法需要精確的風向信息。2003年，Kim 等人[4]根據L波段和C波段沿不同方向傳播的Bragg波分量所占比重近似相等的原理，提出了利用L, C雙波段順軌干涉SAR數據進行凈Bragg波相速度消除的方法。與傳統方法相比，該方法不需要風向信息。但能夠同時接收L和C兩個波段數據的順軌干涉SAR平臺較少，目前知道的僅有JPL的AIRSAR且于2004年已經退役。Romeiser等人[5,6]利用建立的順軌干涉SAR成像仿真模型，通過迭代校正的方法使仿真的干涉相位與實測干涉相位達到最佳匹配，并將此時的流場作為干涉測量的海洋表面流場的最優解，輸出流場的均方根誤差(Root-Mean-Square Error, RMSE)為0.1 m/s，但是文獻中并沒有對流場優化算法進行詳細介紹。

本文首先分析了傳統的速度分離方法的不足，指出在目前波浪運動的相關模型還不完善的情況下，要通過計算波浪運動的精確速度，然后再提取海洋表面流場速度比較困難；其次，根據文獻[5]開發的海面微波成像仿真模型M4S，給出了順軌干涉SAR海洋表面流場速度迭代反演的具體算法，并利用JPL AIRSAR獲得的機載順軌干涉SAR數據進行海洋表面流場反演，將反演結果與普林斯頓海洋模型(Princeton Ocean Model, POM)的輸出流場進行對比，驗證迭代反演算法的有效性。

2 順軌干涉SAR Doppler速度的組成分量

利用順軌干涉SAR直接獲得的Doppler速度實際是給定分辨單元內所有散射體徑向速度的矢量和。在中等入射角情況下，微波后向散射主要為Bragg散射，所以Doppler速度除了包含海表流場速度分量外，還有Bragg相速度、大尺度波軌道速度等的貢獻。Doppler速度的組成可以通過式(1)表示[7]。

其中uc為潮汐流或者洋流等表面流場，uwd為風場引起的表面漂移，uo為大尺度波浪的軌道速度，ub為凈Bragg波相速度。在實際流場反演中可以把前兩項統一看作表面流場，即us=uc+uwd。下面分別對波浪軌道速度和Bragg波相速度進行介紹。

2.1 大尺度波軌道速度

根據深水(水深大于波長的一半)重力波理論，海浪表面的質點在其平衡位置附近做圓周運動。通常情況下，軌道速度分量的去除通常是利用軌道速度沿海浪傳播方向具有周期性的特點，通過空間平均的方法來實現。但這種處理方式實際上存在如下幾個問題：

(1)要利用大尺度波軌道速度的周期性進行空間平均，首先需要知道精確的波浪波長信息，然而從目前的相關文獻看，很少有人首先對海浪波長進行估計，然后再進行平均處理；另外，也很難對整幅干涉相位圖像利用統一的標準進行空間平均。

(2)在實際海洋中，由于風等因素的影響，表面水質點并非進行一個閉合的軌道運動，而是存在一個與海浪傳播方向相同的凈水平速度，稱為 Stokes漂流。該項實際上也可以認為是表面流場us的一部分。

(3)根據 SAR海浪成像機制，朝向雷達一側的波浪的后向散射強度要高于背向雷達一側，從而使得兩側軌道速度分量的權重不一致，這樣就會使Doppler速度產生一個正的速度偏移。在中等海況條件下，該速度偏移可能會達到0.2 m/s[7]。該速度偏移并非真實表面流場的一部分，而是由SAR海浪成像機制導致的，但在傳統速度分離方法中，很少有文獻對此速度偏移進行處理。

由于上述原因，通過空間平均處理后的順軌干涉相位肯定會存在殘余軌道速度分量。

2.2 Bragg波相速度

Bragg波相速度可以通過式(2)計算：

其中g為重力加速度，τs為表面張力，ρ為海水密度，kb為Bragg波波數矢量。通常情況下，存在兩個傳播方向相反的Bragg波，即遠離雷達方向傳播的Bragg波和朝向雷達方向傳播的Bragg波。順軌干涉SAR Doppler速度中包含的凈Bragg相速度實際是這兩個方向Bragg波相速度的加權矢量和，即

其中α(θw)為遠離雷達方向傳播的 Bragg波所占的比重，其大小與風向θw有關。在順風情況下(θw=0o),α=1；在逆風情況下(θw=1 80o),α=0；在側風情況下，Bragg波所占比重與風向的關系還沒有一個精確的理論模型。

從上述分析可知，目前波浪運動的相關模型還不完善，要通過計算波浪運動的精確速度，然后再提取表面流場速度比較困難。下面我們基于文獻[5]開發的海面微波成像仿真模型(M4S)建立了順軌干涉SAR海洋表面流場的迭代反演算法。該算法通過流場的迭代校正，使得仿真得到的順軌干涉相位與實際干涉相位相一致，此時的流場即為“最優流場”。該算法不用對軌道速度或Bragg相速度等進行精確計算，避免了波浪運動模型本身問題引入的誤差。

3 順軌干涉 SAR海洋表面流場迭代反演算法

文獻[5]建立了一種順軌干涉 SAR成像仿真模型，并集成到其開發的微波海洋成像仿真模型?M4S中。其順軌干涉相位計算的主要原理是順軌干涉SAR的相位與后向散射場自相關函數的相位相一致，從而通過計算后向散射場自相關函數的相位來間接計算順軌干涉相位。文獻[8]結合SRTM X-SAR系統參數和 KUSTWAD流場進行順軌干涉相位圖像仿真，并與真實SRTM干涉圖像進行對比，結果發現兩者基本一致，從而驗證了利用M4S進行順軌干涉相位仿真的有效性。文獻[6]基于M4S模型提出了順軌干涉SAR流場迭代反演算法，但對其中的流場優化算法部分并沒有進行介紹。下面基于M4S模型建立順軌干涉SAR流場迭代反演算法，并給出流場校正的具體算法。

利用順軌干涉 SAR進行表面流場迭代反演時需要給定流場的初始值，即初猜流場。初猜流場可以通過順軌干涉相位圖像直接得到，即

其中u0為初猜流場，φ0為實際獲得的順軌干涉相位，λ為雷達波長，B為有效基線長度(單發雙收的雙站工作模式下，有效基線為物理基線長度的一半；單站模式下，有效基線與物理基線相等)，V為平臺速度，θ為入射角。

圖1給出了基于M4S模型的順軌干涉SAR表面流場迭代反演算法的流程圖，具體實現如下：

圖1 順軌干涉SAR海洋表面流場迭代反演算法流程圖

(1)首先需要對流場校正過程需要的部分參數進行初始化，如圖1所示。其中，n為迭代次數；T1為干涉相位均方根誤差閾值。T2為干涉相位偏差閾值；Fij為流場校正標志。

(2)將計算得到的表面流場以及風場、雷達參數、平臺參數等輸入M4S模型中計算順軌干涉相位圖像。

(3)方位向偏移校正。由于流場方位向梯度的存在，部分像素有可能校正到同一個像素單元，而部分像素單元產生“空洞”，這里我們對方位向偏移校正后的干涉相位圖像進行了3×3窗口的均值濾波。

(4)將仿真的干涉相位圖像與實際干涉相位圖像進行對比，如果干涉相位的均方根誤差rmsen＜T1，則停止迭代并認為此時的流場即為“最優流場”；否則，與上次迭代時的均方根誤差相比較，如果rmsen＞r msen-1，則認為迭代已經發散，停止迭代并輸出上一次迭代生成的流場；如果rmsen≤rmsen-1，則認為迭代仍在收斂，跳入步驟(5)。

下面我們利用 JPL AIRSAR獲得的機載順軌干涉SAR數據對該迭代算法進行驗證。

4 算法驗證

4.1 數據說明

下面利用JPL AIRSAR于1998年10月29日(GMT 22:23)獲得的美國Key Largo附近海域的L波段機載順軌干涉 SAR數據進行海洋表面流場提取。圖2給出了數據獲取的具體位置。在成像區域內有一個海上觀測站(MLF1，北緯 25°，西經80.39°)，該觀測站可以記錄波浪參數以及風場參數，但是并沒有流場的相關信息。根據MLF1的歷史記錄，AIRSAR成像時刻觀測站附近海面風速為 7 m/s，風向為北偏東73°。在成像區域內還有部分陸地，在數據處理時可以用來對順軌干涉相位進行定標(理想情況下，陸地的順軌干涉相位應該為零)。表1給出了雷達參數、平臺參數以及風場參數。

4.2 順軌干涉SAR數據處理

圖2 成像區域示意圖

表1 雷達參數、平臺參數及風場信息

圖3給出了地距的順軌干涉相位圖像，其中方位向×地距向像素間隔為5 m×5 m。為了減小噪聲干擾，對干涉相位圖像進行了2次7×7窗口的圓周期均值濾波，濾波后的干涉相位如圖4所示。從圖4可以看出，由于通道不平衡以及平臺俯仰或偏航等因素導致陸地區域(圖像下部)的干涉相位不為零。因為理想情況下，陸地區域的順軌干涉相位為零，所以我們可以利用陸地區域的干涉相位對整幅圖像進行干涉相位定標，定標后結果如圖5所示。

由于 M4S對輸入流場的尺寸有限制(流場最大為1001×1001)，所以我們對獲得的干涉相位圖像進行多視處理，多視數為 20×20，處理后方位向和距離向的像素間隔為100 m，采樣點數分別為258和112(雖然多視處理降低了圖像的空間分辨率，但是由于提高了圖像的信噪比和干涉相位精度，從而有助于提高流場的測速精度)。圖6給出了多視后的干涉相位圖像。根據表1給定的參數可知，表面流場經過速度聚束調制引起的方位向偏移如果要大于100 m(即超出一個像素單元)，其徑向速度分量至少要達到2.2 m/s，而在現實海洋中超過2 m/s的表面流場比較少，所以在下面流場迭代反演時不需要再進行方位向偏移校正。根據式(4)可以獲得初猜流場，結果如圖7所示。

4.3 表面流場反演結果

將圖7的初猜流場代入流場迭代反演模型中，經過5次迭代達到收斂。圖8給出了方位向剖面a和距離向剖面r處(具體位置見圖7)的流場迭代校正情況。圖9給出了收斂時仿真的干涉相位與真實干涉相位的統計對比情況，從圖中可以看出，仿真的干涉相位與真實相位基本一致，兩者的相關系數約為 0.998，均方根誤差 RMSE為 0.027 rad。圖10給出了最終的輸出流場。

因為試驗時沒有進行同步的海洋流場觀測，且觀測站MLF1無法記錄流場信息，附近海面也沒有浮標數據可供參考，所以無法通過實測數據對本文迭代反演流場進行驗證。下面我們通過海洋動力仿真模型 POM 的輸出結果對我們的迭代反演結果進行交叉驗證。實際上文獻[9]在對其流場反演算法進行驗證時也通過與環流數值仿真模型進行對比，他們使用的數值仿真模型 KUSTWAD主要針對的是荷蘭Wadden海域。

4.4 POM海洋表面流場仿真

POM 是由美國普林斯頓大學大氣海洋科學項目組研制的3維海洋動力數值仿真模型。目前POM已經廣泛應用于近海和河口、灣流以及區域海等的流場分布、溫鹽分布、潮汐變化、波浪破碎、污染監控等[10-13]的模擬并取得很大成功。下面我們利用POM來仿真Key Largo附近海域的海洋流場分布。

POM模型采用水位邊界條件，并設置場景初始速度場ux=uy=uz=0 ，其中ux,uy和uz分別為水平方向和垂直方向的速度分量。初始水深選用全球水深數據ETOP05，溫鹽參數采用美國NOAA提供的月平均溫度和鹽度數據，初始風場采用觀測站MLF1提供的風場數據，并假定整個海面風場分布是均勻的。圖11給出表面流場穩定后的輸出結果。為了與 POM 仿真結果相對比，這里我們將迭代反演流場進行地理映射，結果如圖12所示。

將圖12標示出的7個位置點處的徑向流場結果進行對比，對比結果如圖13所示，從圖中可以看出迭代反演流場與 POM 仿真流場的 RMSE為 0.13 m/s，相關系數為0.854，可見兩者基本一致。如果不考慮引起較大誤差的第⑦參考點，則兩者的RMSE為0.06 m/s，相關系數為0.98。

圖3 地距的順軌干涉相位

圖4 濾波后的干涉相位

圖5 定標后的干涉相位

圖6 多視處理后的干涉相位

圖7 初猜流場

圖8 流場迭代校正情況

圖9 收斂時的相位與真實相位對比

圖10 最終輸出的表面流場

圖11 POM仿真流場(表層流場)

圖12 迭代反演流場地理映射后的結果

圖13 迭代反演流場與POM仿真流場的統計對比

5 總結

本文首先對順軌干涉相位直接獲得的Doppler速度組成進行介紹，并分析了傳統速度分離方法的不足，指出在目前波浪速度模型尚不完善的情況下，基于海面成像仿真模型的流場迭代反演算法是進行順軌干涉SAR海洋表面流場反演的一種比較可行的方法。并根據Romersier等提出的海面微波成像仿真模型M4S，給出了順軌干涉SAR海洋表面流場反演的具體迭代校正方法。利用該迭代算法對JPL AIRSAR獲得的美國Key Largo海域的機載順軌干涉SAR數據進行了海洋表面流場反演，并將反演結果與海洋動力仿真模型POM的輸出結果進行了對比，兩者具有很好的一致性，證明了該迭代反演算法的有效性。

[1]Goldstein R M and Zebker H A. Interferometric radar measurement of ocean surface currents[J].Nature, 1987, 328:707-709.

[2]Shemer L, Marom M, and Markman D. Estimates of currents in the nearshore ocean region using interferometric synthetic aperture radar[J].Journal of Geophysical Research, 1993,98(C4): 7001-7010.

[3]Graber H C, Thompson D R, and Carande R E. Ocean surface features and currents measured with synthetic aperature radar interferometry and HF radar[J].Journal of Geophysical Research, 1996, 101(C11): 25813-25832.

[4]Kim D, Moon W M, Moller D,et al.. Measurements of ocean surface waves and currents using L-and C-band along-track interferometric SAR[J].IEEE Transactions on Geoscience and Remote Sensing, 2003, 41(12): 2821-2832.

[5]Romeiser R and Thompson D R. Numerical study on the along-track interferometric radar imaging mechanism of oceanic surface currents[J].IEEE Transactions on Geoscience and Remote Sensing, 2000, 38(1): 446-485.

[6]Romeiser R. Current measurements by airborne along-track InSAR: measuring technique and experimental results [J].IEEE Journal of Oceanic Engineering, 2005, 30(3): 552-569.

[7]Moller D and Frasier S J. Radar-derived interferometric surface currents and their relationship to subsurface current structure[J].Journal of Geophysical Research, 1998, 103(C6):12839-12852.

[8]Romeiser R. On the potential of current measurements by spaceborne along-track InSAR for river runoff monitoring[C].Oceans-2005 Europe, Brest, France, 2005: 787-792.

[9]Romeiser R, Breit H, Eineder M,et al.. Current measurements by SAR along-track interferometry from a space shuttle[J].IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(10): 2315-2324.

[10]Chang Y, Li X, Xu Q,et al.. SAR observation and model tracking of an oil spill event in coastal waters[J].Marine Pollution Bulletin, 2011, 62(2): 350-363.

[11]Huang C J, Qiao F, Song Z,et al.. Improving simulations of the upper ocean by inclusion of surface waves in the Mellor-Yamada turbulence scheme[J].Journal of Geophysical Research, 2011, 116, C1007, doi:10.1029/2010JC006320.

[12]Zhang X, Han G, Wang D,et al.. Effect of surface wave breaking on the surface boundary layor of temperture in the Yellow sea in the summer[J].Ocean Modelling, 2011, 38(3-4):267-279.

[13]Santoki M, Ratheesh S, Sharma R,et al.. Assimilation of drifter data in a circulation model of the Indian ocean[J].IEEE Geoscience and Remote Sensing Letters, 2012, 9(1):100-103.