CN110532615B - 一种浅海复杂地形的逐级分解方法 - Google Patents

Abstract

本发明公开了一种浅海复杂地形的逐级分解方法。在浅海陆架上，多种尺度的海底地形叠加形成复杂的地貌体系，导致难以进行量化研究，本发明利用傅里叶分析方法，结合Butterworth滤波器，实现地形信号由空间域到频率域的转换，通过分频将叠加地形分解成几组具有不同频率的单一类型地形。实例测试结果表明：该方法可快速、准确地将复杂海底叠加地形分解成不同空间尺度的地形。本发明在海洋测绘、海洋工程的安全评估和海洋科学研究等方面具有重要实用价值。

Description

一种浅海复杂地形的逐级分解方法

技术领域

本发明涉及一种浅海复杂地形的逐级分解方法。

背景技术

陆架浅海是陆地向外海延伸的海域，其中常见的地形有浅滩、沙脊、沙波、海底水道等，在复杂的海底水动力的作用之下，海底地形的形态往往以多尺度糅杂在一起为主，给不同尺度地形的定量化分析带来极大困难。

随着多波束测深技术在海底地貌探测中的应用，实现了由点到面的水下测量方式，推动了海底地貌学由定性向定量化研究的发展，在陆架残留沙脊和地貌界限研究等方面发挥了重要作用。Ronald R.Gutierrez等提出将海底沙波的命名和符号标准化，采用连续小波变换和鲁棒样条滤波器相结合的方法对浅海沙波地形进行了定性分类。Knaapen尝试使用低通滤波消除沙波地貌中的沙波纹信号，仅定量计算大型沙波的形态特征和迁移速率。由于多尺度底床形态共生的特性，高精度多波束测深数据中叠加有丰富的不同尺度地形的形态参数信息，如波高、波长、走向等，然而现有的分析方法无法对这些信息进行充分利用。

目前，基于实测资料的复合型地貌研究，仍以多尺度地形复合在一起作为整体进行的定性分析为主，对每种尺度的具体量化尚显不足。常用的频谱分析的方法有傅里叶变换，小波变换等，应用在海底地形信号处理方面多以降噪为主，而对于多尺度地形的分解研究鲜有出现。

发明内容

为了解决海底复杂地形的难以量化和研究的问题，本发明公开了一种浅海复杂地形的逐级分解方法。

本发明通过以下具体技术方案实现：

1)数据预处理：对勘测获取的原始多波束水深数据进行CUBE滤波、声速改正、潮位改正数据处理过程，构建出高精度数字水深模型；

2)频率分析：采用一维和二维傅里叶变换将叠加地形由空间域转换到频率域；

一维傅里叶变换：基于所构建的高精度数字水深模型，选取特征地形剖面，形成地形剖面水深数据集，对离散水深信号做一维傅里叶变换，得到地形剖面的一维频率信号数据集；

二维傅里叶变换：基于高精度数字水深模型构建二维水深矩阵，进行二维离散傅里叶变换，得到叠加地形的二维频率信号矩阵；

3)分频：对步骤2)中转换得到的频率域地形信号进行频率分析，确定信号截止频率，通过Butterworth滤波器进行低通滤波，分离出包含不同频率的地形信号，经傅里叶反变换，得到包含不同尺度地形的一维水深剖面和二维灰度矩阵，从而分解出不同频率的地形；

4)灰度变换：将步骤3)中分解得到的不同频率地形的水深值映射为灰度值进行显示，采用直方图均衡算法进行灰度变换，使分解后的地形具有更加的显示效果。

步骤2)中，所述的一维傅里叶变换，空间域地形信号中，空间频率为单位距离上周期发生的次数，单位为1/m，一个有限区间长度的离散信号序列有n个采样值f(i)，i＜n，这些值由傅里叶级数表达为公式(1)：

式中F(k)是傅里叶变换计算得到的傅里叶级数，e为自然常数。

步骤3)中，所述的分频：

地形信号通过Buterworth滤波器进行低通滤波，将三种频率信号进行分解，传递函数为公式(2)：

其中，(i,j)为水深矩阵中的点位置坐标，n为滤波器的阶数，d₀为低通滤波器的截止频率，d(i,j)为到中心原点的距离；

傅里叶变换将空间域转换为频率域后，频率信号v(i,j)由不同频率的地形信号组成，如公式(3)所示：

v(i,j)＝v_l(i,j)+v_m(i,j)+v_h(i,j) (3)

其中，v_l(i,j)、v_m(i,j)和v_h(i,j)分别为低频、中频和高频分量；

将转换后的地形信号通过Butterworth滤波器h(i,j)，即传递函数与傅里叶变换后的频率域信号v(i,j)做卷积运算，得到低频、中频和高频信号，从而实现不同尺度地形的分离。

步骤4)中，所述的直方图均衡算法：

将分解得到的地形水深值映射到灰度值后，得到地形的灰度图像，计算图像的直方图分布概率，进而计算出累计灰度级直方图各项，如公式(5)所示：

其中，S_k为累计直方图各项，n_j表示第j个灰度在整个地形图像中像素点个数，n为像素点的总数；

根据累计直方图概率值求得灰度转换映射关系，如公式(6)所示：

V(k)＝int[S_k(max(Z_k)-min(Z_k))+0.5] (6)

其中，Z_k为原图像灰度，V_k为均衡化图像灰度；

基于上述映射关系将原图像转化到均衡化图像灰度，使图像得到增强，直方图均衡算法使图像直方图得到延展，均衡化后图像的灰度级范围更宽，有效地增强了图像的对比度。

本发明有益效果是：本发明基于傅里叶变换和Buterworth滤波器相结合的方法，快速、准确地将浅海复杂叠加地形分解成不同尺度的单一类型地形，为后续进行地形特征的分析研究提供了较好的数据基础。

附图说明

图1是本发明的地形逐级分解方法流程图；

图2(a)是原始复杂地形平面图(2-D)；

图2(b)是原始复杂地形剖面图(1-D)

图3(a)是低频分解结果；

图3(b)是中频分解结果；

图3(c)是高频分解结果；

图4是沙波地形特征参数定义方式；

图5是沙波峰线散点图。

具体实施方式

以下结合附图对本发明做进一步的阐述。

实施例1一种浅海复杂地形的逐级分解方法实施流程

如图1所示，本实施例描述的是基于傅里叶变换和Butterworth滤波器原理相结合的一种浅海复杂地形逐级分解方法，包括以下主要步骤：

1)数据预处理：对勘测获取的原始多波束水深数据进行CUBE滤波、声速改正、潮位改正等数据处理过程，构建出高精度数字水深模型；

2)频率分析：采用一维和二维傅里叶变换将叠加地形由空间域转换到频率域。

空间域地形信号中，空间频率为单位距离上周期发生的次数，单位为1/m，一个有限区间长度的离散信号序列有n个采样值f(i)，i＜n，这些值由傅里叶级数表达为公式(1)：

式中F(k)是傅里叶变换计算得到的傅里叶级数，e为自然常数。

一维傅里叶变换：基于所构建的高精度数字水深模型，选取特征地形剖面，形成地形剖面水深数据集，对离散水深信号做一维傅里叶变换，得到地形剖面的一维频率信号数据集；

二维傅里叶变换：基于高精度数字水深模型构建二维水深矩阵，进行二维离散傅里叶变换，得到叠加地形的二维频率信号矩阵；

3)分频：对一维和二维的频率信息进行分析，确定信号截止频率，通过Butterworth滤波器的低通滤波分离出包含不同频率(尺度)的沙波信号，对其分别进行傅里叶反变换，得到包含不同尺度沙波的一维水深剖面和二维灰度矩阵，从而分解出不同频率的沙波；

地形信号通过Buterworth滤波器进行低通滤波后，将三种频率信号进行分解，传递函数为公式(2)：

其中，(i,j)为水深点坐标，n为滤波器的阶数，d₀为低通滤波器的截止频率，d(i,j)为到中心原点的距离；

傅里叶变换将空间域转换为频率域后，频率信号v(i,j)由不同频率的地形信号组成，如公式(3)所示：

v(i,j)＝v_l(i,j)+v_m(i,j)+v_h(i,j) (3)

其中，v_l(i,j)、v_m(i,j)和v_h(i,j)分别为低频、中频和高频分量；

将图像通过Butterworth滤波器h(i,j)，即传递函数与傅里叶变换后的频率域信号v(i,j)做卷积运算，得到低频、中频和高频信号，从而实现不同尺度地形的分离。

4)灰度变换：将步骤3)中分解得到的不同频率地形的水深值映射为灰度值进行显示，得到地形的灰度图像，计算图像的直方图分布概率，进而计算出累计灰度级直方图各项，如公式(5)所示：

其中，S_k为累计直方图各项，n_j表示第j个灰度在整个地形图像中像素个数，n为整个图像像素点的总数；

根据累计直方图概率值求得灰度转换映射关系，如公式(6)所示：

V(k)＝int[S_k(max(Z_k)-min(Z_k))+0.5] (6)

其中，Z_k为原图像灰度，V_k为均衡化图像灰度。

采用直方图均衡算法进行灰度变换，使分解后的地形具有更加的显示效果。

实施例2以浅海典型沙波地形为例的具体应用

以图2(a)、2(b)所示数据为例，对该方法进行验证分析。所示数据为多波束勘测得到的全覆盖水深数据，分辨率为1×1m，图2(a)为平面地形，图2(b)截取的一条沙波剖面，图中显示出多种尺度的沙波地形叠加在一起形成复杂的浅海地形。按照实施例1中的方法流程，经傅里叶变换将地形信号转换到频率域，通过分频处理逐级分解成具有不用频率的沙波地形。经灰度变换后得到地形灰度图像，图3(a)-(c)为最终分解得到的三中类型的地形灰度图像。

图4为沙波形态参数的定义方式，波长为相邻波峰之间的水平距离，波高为波峰到两个相邻波谷连线的垂直距离。基于分离出的不同尺度的沙波地形数据，沿剖面方向求取极值点，全面识别了研究区的沙波峰值点，参见附图5，图中深色点为提取的巨型沙波波峰点，浅色点为中型沙波波峰点，此外通过线性拟合可计算出沙波的平均走向。

基于分解出的地形和提取的波峰信息，统计出了试验区沙波地形的形态参数特征。表1显示了分离后三种尺度沙波主要几何参数的统计结果，需要说明的是，小型沙波纹的数目众多且走向规律性不太明显，因此未对其进行统计。其中在实验区分离出的尺度最大的沙波，平均波高超过5m，属于巨型沙波；叠加在巨型沙波上尺度次一级的中型沙波平均波高3m左右，波长在100m以内；傅里叶变换产生大量傅里叶级数，对原始水深剖面分离结果表现出较高的分辨率，分离出的小尺度沙波和波纹数目众多，一般波高小于0.4m，平均波长在5m以内。

表1三种尺度沙波主要几何参数的统计结果

Claims (5)

1.一种浅海复杂地形的逐级分解方法，其特征在于包括以下步骤：

1)数据预处理：对勘测获取的原始多波束水深数据进行CUBE滤波、声速改正、潮位改正数据处理过程，构建出高精度数字水深模型；

2)频率分析：采用一维和二维傅里叶变换将叠加地形由空间域转换到频率域；

所述的一维傅里叶变换：基于所构建的高精度数字水深模型，选取特征地形剖面，形成地形剖面水深数据集，对离散水深信号做一维傅里叶变换，得到地形剖面的一维频率信号数据集；

所述的二维傅里叶变换：基于高精度数字水深模型构建二维水深矩阵，进行二维离散傅里叶变换，得到叠加地形的二维频率信号矩阵；

3)分频：对步骤2)中转换得到的频率域地形信号进行频率分析，确定信号截止频率，通过Butterworth滤波器进行低通滤波，分离出包含不同频率的地形信号，经傅里叶反变换，得到包含不同尺度地形的一维水深剖面和二维灰度矩阵，从而分解出不同频率的地形。

2.如权利要求1所述的一种浅海复杂地形的逐级分解方法，其特征在于，步骤2)中，所述的一维傅里叶变换，空间域地形信号中，空间频率为单位距离上周期发生的次数，单位为1/m，一个有限区间长度的离散信号序列有n个采样值f(i)，i＜n，这些值由傅里叶级数表达为公式(1)：

式中F(k)是傅里叶变换计算得到的傅里叶级数，e为自然常数。

3.如权利要求1所述的一种浅海复杂地形的逐级分解方法，其特征在于，步骤3)中，所述的分频：

地形信号通过Buterworth滤波器进行低通滤波，将三种频率信号进行分解，传递函数为公式(2)：

其中，(i,j)为水深矩阵中的点位置坐标，n为滤波器的阶数，d₀为低通滤波器的截止频率，d(i,j)为到中心原点的距离；

傅里叶变换将空间域转换为频率域后，频率信号v(i,j)由不同频率的地形信号组成，如公式(3)所示：

v(i,j)＝v_l(i,j)+v_m(i,j)+v_h(i,j) (3)

其中，v_l(i,j)、v_m(i,j)和v_h(i,j)分别为低频、中频和高频分量；

将转换后的地形信号通过Butterworth滤波器h(i,j)，即传递函数与傅里叶变换后的频率域信号v(i,j)做卷积运算，得到低频、中频和高频信号，从而实现不同尺度地形的分离。

4.如权利要求1所述的一种浅海复杂地形的逐级分解方法，其特征在于，

4)灰度变换：将步骤3)中分解得到的不同频率地形的水深值映射为灰度值进行显示，采用直方图均衡算法进行灰度变换，使分解后的地形具有更加的显示效果。

5.如权利要求4所述的一种浅海复杂地形的逐级分解方法，其特征在于，

步骤4)中，所述的直方图均衡算法：

将分解得到的地形水深值映射到灰度值后，得到地形的灰度图像，计算图像的直方图分布概率，进而计算出累计灰度级直方图各项，如公式(5)所示：

其中，S_k为累计直方图各项，n_j表示第j个灰度在整个地形图像中像素点个数，n为像素点的总数；

根据累计直方图概率值求得灰度转换映射关系，如公式(6)所示：

V(k)＝int[S_k(max(Z_k)-min(Z_k))+0.5] (6)

其中，Z_k为原图像灰度，V_k为均衡化图像灰度；

基于上述映射关系将原图像转化到均衡化图像灰度，使图像得到增强，直方图均衡算法使图像直方图得到延展，均衡化后图像的灰度级范围更宽，有效地增强了图像的对比度。

Images