基于VMD對SAR海洋內波參數的自動反演

張洪生，鄭應剛，王有強

(上海海事大學 海洋科學與工程學院，上海 201306)

海洋內波是在海水密度穩(wěn)定分層的海洋中產生的、最大振幅發(fā)生在海洋深處的波動現象。它通常是由海洋內部擾動源而引起的[1]，對海洋的質量、動量和能量的傳輸有重要作用[2]。對包括內孤立波在內的內波進行研究越來越引起人們的關注。直接觀測是研究內孤立波的最佳手段[3]，可采用現場測量海水參數(溫鹽密、海水流速等)的方法對內波進行直接精確測量。但是由于直接觀測法造價昂貴，儀器容易損害和丟失等，除了對一些重要的海域需要進行觀測之外，開展大范圍長時間的觀測是不現實的。

目前，利用合成孔徑雷達(synthetic aperture radar，簡稱SAR)圖像進行內波研究是一個亟待發(fā)展的方向[4]。內波在傳播過程中引起了海洋表面的輻聚輻散，使得海洋表面的粗糙度發(fā)生改變，因此內波在SAR圖像上表現為明暗相間的條帶[5]。SAR具有全天時、全天候、距離遠、范圍大、分辨率高等優(yōu)點，已成為觀測內波不可或缺的手段[6]。利用SAR圖像研究內波一方面可以從SAR圖像中反演內波的波長、波包間距等參數，另一方面可以利用實測資料進而反演內波的振幅和混合層深度等參數[7]。

基于SAR圖像進行內波參數反演的方法主要有傅里葉變換法、小波變換法和經驗模態(tài)分解(empirical mode decomposition，簡稱EMD)法等。傅里葉變換容易獲取內波的主波信息[8]，而小波變換能夠獲取內波波形信息，但難在根據具體問題選擇出最合適的小波基函數[5]。EMD方法是由Huang等[9]在1998年提出的，主要應用于物理學、生物學和海洋遙感等領域。EMD將信號序列分解成許多個本征模態(tài)函數(intrinsic mode function，簡稱IMF)分量，這些分量具有不同的特征尺度，與分解前的原始信號數據序列相比更有規(guī)律，但其存在模態(tài)混疊問題，影響到后續(xù)IMF的分離[10]，最終影響到內波參數的反演效果。

楊勁松等[10]分別利用傅里葉變換和小波變換提取出了SAR圖像中的內波波長和波向等參數。甘錫林等[11]通過比較小波變換、傅里葉變換和EMD三種方法所提取的內波波長參數，驗證了EMD方法的有效性。陳捷等[12]通過圓擬合算法估計SAR圖像中內波傳播方向，進而選取內波斷面并利用 EMD 反演內波波長和振幅。Wu等[13]通過添加一種噪聲輔助數據分析方法，提出了集合經驗模態(tài)分解(ensemble empirical mode decomposition，簡稱EEMD)法，有效地克服了EMD模態(tài)混疊的問題。汪雄良等[5]運用EEMD方法所反演的前導波振幅與實測資料非常接近，結果表明EEMD方法較好地克服了EMD所存在的模態(tài)混疊問題。Yeh等[14]針對EEMD加入的白噪聲不能完全中和的問題，提出了一種新的噪聲增強數據分析方法——互補集合經驗模態(tài)分解(complementary ensemble empirical mode decomposition，簡稱CEEMD)法，有效地消除了IMF中的殘留噪聲。該方法廣泛地應用于信號的去噪[15]和故障特征的提取[16]等。

變分模態(tài)分解(variational mode decomposition，簡稱VMD)是基于嚴謹的數學理論而提出的一種新的信號分解方法，具有良好的魯棒性和信號分解精度[17]，目前廣泛應用于故障診斷等方面[18]。VMD分解后的固有模態(tài)函數的個數K需預先設置，通過不斷地調整參數K就能解決模態(tài)混疊的問題。

目前，在對SAR影像中的海洋內波參數進行反演研究時，大多是人為選取內波斷面數據或者通過輔助點進行擬合選取進而分析斷面數據。文中通過獲取SAR圖像在Canny預處理后的條紋信息，再根據內波傳播方向自動選取灰度剖面，然后利用EEMD信號特征自適應分解模態(tài)函數的特點，再將分解得到的有效模態(tài)數作為VMD中參數K的參考值，最后利用VMD分解后的數據進行內波參數反演。該方法解決了人為選取剖面可能導致的誤差。通過對比分析EMD方法、EMMD方法和VMD方法的試驗結果，并將VMD方法和EMMD方法反演的結果與實測數據資料進行比較，論證了VMD方法的可行性。

1 方法原理

1.1 基于KdV方程內波參數的反演原理

在內波參數反演中，最常用的是內孤立波理論。內孤立波理論是以KdV方程為基礎的淺水波理論。一維非線性內孤立波方程為：

(1)

式中：η為內波波面位移，t為時間，c0為線性速度，α為非線性參數，β為頻散參數。

假定海洋由兩層密度不同的海水組成，即可概化為兩層模型[19]。假設上層流體的深度為h1、密度為ρ1，下層流體的深度為h2、密度為ρ2，總水深H為：

H=h1+h2

(2)

求解式(1)可求得穩(wěn)態(tài)的內孤立波解：

(3)

式中：η0為內孤立波的最大振幅，C為內波相速度，l為內孤立波半振幅寬度。C，c0和l的表達式分別為：

(4)

(5)

(6)

在海洋內波SAR圖像中，由于內波的存在改變了海面的粗糙度，由此引起SAR圖像灰度值的相對變化[20]。這種相對變化可以定義為：

(7)

(8)

式中：I0為海面的背景強度值，ΔI為含有內波區(qū)域的強度值與海面的背景強度值之差。x′=x-Ct為與內波一起運動的坐標，μ為松弛率，γ為Bragg波群速度與相速度之比，φ為內波傳播方向(x方向)與雷達視向的夾角，B是與內波相速度、振幅、松弛率、內波深度和雷達視向等相關的系數。

由式(7)計算內波明暗條紋關系：

(9)

可得式(9)的解為：

x′=±0.66l

(10)

單個內波最亮點和最暗點之間的距離D可表示為：

D=1.32l

(11)

由式(6)可知：

(12)

由式(11)可知要反演出半振幅寬度l需要首先知道最亮點和最暗點之間距離D，由式(12)可知要反演出內波振幅除了需要已知半振幅寬度l外，還需要知道上下層流體深度h1和h2。

1.2 變分模態(tài)分解的原理

VMD實際上是對構造的變分模型進行尋優(yōu)的一種算法，對所分解模態(tài)函數的中心頻率和帶寬采用交替方向乘子法(altarmate direction method of multipliers，簡稱ADMM)更新。首先保持中心頻率不變更新帶寬，然后再保持帶寬不變更新中心頻率，最后設置迭代結束的閾值。當達到預定的閾值則迭代結束并得到分解的信號，根據實際信號的特征頻帶自適應分解成多個IMF分量[17]。

VMD使用了EMD中IMF的定義，根據調制準則重新將IMF定義為調幅調頻(amplitude modulated frequency modulated，簡稱AM-FM)的信號：

uk(t)=Ak(t)cosφk(t)

(13)

1) 首先對模態(tài)函數uk(t)進行希爾伯特變換(Hilbert transform)以求得解析信號：

(14)

2)得到的解析信號乘以e-jωkt，以使各模態(tài)函數移動到中心頻率ωk：

(15)

3) 利用H1高斯平滑求解各模態(tài)函數的帶寬之和，即L2范數梯度的平方：

(16)

4) 構造帶約束的變分模型，在該約束下求各個模態(tài)函數的最優(yōu)解：

(17)

(18)

5) 為了求解上述約束問題的解，引入二次懲罰項和拉格朗日(Lagrange)乘子將約束性問題轉化為非約束性問題并求解，構造如下形式的增廣拉格朗日函數：

(19)

式中：α為懲罰因子，取默認值2 000，λ為拉格朗日乘子，L為增廣拉格朗日算子，〈〉為內積。

使用交替方向乘子法處理式(19)，得到模型的最優(yōu)解，從而將原信號分解成K個IMF信號。具體步驟如下：

2)n=n+1；

VMD按照信號特征進行分離，形成多個信號分量。其中，參數K的選取非常關鍵，并直接影響到參數反演的準確性。利用EEMD方法不需要事先設置分解模態(tài)個數且能夠自適應分解的特點，將EEMD分解的層數作為VMD中參數K的取值。

2 圖像預處理和數據分析

2.1 圖像的預處理

研究所使用的SAR圖像，為在巴士海峽附近于1995年6月16日02時29分發(fā)生海洋內波時的圖像，選取的圖像是歐空局ERS-1星載衛(wèi)星所拍攝的SAR圖像，圖像像素大小為9 487×6 130，成像波段為C波段，極化方式為VV，像素之間的距離為12.5 m，經緯度范圍為19°46′ N～23°03′ N，119°51′ E～121°01′ E。

ERS-1原始SAR圖像如圖1所示。在提取內波信號之前，需要對SAR圖像進行預處理。由于下載的產品是單視一級產品，這是一個傾斜的范圍圖像，所包含復雜的數據尚未進行多視處理。首先需要對圖像進行輻射定標(calibration radiometrically)，使像素坐標真實地反應雷達反射面的后向散射強度；其次將圖像進行多視處理(multilook processing)，將圖像從傾斜范圍轉化為地面范圍后，圖像具有較小的噪聲和近似的平方像素間距；最后采用Gamma map[21]濾波對圖像斑點噪聲抑制，以減少對圖像分辨率和信號分離的影響。采用歐空局(ESA)提供的SNAP軟件對以上圖像進行預處理，然后通過Canny處理得到篩選后的內波條紋信息[22]，最終得到的圖像如圖2所示。

圖1 ERS-1原始SAR圖像Fig. 1 ERS-1 original SAR image

圖2 海洋內波SAR圖像預處理結果Fig. 2 Image preprocessing results of oceanic internal wave with SAR

基于最小二乘法對圖2中每一條內波條紋所構成的點集進行直線擬合K，得到每一條內波條紋的大致方向N，然后將N作為初始值，根據兩點間距離最小性質進行迭代計算出每一條紋上每一點的真實方向N*，最后根據真實方向自動選取通過該點的灰度剖面，剖面示意見圖3。此外，也可以利用該點所在的局部曲線段進行垂線求解得到該點剖面的大致方向。

圖3 迭代法選擇剖面示意Fig. 3 Sketch of selection profile by iterative method

為了便于驗證反演算法的有效性，選擇了與文獻[5]進行研究所采用的大致相同的區(qū)域。結合該地區(qū)已有的實測數據，對海洋內波振幅進行反演。因此僅對黑色線框區(qū)域進行研究，如圖4所示。為了提取內波信號，選擇了A點處的灰度剖面數據作為研究對象。

圖4 灰度剖面數據位置Fig. 4 Location of gray profile data

2.2 數據處理結果與分析

以A點處的灰度剖面數據作為原始信號，分別做EMD、EEMD和VMD分解。其中EMD和EEMD均自適應分解為8層模態(tài)分量。分解結果如圖5和6所示。

圖5 內波剖面數據的經驗模態(tài)(EMD)分解Fig. 5 Empirical mode decomposition of internal wave profile data

圖6 內波剖面數據的集合經驗模態(tài)(EEMD)分解Fig. 6 Ensemble empirical mode decomposition of internal wave profile data

楊洪柏等[23]針對EMD不需預先設定模態(tài)數自適應的分解特點，通過對EMD分解結果的分析，對VMD分解模態(tài)數K進行估計。針對EEMD同樣具有自適應分解為一定數量模態(tài)數的特點，對分解的模態(tài)函數進行快速傅里葉變換(fast fourier transformation，簡稱FFT)頻譜分析，結果如圖7所示，進而近似求得VMD的分解模態(tài)數K。由圖7可見，模態(tài)1頻譜頻帶很寬，可視為噪聲，無頻率中心。模態(tài)4和模態(tài)5的中心頻率兩者都在1.672 Hz左右，可認為這兩個分量具有一個頻率中心，即可視為一個模態(tài)數。模態(tài)2、3、6、7和8可以分別作為一個獨立的模態(tài)。因此VMD的分解模態(tài)數K可取6。VMD分解的結果如圖8所示。

圖7 集合經驗模態(tài)分解各模態(tài)分量的頻譜Fig. 7 Ensemble empirical mode decomposition of the spectrum of each modal component

圖8 內波剖面數據的變分模態(tài)(VMD)分解Fig. 8 Variational mode decomposition of internal wave profile data

對原始信號進行分解的結果必有一層所含的內波信息較豐富?？梢岳糜嬎惴纸獾玫降拿繉有盘枤w一化方差，來表示其內波相對能量的大小。因為內波蘊含著巨大的能量，可以用數值最大的歸一化方差分量代表內波信息[10]。歸一化方差表達為：

(20)

式中：σi為各層模態(tài)分量的方差，n為分解的層數。

分別計算每一層的歸一化方差，計算結果如表1所示。由表1可知利用EMD方法分解之后的模態(tài)函數IMF5(EMD_IMF5)最大歸一化方差為0.710 0。利用EEMD方法分解之后的模態(tài)函數IMF5(EEMD_IMF5)最大歸一化方差為0.596 9。利用VMD方法分解之后的模態(tài)函數U2(VMD_U2)最大歸一化方差為0.938 1。歸一化方差的最大值對應的模態(tài)函數所含的內波能量也最大，因此可代表內波信號。

表1 經驗模態(tài)分解、集合經驗模態(tài)分解和變分模態(tài)分解各模態(tài)分量的歸一化方差

將原始信號、EMD_IMF5、EEMD_IMF5和VMD_U2繪制在一起，如圖9所示。從圖9可看出，EEMD_IMF5和VMD_U2能夠真實反映出原始信號上下起伏的周期性特征，同時驗證了歸一化方差最大的分量可代表內波信息這一準則。從圖9也能看出，EMD_IMF5在區(qū)間[300，430]出現了模態(tài)混疊的現象，其中區(qū)間[350，400]中的內波信號出現在EMD分解的第四個模態(tài)函數IMF4中，模態(tài)混疊問題的出現影響到了下一步的分解。這就不能真實地反映出原始信號的特征。然而EEMD_IMF5和VMD_U2信號特征相近、其周期性強、上下起伏穩(wěn)定，能夠反映出真實的內波物理特性，因此這兩者適合做內波參數的反演。

圖9 原始信號、EMD_IMF5、EEMD_IMF5和VMD_U2Fig. 9 Original signal， EMD_IMF5， EEMD_IMF5 and VMD_U2

選取EEMD_IMF5和VMD_U2的數據。從左端開始算起，得到最亮點和最暗點之間的橫坐標，進而計算最亮點和最暗點之間的像素距離差；然后根據相鄰像素之間的距離為12.5 m，求得最亮點和最暗點的實際間距D；再根據式(11)可計算得到內波的半振幅寬度l。計算結果如表2和3所示。

表2 EEMD_IMF5中各孤立子波亮暗條紋之間的距離Di和半振幅寬度liTab. 2 The distance between the light-dark strip and the half-amplitude width of each solitary wave in EEMD_IMF5

表3 VMD_U2中各孤立子波亮暗條紋之間的距離Di和半振幅寬度liTab. 3 The distance between the light-dark strip and the half-amplitude width of each solitary wave in VMD_5

借助研究區(qū)域的實測數據，所研究的遙感圖像區(qū)域位置混合層深度大約為46 m，總水深大約為4 000 m，前導波振幅大約為100 m左右[5， 11]。分別利用表2和3中的序號1半振幅寬度l，根據式(12)反演出該地區(qū)前導波的振幅分別為：

(21)

(22)

由式(21)和(22)可知，運用EEMD方法反演的結果與文獻[5]中EEMD反演的結果97 m相差無幾，這不僅說明了選取的剖面數據的可靠性和利用EEMD方法反演結果的有效性，同時也說明了VMD方法反演出的前導波振幅與實測資料十分吻合，且相比于EEMD方法更加接近實測數據，進而驗證了VMD反演內波參數的可行性。

3 結 語

針對EMD方法在反演內波參數方面存在模態(tài)混疊的問題，提出了一種基于VMD對SAR遙感圖像中的內波參數進行自動反演的方法。該方法先對SAR圖像進行Canny處理，獲取圖像中的內波條紋信息，再根據內波傳播方向自動選取灰度剖面；然后利用EEMD信號特征自適應分解模態(tài)函數的特點，將分解得到的有效模態(tài)數作為VMD中參數K的參考值；最后利用VMD分解后的數據進行內波參數反演。結果表明VMD方法不僅能夠有效地克服EMD方法模態(tài)混疊的問題，提取的波形符合內波物理特性，而且反演的前導波振幅與實測資料相比基本吻合，這進一步說明該反演方法是有效的。目前僅使用一種方法對K值進行估計，運用多種方法來估計K值，并對比在不同方法下參數的反演結果是下一步要繼續(xù)探討的內容。