二維EMD方法在微尺度波斜率波數(shù)譜提取中的應(yīng)用

宋平艦,張 杰,賈永君,王巖峰

(國家海洋局 第一海洋研究所,山東 青島 266061)

微尺度波(Micro–scale wave)是反應(yīng)海面粗糙度的主要特征之一,海面粗糙度直接影響海氣界面面積等要素。因此,對海面微尺度波的研究,在海洋科學(xué)和海洋技術(shù)應(yīng)用方面有重要的現(xiàn)實意義。這里所謂的微尺度波是指毛細(xì)波和毛細(xì)重力波,其波長范圍從 1毫米到幾厘米,在實際海洋現(xiàn)象中通常疊加于重力長波之上。因此,海面微結(jié)構(gòu)光學(xué)測量裝置得到的圖像中不僅有毛細(xì)波和重力毛細(xì)波的信息,同時也包含有重力波的信息。微尺度波的主要幾何特征是波面斜率大(陡)[1-2],并且在空間尺度上由于重力波的影響,其斜率變化很大。由于這個原因,在提取微尺度波斜率波數(shù)譜時必須將重力波濾掉,這樣可以減少重力波對微尺度波譜的擾動。

由于圖像的空間范圍大小僅為 27.8 cm×21.5 cm,圖像中只包含了海面局部區(qū)域的波面信息。區(qū)域的邊界將重力波有限截斷,不能包含一個完整周期內(nèi)的重力波信息。這樣在利用Fourier變換進(jìn)行譜分析時,重力波的頻譜信息混疊到其他頻率的波譜中,在數(shù)學(xué)上表現(xiàn)為Gibbs現(xiàn)象。張福元[3]利用小波分析方法對實驗室獲取到的微尺度圖像進(jìn)行重力波濾除,然后對濾除后的圖像進(jìn)行斜率波數(shù)譜提取,但其處理結(jié)果有一個較大的先下沉后上升的明抖動,處理效果有待提高。

二維EMD(經(jīng)驗?zāi)B(tài)分解)方法是由宋平艦、張杰[4]在 Huang[5]提出的一維 EMD圖像處理技術(shù)基礎(chǔ)之上進(jìn)行擴(kuò)展后得到的,可很好地應(yīng)用于二維圖像的處理。這種方法能直接提取各種內(nèi)部時空尺度的能量,消除數(shù)據(jù)漂移,把一幅復(fù)雜的非線性、非平穩(wěn)遙感圖像數(shù)據(jù)分解成若干線性、平穩(wěn)的數(shù)據(jù)層,有效提取信號的不同尺度信息。本文嘗試采用二維EMD方法對獲取的微尺度圖像進(jìn)行處理,以期能夠有效去除或減少長重力波對微尺度波信息提取所造成的影響,在提取斜率波數(shù)譜前將長重力波濾除。

1 方法

1.1 二維EMD方法

二維EMD圖像處理方法包括四個方面的內(nèi)容:尋找平面內(nèi)局部極值點(diǎn)、平面剖分、包絡(luò)面擬合和尺度分離。

1.1.1 尋找平面內(nèi)局部極值點(diǎn)

即求平面上所有的值比周圍緊臨點(diǎn)都大或都小的點(diǎn)。對于邊界數(shù)據(jù)的處理,因為其只有一半鄰域有數(shù)據(jù),所以只能在1/2的鄰域區(qū)間中尋找極值點(diǎn)。

1.1.2 平面剖分

極值點(diǎn)找出來后,它們在平面上是散亂分布的,需要把它們按一定的鄰接關(guān)系有序地組織起來,以便在空間上進(jìn)行曲面擬合。本文采用的是 Delaunay三角剖分。

1.1.3 包絡(luò)面的擬合

對于曲面極大值點(diǎn)與極小值點(diǎn)進(jìn)行擬合,形成其包絡(luò)面。這是二維EMD方法的關(guān)鍵。本文采用的是三次 B樣條插值法。這是在剖分所提供的離散點(diǎn)鄰接表的基礎(chǔ)上進(jìn)行的,此插值法可以保證在三角形每一條鄰接邊上二階光滑,曲面擬合效果好。

1.1.4 圖像尺度分離

在解決了極值包絡(luò)問題后,參照一維EMD圖像處理方法,對圖像信息進(jìn)行尺度分離,處理后得到的表達(dá)式為:

式中,Di是尺度分離后得到的不同尺度圖層,Rn是得到的最終趨勢項。

考慮到微結(jié)構(gòu)圖像是彩色圖像,故將圖像分層為三個通道(r,g,b)分別實施二維 EMD處理。每個通道圖像在經(jīng)過 EMD 處理后,將得到的Di合并,得到每個通道去除長重力波影響后的微尺度圖像信息,然后再將這三個通道的微尺度圖像合成為彩色圖像,供斜率波數(shù)譜計算。

1.2 斜率波數(shù)譜的計算

作為隨機(jī)海浪的一個重要統(tǒng)計參數(shù),譜可以直接給出海浪組成波的能量相對于頻率、方向以及空間的分布。斜率波數(shù)譜表示波能量在時間和空間上的平均分布。海面微結(jié)構(gòu)光學(xué)測量裝置實驗獲取的圖像包含了海面微尺度波的斜率信息。對圖像進(jìn)行處理,得到波面上每一點(diǎn)的順風(fēng)方向和側(cè)風(fēng)方向的兩個斜率分量,進(jìn)一步利用 Fourier變換及其性質(zhì)可以推導(dǎo)出微尺度波的斜率波數(shù)譜。

設(shè)波面位移的方程為 z=f(x,y),關(guān)于 x和 y偏導(dǎo)數(shù)分別為 fx(x,y)和fy(x,y),波面位移 f(x,y)的二維Fourier變換是

其中kx,ky是波數(shù)向量 K的兩個正交的分量,且k=|K|=。

根據(jù)Fourier變換的微分性質(zhì)可得兩個垂直方向上的波面斜率 fx(x,y)和fy(x,y)的二維Fourier變換為

若fx(x,y)和fy(x,y)已知,Fx(kx,ky)和Fy(kx,ky)則可以通過對fx(x,y)和fy(x,y)做 Fourier變換得到,由(6)式可以導(dǎo)出F(kx,ky)。對F(kx,ky)做二維Fourier逆變換則可得到波面位移方程為

其中C為常數(shù)。

根據(jù)均勻穩(wěn)定域內(nèi)的波譜計算公式

其中kx,ky分別是波數(shù)矢量兩個互相垂直的分量,ω是波的角頻率,譜的定義域是-∞＜kx,ky,ω＜+∞。定標(biāo)以后,不同位置不同時間的波面位移的自相關(guān)函數(shù)定義為

其中,上橫杠表示時間平均。波數(shù)譜定義為

斜率波數(shù)譜S(k)定義為兩個正交斜率波數(shù)譜的和,即

其中Sx(k,θx),Sy(k,θy),Φ(k,θx)和Φ(k,θy)是極坐標(biāo)下的斜率波數(shù)譜和波數(shù)譜。

2 實驗

本文實驗所用測量設(shè)備為國家海洋局第一海洋研究所研制的以顏色編碼原理為基礎(chǔ)的“海面微結(jié)構(gòu)光學(xué)測量裝置”。實驗數(shù)據(jù)是在中國海洋大學(xué)風(fēng)浪水槽進(jìn)行動力水槽實驗時得到的含有海面微尺度波二維斜率信息的圖像。

圖1是獲取的實驗圖像,此時風(fēng)速 6.6 m/s;在后繼處理之前需要對圖像進(jìn)行定標(biāo)誤差校正[6]。圖2為使用二維 EMD處理后得到的濾除長重力波后的海面微結(jié)構(gòu)圖像(EMD分層至第三層),圖 3是經(jīng)二維EMD處理后得到的斜率波數(shù)譜。

圖1 風(fēng)速6.6 m/s時獲取的海面微結(jié)構(gòu)圖像Fig.1 Original experiment image,wind speed at 6.6 m/s

圖2 經(jīng)二維EMD分層處理后得到的前三層合成圖Fig.2 Synthesis of the first three layers by 2D-EMD

圖3 6.6 m/s風(fēng)速時毛細(xì)波斜率波數(shù)譜Fig.3 The slope wave number spectrum of capillary wave at 6.6 m/s wind speed

從圖3中可以看出,3 rad/cm左右有明顯的“dip”現(xiàn)象,從圖像中還可以看到,在 7 rad/cm 左右這種“dip”現(xiàn)象消失,也就是大約 7 rad/cm是斜率波數(shù)譜從“dip”現(xiàn)象到以波數(shù)的冪指數(shù)慢變過程的拐點(diǎn)。比較國外類似研究,發(fā)現(xiàn)和美國 Scripps海洋研究所Zhang Xin的結(jié)果[7]完全一致。另外,大約4 rad/cm時出現(xiàn)以大約波數(shù)的–3次方下降現(xiàn)象(圖 3中的虛線所示斜率)。這也與已有研究成果相似[7]。

若對定標(biāo)誤差校正后的圖像不進(jìn)行二維 EMD處理,而直接進(jìn)行計算,計算結(jié)果如圖4所示。

圖4 使用二維EMD后的結(jié)果變化Fig.4 The effect of 2D-EMD

從圖中可以看出,由于實際海洋現(xiàn)象中毛細(xì)波、毛細(xì)重力波和重力波同時存在,長波影響明顯,且其中混疊 Gibbs現(xiàn)象的影響,所以未經(jīng)過二維 EMD濾波處理的結(jié)果與理論推導(dǎo)結(jié)果相比在波數(shù)小于約5 rad/cm時加速發(fā)散,偏差較大,且沒有反映出 dip現(xiàn)象。而經(jīng)過二維EMD處理后的數(shù)據(jù)則與半經(jīng)驗?zāi)Ｐ退鶎?dǎo)出的理論結(jié)果較為接近[8],dip現(xiàn)象較為明顯。在波數(shù)較高的區(qū)域(大于約 11 rad/cm),這兩種結(jié)果的差異較小,這也反映出二維EMD方法在消除長波影響的同時較好地保留了微尺度波的有效信息。

為了進(jìn)一步比對二維EMD方法的處理效果,本文采用小波處理算法與其進(jìn)行比較。主要方法為利用多分辨小波分析對圖像進(jìn)行分解,可得到不同尺度信息分量。去掉低頻部分分量,利用剩余的高頻分量進(jìn)行圖像重構(gòu),最后由重構(gòu)的圖像計算微尺度波數(shù)譜。圖5給出風(fēng)速6.6 m/s下,經(jīng)過四層小波變換濾波后得到的一維斜率波數(shù)譜。

在圖 5中,小波處理結(jié)果有一個非常明顯的特征,當(dāng)波數(shù)較小的時候(波數(shù)小于3 rad/cm左右)斜率波數(shù)譜總有一個較大的先下沉后上升的抖動,且得到的斜率波數(shù)譜沒有得到明顯的“dip”現(xiàn)象。這與Zhang Xin給出的3.6 rad/cm處發(fā)生dip現(xiàn)象的結(jié)論是不對應(yīng)的,而二維EMD方法處理后的結(jié)果就與此結(jié)論吻合較好。但在大于約 20 rad/cm 后,小波與EMD處理的結(jié)果相差不大,兩種不同小波基在大于約 4 rad/cm也差異不大,這反映出存在結(jié)果差異的一個主要原因就是長波對圖像信息的影響。

圖5 風(fēng)速6.6 m/s時的小波與二維EMD處理結(jié)果的比對Fig.5 Comparison between the results obtained by wavelet and 2D-EMD,wind speed at 6.6 m/s

3 結(jié)論與討論

為了解決有效截斷以及避免出現(xiàn)小波數(shù)情形下的偏差,本文采用二維EMD方法對圖像中包含的重力波信息進(jìn)行濾除,然后對濾除后的圖像進(jìn)行微尺度波斜率波數(shù)譜計算,計算結(jié)果與國內(nèi)外其他研究者的結(jié)果基本一致,說明采用二維EMD方法濾除圖像中混疊的長波信息,在微尺度斜率波譜計算中是有意義的。文中將二維EMD方法與小波方法的處理結(jié)果進(jìn)行了比對,結(jié)果表明在濾除長波的影響時,二維EMD方法具有其較明顯的優(yōu)勢。

二維 EMD方法對長重力波的濾除效果目前還主要從斜率波數(shù)譜計算結(jié)果圖中定性得到,隨著實驗裝置能力的提升和研究的深入,需要提出一種更直接的評價方法,定量地反映二維 EMD的處理能力。

[1]Cox S,Munk W.Measurements of the roughness of the sea surface from photographs of the sun’s glitter[J].J Optical Soc Amer,1954,44: 838-850.

[2]Zhang Xin,Charles S C.Measuring the two-dimensional structure of a wavy water surface optically: a surface gradient detector[J].Experiments in Fluids,1994,17:225-237.

[3]張福元,張杰,王巖峰,等.基于小波分析的微尺度波斜率波數(shù)譜計算[J].電子學(xué)報,2009,37(7):1532-1535.

[4]宋平艦,張杰.二維經(jīng)驗?zāi)７纸庠诤Ｑ筮b感圖像信息分離中的應(yīng)用[J].高技術(shù)通訊,2001,11(9): 62-67.

[5]Huang N E,Shen Zheng,Long S R,et al.The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis[J].Proc R Soc Lond A ,1998,454: 903-995.

[6]賈永君,張杰,任啟峰,等.顏色編碼方法中波面斜率的校準(zhǔn)[J].高技術(shù)通訊,2008,18(1): 65-70.

[7]Zhang Xin.Wavenumber spectrum of short wind waves:an application of two-dimensional slepian windows to spectral estimation[J].Journal of Atmospheric and Oceanic Technology,1994,11: 489-505.

[8]Toba Y.Local balance in the air-sea boundary processes I.On the process of wind waves[J].Journal of the Oceanographical Society of Japan,1972,28: 109-121.