宋平舰,张 杰,贾永君,王岩峰

(国家海洋局 第一海洋研究所,山东 青岛 266061)

微尺度波(Micro–scale wave)是反应海面粗糙度的主要特征之一,海面粗糙度直接影响海气界面面积等要素。因此,对海面微尺度波的研究,在海洋科学和海洋技术应用方面有重要的现实意义。这里所谓的微尺度波是指毛细波和毛细重力波,其波长范围从 1毫米到几厘米,在实际海洋现象中通常叠加于重力长波之上。因此,海面微结构光学测量装置得到的图像中不仅有毛细波和重力毛细波的信息,同时也包含有重力波的信息。微尺度波的主要几何特征是波面斜率大(陡)[1-2],并且在空间尺度上由于重力波的影响,其斜率变化很大。由于这个原因,在提取微尺度波斜率波数谱时必须将重力波滤掉,这样可以减少重力波对微尺度波谱的扰动。

由于图像的空间范围大小仅为 27.8 cm×21.5 cm,图像中只包含了海面局部区域的波面信息。区域的边界将重力波有限截断,不能包含一个完整周期内的重力波信息。这样在利用Fourier变换进行谱分析时,重力波的频谱信息混叠到其他频率的波谱中,在数学上表现为Gibbs现象。张福元[3]利用小波分析方法对实验室获取到的微尺度图像进行重力波滤除,然后对滤除后的图像进行斜率波数谱提取,但其处理结果有一个较大的先下沉后上升的明抖动,处理效果有待提高。

二维EMD(经验模态分解)方法是由宋平舰、张杰[4]在 Huang[5]提出的一维 EMD图像处理技术基础之上进行扩展后得到的,可很好地应用于二维图像的处理。这种方法能直接提取各种内部时空尺度的能量,消除数据漂移,把一幅复杂的非线性、非平稳遥感图像数据分解成若干线性、平稳的数据层,有效提取信号的不同尺度信息。本文尝试采用二维EMD方法对获取的微尺度图像进行处理,以期能够有效去除或减少长重力波对微尺度波信息提取所造成的影响,在提取斜率波数谱前将长重力波滤除。

1 方法

1.1 二维EMD方法

二维EMD图像处理方法包括四个方面的内容:寻找平面内局部极值点、平面剖分、包络面拟合和尺度分离。

1.1.1 寻找平面内局部极值点

即求平面上所有的值比周围紧临点都大或都小的点。对于边界数据的处理,因为其只有一半邻域有数据,所以只能在1/2的邻域区间中寻找极值点。

1.1.2 平面剖分

极值点找出来后,它们在平面上是散乱分布的,需要把它们按一定的邻接关系有序地组织起来,以便在空间上进行曲面拟合。本文采用的是 Delaunay三角剖分。

1.1.3 包络面的拟合

对于曲面极大值点与极小值点进行拟合,形成其包络面。这是二维EMD方法的关键。本文采用的是三次 B样条插值法。这是在剖分所提供的离散点邻接表的基础上进行的,此插值法可以保证在三角形每一条邻接边上二阶光滑,曲面拟合效果好。

1.1.4 图像尺度分离

在解决了极值包络问题后,参照一维EMD图像处理方法,对图像信息进行尺度分离,处理后得到的表达式为:

式中,Di是尺度分离后得到的不同尺度图层,Rn是得到的最终趋势项。

考虑到微结构图像是彩色图像,故将图像分层为三个通道(r,g,b)分别实施二维 EMD处理。每个通道图像在经过 EMD 处理后,将得到的Di合并,得到每个通道去除长重力波影响后的微尺度图像信息,然后再将这三个通道的微尺度图像合成为彩色图像,供斜率波数谱计算。

1.2 斜率波数谱的计算

作为随机海浪的一个重要统计参数,谱可以直接给出海浪组成波的能量相对于频率、方向以及空间的分布。斜率波数谱表示波能量在时间和空间上的平均分布。海面微结构光学测量装置实验获取的图像包含了海面微尺度波的斜率信息。对图像进行处理,得到波面上每一点的顺风方向和侧风方向的两个斜率分量,进一步利用 Fourier变换及其性质可以推导出微尺度波的斜率波数谱。

设波面位移的方程为 z=f(x,y),关于 x和 y偏导数分别为 fx(x,y)和fy(x,y),波面位移 f(x,y)的二维Fourier变换是

其中kx,ky是波数向量 K的两个正交的分量,且k=|K|=。

根据Fourier变换的微分性质可得两个垂直方向上的波面斜率 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)式可以导出F(kx,ky)。对F(kx,ky)做二维Fourier逆变换则可得到波面位移方程为

其中C为常数。

根据均匀稳定域内的波谱计算公式

其中kx,ky分别是波数矢量两个互相垂直的分量,ω是波的角频率,谱的定义域是-∞＜kx,ky,ω＜+∞。定标以后,不同位置不同时间的波面位移的自相关函数定义为

其中,上横杠表示时间平均。波数谱定义为

斜率波数谱S(k)定义为两个正交斜率波数谱的和,即

其中Sx(k,θx),Sy(k,θy),Φ(k,θx)和Φ(k,θy)是极坐标下的斜率波数谱和波数谱。

2 实验

本文实验所用测量设备为国家海洋局第一海洋研究所研制的以颜色编码原理为基础的“海面微结构光学测量装置”。实验数据是在中国海洋大学风浪水槽进行动力水槽实验时得到的含有海面微尺度波二维斜率信息的图像。

图1是获取的实验图像,此时风速 6.6 m/s;在后继处理之前需要对图像进行定标误差校正[6]。图2为使用二维 EMD处理后得到的滤除长重力波后的海面微结构图像(EMD分层至第三层),图 3是经二维EMD处理后得到的斜率波数谱。

图1 风速6.6 m/s时获取的海面微结构图像Fig.1 Original experiment image,wind speed at 6.6 m/s

图2 经二维EMD分层处理后得到的前三层合成图Fig.2 Synthesis of the first three layers by 2D-EMD

图3 6.6 m/s风速时毛细波斜率波数谱Fig.3 The slope wave number spectrum of capillary wave at 6.6 m/s wind speed

从图3中可以看出,3 rad/cm左右有明显的“dip”现象,从图像中还可以看到,在 7 rad/cm 左右这种“dip”现象消失,也就是大约 7 rad/cm是斜率波数谱从“dip”现象到以波数的幂指数慢变过程的拐点。比较国外类似研究,发现和美国 Scripps海洋研究所Zhang Xin的结果[7]完全一致。另外,大约4 rad/cm时出现以大约波数的–3次方下降现象(图 3中的虚线所示斜率)。这也与已有研究成果相似[7]。

若对定标误差校正后的图像不进行二维 EMD处理,而直接进行计算,计算结果如图4所示。

图4 使用二维EMD后的结果变化Fig.4 The effect of 2D-EMD

从图中可以看出,由于实际海洋现象中毛细波、毛细重力波和重力波同时存在,长波影响明显,且其中混叠 Gibbs现象的影响,所以未经过二维 EMD滤波处理的结果与理论推导结果相比在波数小于约5 rad/cm时加速发散,偏差较大,且没有反映出 dip现象。而经过二维EMD处理后的数据则与半经验模型所导出的理论结果较为接近[8],dip现象较为明显。在波数较高的区域(大于约 11 rad/cm),这两种结果的差异较小,这也反映出二维EMD方法在消除长波影响的同时较好地保留了微尺度波的有效信息。

为了进一步比对二维EMD方法的处理效果,本文采用小波处理算法与其进行比较。主要方法为利用多分辨小波分析对图像进行分解,可得到不同尺度信息分量。去掉低频部分分量,利用剩余的高频分量进行图像重构,最后由重构的图像计算微尺度波数谱。图5给出风速6.6 m/s下,经过四层小波变换滤波后得到的一维斜率波数谱。

在图 5中,小波处理结果有一个非常明显的特征,当波数较小的时候(波数小于3 rad/cm左右)斜率波数谱总有一个较大的先下沉后上升的抖动,且得到的斜率波数谱没有得到明显的“dip”现象。这与Zhang Xin给出的3.6 rad/cm处发生dip现象的结论是不对应的,而二维EMD方法处理后的结果就与此结论吻合较好。但在大于约 20 rad/cm 后,小波与EMD处理的结果相差不大,两种不同小波基在大于约 4 rad/cm也差异不大,这反映出存在结果差异的一个主要原因就是长波对图像信息的影响。

图5 风速6.6 m/s时的小波与二维EMD处理结果的比对Fig.5 Comparison between the results obtained by wavelet and 2D-EMD,wind speed at 6.6 m/s

3 结论与讨论

为了解决有效截断以及避免出现小波数情形下的偏差,本文采用二维EMD方法对图像中包含的重力波信息进行滤除,然后对滤除后的图像进行微尺度波斜率波数谱计算,计算结果与国内外其他研究者的结果基本一致,说明采用二维EMD方法滤除图像中混叠的长波信息,在微尺度斜率波谱计算中是有意义的。文中将二维EMD方法与小波方法的处理结果进行了比对,结果表明在滤除长波的影响时,二维EMD方法具有其较明显的优势。

二维 EMD方法对长重力波的滤除效果目前还主要从斜率波数谱计算结果图中定性得到,随着实验装置能力的提升和研究的深入,需要提出一种更直接的评价方法,定量地反映二维 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]张福元,张杰,王岩峰,等.基于小波分析的微尺度波斜率波数谱计算[J].电子学报,2009,37(7):1532-1535.

[4]宋平舰,张杰.二维经验模分解在海洋遥感图像信息分离中的应用[J].高技术通讯,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]贾永君,张杰,任启峰,等.颜色编码方法中波面斜率的校准[J].高技术通讯,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.