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

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

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

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

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

1 方法

1.1 二維EMD方法

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

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

即求平面上所有的值比周圍緊臨點(diǎn)都大或都小的點(diǎn)。對(duì)于邊界數(shù)據(jù)的處理,因?yàn)槠渲挥幸话豚徲蛴袛?shù)據(jù),所以只能在1/2的鄰域區(qū)間中尋找極值點(diǎn)。

1.1.2 平面剖分

極值點(diǎn)找出來(lái)后,它們?cè)谄矫嫔鲜巧y分布的,需要把它們按一定的鄰接關(guān)系有序地組織起來(lái),以便在空間上進(jìn)行曲面擬合。本文采用的是 Delaunay三角剖分。

1.1.3 包絡(luò)面的擬合

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

1.1.4 圖像尺度分離

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

式中,Di是尺度分離后得到的不同尺度圖層,Rn是得到的最終趨勢(shì)項(xiàng)。

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

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

作為隨機(jī)海浪的一個(gè)重要統(tǒng)計(jì)參數(shù),譜可以直接給出海浪組成波的能量相對(duì)于頻率、方向以及空間的分布。斜率波數(shù)譜表示波能量在時(shí)間和空間上的平均分布。海面微結(jié)構(gòu)光學(xué)測(cè)量裝置實(shí)驗(yàn)獲取的圖像包含了海面微尺度波的斜率信息。對(duì)圖像進(jìn)行處理,得到波面上每一點(diǎn)的順風(fēng)方向和側(cè)風(fēng)方向的兩個(gè)斜率分量,進(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的兩個(gè)正交的分量,且k=|K|=。

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

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

其中C為常數(shù)。

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

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

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

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

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

2 實(shí)驗(yàn)

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

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

圖1 風(fēng)速6.6 m/s時(shí)獲取的海面微結(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)速時(shí)毛細(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ù)慢變過(guò)程的拐點(diǎn)。比較國(guó)外類似研究,發(fā)現(xiàn)和美國(guó) Scripps海洋研究所Zhang Xin的結(jié)果[7]完全一致。另外,大約4 rad/cm時(shí)出現(xiàn)以大約波數(shù)的–3次方下降現(xiàn)象(圖 3中的虛線所示斜率)。這也與已有研究成果相似[7]。

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

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

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

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

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

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

3 結(jié)論與討論

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

二維 EMD方法對(duì)長(zhǎng)重力波的濾除效果目前還主要從斜率波數(shù)譜計(jì)算結(jié)果圖中定性得到,隨著實(shí)驗(yàn)裝置能力的提升和研究的深入,需要提出一種更直接的評(píng)價(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ì)算[J].電子學(xué)報(bào),2009,37(7):1532-1535.

[4]宋平艦,張杰.二維經(jīng)驗(yàn)?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.