本发明公开了一种基于ODP的海底沙波特征自动识别方法。该方法根据水深曲面的最优剖面方向求导并判定极值，以提取海底沙波脊线和谷线特征点，用于海底沙波特征线的自动识别与提取，并通过1）构建数字水深矩阵，2）构建最优方向矩阵，3）提取海底沙波特征等三大步骤实现了海底沙波的自动识别。经实测多波束水深数据测试验证，采用该方法自动识别的海底沙波特征线与人工识别特征线的相关系数可达80%以上，与传统方法相比自动识别准确率平均提升约30%，且勿需设置阈值，大幅提升了工作效率。该技术方法可推广应用到其他类型海底地形以及陆地地形的自动识别，在海洋测绘、海洋地理信息系统、计算机图形学和海底科学研究等方面具有重要的实际应用价值。
1.一种基于格网的海底沙波特征线的最优方向剖面自动识别方法，其特征在于，包括下列步骤：
1.1.1)若无多波束离散水深数据，采用全覆盖测量方式获取多波束水深数据，经数据编辑和改正处理后获取多波束离散数据集合MBES，进入步骤1.1.2)；
1.1.2)若有多波束离散水深数据集合MBES，采用距离反比加权方法构建水深矩阵其中G(i，j)为第i行第j列水深矩阵点，nx为矩阵总行数，ny为矩阵总列数，i、j、nx、ny均为自然数；
以水深矩阵 的行i为顺序，依次建立每行的数据索引Indexi＝{id[m]}m＝1，n，id[m]中存储着水深矩阵Depth第i行的第m个有水深数值的格网所在列，n表示水深矩阵Depth第i行有水深数值的格网的个数，其中m和n均为自然数；
2.1)建立最优剖面方向矩阵 令a(i，j)＝0；
2.2)建立步骤1.1.2)所述的G(i，j)的四邻域窗口矩阵
2.3)计算G(i，j)的最优剖面方向值 其中dx＝[G(i，j-1)-G(i，j+1)]，dy＝[G(i+1，j)-G(i-1，j)]，根据步骤2.3.1)至2.3.3)三种情况计算θ值：
2.3.4)将步骤2.3.1)至2.3.3)中计算的最优剖面方向值θ赋予a(i，j)＝θ；
2.4)根据步骤1.2)所述的数据索引Indexi，提取对应的G(i，j)，循环步骤2.2)至
2.3)，建立完整的最优剖面方向矩阵Aspect；
提取G(i，j)对应的最优剖面方向a(i，j)，a(i，j)如步骤2.3)所述，根据步骤3.1.1)至3.1.4)确定对应G(i，j)的求导方向：
3.1.1)若a(i，j)在区间I或区间V，则沿y＝0方向求导；
3.1.2)若a(i，j)在区间II或区间VI，则沿y＝x方向求导；
3.1.3)若a(i，j)在区间III或区间VII，则沿x＝0方向求导；
3.1.4)若a(i，j)在区间IV或区间VIII，则沿y＝-x方向求导；
其中，区间I＝(0°，22.5°]∪(337.5°，360°]，区间II＝(22.5°，67.5°]，区间III＝(67.5°，112.5°]，区间IV＝(112.5°，157.5°]，区间V＝(157.5°，202.5°]，区间VI＝(202.5°，247.5°]，区间VII＝(247.5°，292.5°]，区间VIII＝(292.5°，
3.2.2)建立极小值矩阵 令c(i，j)＝0，Crest也称为沙波脊线矩阵；
3.2.3)建立极大值矩阵 令t(i，j)＝0，Trough也称为沙波谷线矩阵；
3.2.4)根据下述步骤a)至d)来判断沙波特征：a)当G(i，j)的沿y＝0方向求导时：i)若[G(i，j)-G(i-1，j)]＜0且[G(i，j)-G(i+1，j)]＜0，则G(i，j)为极小值点，c(i，j)＝1；
ii)b)若[G(i，j)-G(i-1，j)]＞0且[G(i，j)-G(i+1，j)]＞0，则G(i，j)为极大值点，t(i，j)＝1；
b)当G(i，j)的沿y＝x方向求导时：i)若[G(i，j)-G(i-1，j+1)]＜0且[G(i，j)-G(i+1，j-1)]＜0，则G(i，j)为极小值点，c(i，j)＝1；
ii)b)若[G(i，j)-G(i-1，j+1)]＞0且[G(i，j)-G(i+1，j-1)]＞0，则G(i，j)为极大值点，t(i，j)＝1；
c)当G(i，j)的沿x＝0方向求导时：i)若[G(i，j)-G(i，j-1)]＜0且[G(i，j)-G(i，j+1)]＜0，则G(i，j)为极小值点，c(i，j)＝1；
ii)b)若[G(i，j)-G(i，j-1)]＞0且[G(i，j)-G(i，j+1)]＞0，则G(i，j)为极大值点，t(i，j)＝1；
d)当G(i，j)的沿y＝-x方向求导时：i)若[G(i，j)-G(i-1，j-1)]＜0且[G(i，j)-G(i+1，j+1)]＜0，则G(i，j)为极小值点，c(i，j)＝1；
ii)b)若[G(i，j)-G(i-1，j-1)]＞0且[G(i，j)-G(i+1，j+1)]＞0，则G(i，j)为极大值点，t(i，j)＝1；
3.3)根据步骤1.2)所述的数据索引Indexi，提取对应的G(i，j)，依次循环步骤3.1)至3.2)；
3.4)输出极小值矩阵Crest和极大值矩阵Trough。
一种基于ODP的海底沙波特征自动识别方法
技术领域
[0001] 本发明涉及海底地形地貌制图、海洋测绘、海洋地理信息系统、计算机图形学和海底科学等技术领域。
背景技术
[0002] 海底沙波是发育在近海陆架上的一种极常见的海洋地貌，沙波的迁移和运动可表征海洋沉积环境的特征，也可能影响甚至妨碍人类活动，如可能掏空或掩埋铺设的海底管线，淤浅航道港口，导致风电场的构筑物基础失稳等，可造成严重的海底灾害。因此，海底沙波的研究对了解浅海沉积物输运，合理设计规划海洋工程提供有力的信息和支持。沙波具有似波的形态，波峰和波谷以及由其在平面方向上连结而成的特征线脊线和谷线是描述其形态特征的最基本参量，也是研究沙波迁移运动规律的基础。
[0003] 陆地山谷地形自动识别方面有多重方法，如：基于地表流水物理模拟的水文分析法，基于地表几何形态的坡向分析法，以及综合运用分析法等。陆地山谷地形发育多受控于构造运动，而海底沙波发育的主因是潮流和波浪，二者在几何尺度不可类比。因此难以简单照搬这些方法进行海底沙波的识别。试验表明，应用上述方法提取海底沙波形态特征时对于地形变化平缓和局部起伏差异大的区域，特征信息提取效果不佳。此外，分类阈值是这些方法识别陆地地形特征的基础，但是目前阈值的选取还存在一些局限性，需要反复尝试且易受主观因素影响，导致地形特征识别的自动化程度和工作效率降低。
[0004] 在海底沙波研究方面，已有一些研究成果，如：“一种基于MBES的海底沙波地貌运动探测方法(ZL.201310317429.9)”公开了一种基于多波束测深数据来识别海底沙波运动的方法，“一种海底大型复杂沙波地貌的精确探测方法(ZL.201310317430.1)”公开了一种基于剖面的海底沙波FFT分解与拟合方法。但这些方法均没提及如何自动识别海底沙波的特征线，也就是海底沙波的脊线和谷线。
[0005] 因此，目前的陆地地形特征识别方法难以简单移植到海底沙波特征自动识别，对于海底沙波的形态自动识别尚缺乏有针对性的方法。
发明内容
[0006] 针对现有地形识别方法强依赖阈值选取的局限性和识别率不高的问题，本发明公开了一种基于ODP的海底沙波特征自动识别方法，具体而言是一种基于格网的海底沙波特征线的最优方向剖面自动识别方法(ODP，Optimally-Directional Profiling method)：首先基于水深曲面的最优剖面方向求导并判定极值，以提取海底沙波脊线和谷线特征点。
[0007] 一种基于ODP的海底沙波特征自动识别方法，包括下列步骤：
[0008] 1)测绘数字水深矩阵
[0009] 1.1)构建水深矩阵
[0010] 1.1.1)若无多波束离散水深数据，采用全覆盖测量方式获取多波束水深数据，经数据编辑和改正处理后获取多波束离散数据集合MBES，进入步骤1.1.2)；
[0011] 1.1.2)若有多波束离散水深数据集MBES，采用距离反比加权方法构建水深矩阵其中G(i，j)为第i行第j列水深矩阵点，nx为矩阵总行数，ny为矩阵总列数，i、j、nx、ny均为自然数；
[0012] 1.2)构建水深矩阵的数据索引
[0013] 以水深矩阵 的行i为顺序，依次建立每行的数据索引Indexi＝{id[m]}m＝1，n，id[m]中存储着水深矩阵Depth第i行的第m个有水深数值的格网所在列，n表示水深矩阵Depth第i行有水深数值的格网的个数，其中m和n均为自然数；
[0014] 该步骤的目的在于建立Indexi与水深矩阵Depth间的映射关系，以便在Depth中快速检索到有水深值的矩阵点；当通过Indexi提取水深矩阵Depth水深点时，令j＝id[m]，则可直接提取矩阵点G(i，j)；
[0015] 2)构建最优剖面方向矩阵
[0016] 2.1)建立最优剖面方向矩阵 令a(i，j)＝0；
[0017] 2.2) 建立如步骤 1.1.2) 所述的 G(i，j) 的四邻域窗口矩阵
[0018] 2.3)计算G(i，j)的最优剖面方向值 其中dx＝[G(i，j-1)-G(i，j+1)]，dy＝[G(i+1，j)-G(i-1，j)]，根据步骤2.3.1)至2.3.3)三种情况计算θ值：
[0019] 2.3.1)当dy＝0时：
[0020] a)若dx＞0，θ＝90°；
[0021] b)若dx＜0，θ＝270°；
[0022] 2.3.2)当dy＞0时：
[0023] a)若dx≥0，θ＝θ；
[0024] b)若dx＜0，θ＝360°-abs(θ)；
[0025] 2.3.3)当dy＜0时：
[0026] a)若dx≥0，θ＝180°-θ；
[0027] b)若dx＜0，θ＝180°+abs(θ)；
[0028] 2.3.4)将步骤2.3.1)至2.3.3)中计算的最优剖面方向值θ赋予a(i，j)＝θ；
[0029] 2.4)根据如步骤1.2)所述的数据索引Indexi，提取对应的G(i，j)，循环步骤2.2)至2.3)，建立完整的最优剖面方向矩阵Aspect；
[0030] 3)提取海底沙波特征线
[0031] 3.1)确定求导方向
[0032] 提取G(i，j)对应的最优剖面方向a(i，j)，根据步骤3.1.1)至3.1.4)确定对应G(i，j)的求导方向：
[0033] 3.1.1)若a(i，j)在区间I或区间V，则沿y＝0方向求导；
[0034] 3.1.2)若a(i，j)在区间II或区间VI，则沿y＝x方向求导；
[0035] 3.1.3)若a(i，j)在区间III或区间VII，则沿x＝0方向求导；
[0036] 3.1.4)若a(i，j)在区间IV或区间VIII，则沿y＝-x方向求导；
[0037] 其中，区间I＝(0°，22.5°]∪(337.5°，360°]，区间II＝(22.5°，67.5°]，区间III＝(67.5°，112.5°]，区间IV＝(112.5°，157.5°]，区间V＝(157.5°，202.5°]，区间VI＝(202.5°，247.5°]，区间VII＝(247.5°，292.5°]，区间VIII＝(292.5°，337.5°]；
[0038] 3.2)判定极值，即判定沙波脊线和谷线特征点
[0039] 3.2.1)构建G(i，j)的八邻域窗口矩阵B，：
[0040]
[0041] 3.2.2)建立极小值矩阵 令c(i，j)＝0，Crest也称为沙波脊线矩阵；
[0042] 3.2.3)建立极大值矩阵 令t(i，j)＝0，Trough也称为沙波谷线矩阵；
[0043] 3.2.4)根据下述步骤a)至d)来判断沙波特征：
[0044] a)当G(i，j)的沿y＝0方向求导时：
[0045] i)若[G(i，j)-G(i-1，j)]＜0且[G(i，j)-G(i+1，j)]＜0，则G(i，j)为极小值，即沙波脊线特征点，c(i，j)＝1；
[0046] ii)b)若[G(i，j)-G(i-1，j)]＞0且[G(i，j)-G(i+1，j)]＞0，则G(i，j)为极大值，即沙波谷线特征点，t(i，j)＝1；
[0047] b)当G(i，j)的沿y＝x方向求导时：
[0048] i)若[G(i，j)-G(i-1，j+1)]＜0且[G(i，j)-G(i+1，j-1)]＜0，则G(i，j)为极小值，即沙波脊线特征点，c(i，j)＝1；
[0049] ii)b)若[G(i，j)-G(i-1，j+1)]＞0且[G(i，j)-G(i+1，j-1)]＞0，则G(i，j)为极大值，即沙波谷线特征点，t(i，j)＝1；
[0050] c)当G(i，j)的沿x＝0方向求导时：
[0051] i)若[G(i，j)-G(i，j-1)]＜0且[G(i，j)-G(i，j+1)]＜0，则G(i，j)为极小值，即沙波脊线特征点，c(i，j)＝1；
[0052] ii)b)若[G(i，j)-G(i，j-1)]＞0且[G(i，j)-G(i，j+1)]＞0，则G(i，j)为极大值，即沙波谷线特征点，t(i，j)＝1；
[0053] d)当G(i，j)的沿y＝-x方向求导时：
[0054] i)若[G(i，j)-G(i-1，j-1)]＜0且[G(i，j)-G(i+1，j+1)]＜0，则G(i，j)为极小值，即沙波脊线特征点，c(i，j)＝1；
[0055] ii)b)若[G(i，j)-G(i-1，j-1)]＞0且[G(i，j)-G(i+1，j+1)]＞0，则G(i，j)为极大值，即沙波谷线特征点，t(i，j)＝1；
[0056] 3.3)根据如步骤1.2)所述的数据索引Indexi，提取对应的G(i，j)，依次循环步骤3.1)至3.2)；
[0057] 3.4)输出极小值矩阵Crest和极大值矩阵Trough，也是沙波脊线矩阵Crest和沙波谷线矩阵Trough。
[0058] 本发明的有益效果
[0059] 本发明提出并实现了一种基于ODP的海底沙波特征自动识别方法。通过实际测试，该方法能准确地、有效地、自动地识别海底沙波形态特征线。通过与人工识别的特征线的对比，该方法识别的海底沙波特征线与之相关系数均在80％以上，与传统方法相比自动识别准确率平均提升约30％。因此，该方法最大的优势是在保证识别的特征线准确性的前提下，避免了分类阈值选取过程的重复操作，地形自动化识别程度得到进一步提升，具有重要的实际应用价值。
附图说明
[0060] 图1本发明实施例的流程示意图；
[0061] 图2是图1中的构建最优剖面方向矩阵流程图；
[0062] 图3是图1中的提取海底沙波特征线流程图；
[0063] 图4本发明实施例中沙波区的水深矩阵；
[0064] 图5本发明实施例中沙波区的最优剖面方向矩阵；
[0065] 图6本发明实施例中提取沙波区脊线和谷线特征点与地形的映射；
[0066] 图7本发明实施例中水文分析方法提取沙波区特征点个数，以及与人工识别特征线相关系数统计图；
[0067] 图8本发明实施例中坡向变率方法提取沙波区特征点个数，以及与人工识别特征线相关系数统计图；
具体实施方式
[0068] 下面结合附图和实施例对本发明做进一步的说明。
[0069] 实施例1
[0070] 一种基于ODP的海底沙波特征自动识别方法，包括下列步骤：
[0071] 本发明的流程示意图见图1。
[0072] 步骤1)测绘数字水深矩阵
[0073] 1.1)构建水深矩阵
[0074] 1.1.1)若无多波束离散水深数据，采用全覆盖测量方式获取多波束水深数据，经数据编辑和改正处理后获取多波束离散数据集合MBES，进入步骤1.1.2)；
[0075] 1.1.2)若有多波束离散水深数据集MBES，采用距离反比加权方法构建水深矩阵其中G(i，j)为第i行第j列水深矩阵点，nx为矩阵总行数，ny为矩阵总列数，i、j、nx、ny均为自然数；
[0076] 1.2)构建水深矩阵的数据索引
[0077] 以水深矩阵 的行i为顺序，依次建立每行的数据索引Indexi＝{id[m]}m＝1，n，id[m]中存储着水深矩阵Depth第i行的第m个有水深数值的格网所在列，n表示水深矩阵Depth第i行有水深数值的格网的个数，其中m和n均为自然数；
[0078] 步骤2)构建最优剖面方向矩阵
[0079] 2.1)建立最优剖面方向矩阵 令a(i，j)＝0；
[0080] 2.2) 建立如步骤 1.1.2) 所述的 G(i，j) 的四邻域窗口矩阵
[0081] 2.3)计算G(i，j)的最优剖面方向值 其中dx＝[G(i，j-1)-G(i，j+1)]，dy＝[G(i+1，j)-G(i-1，j)]，根据步骤2.3.1)至2.3.3)三种情况计算θ值：
[0082] 2.3.1)当dy＝0时：
[0083] a)若dx＞0，θ＝90°；
[0084] b)若dx＜0，θ＝270°；
[0085] 2.3.2)当dy＞0时：
[0086] a)若dx≥0，θ＝θ；
[0087] b)若dx＜0，θ＝360°-abs(θ)；
[0088] 2.3.3)当dy＜0时：
[0089] a)若dx≥0，θ＝180°-θ；
[0090] b)若dx＜0，θ＝180°+abs(θ)；
[0091] 2.3.4)将步骤2.3.1)至2.3.3)中计算的最优剖面方向值θ赋予a(i，j)＝θ；
[0092] 2.4)根据如步骤1.2)所述的数据索引Indexi，提取对应的G(i，j)，循环步骤2.2)至2.3)，建立完整的最优剖面方向矩阵Aspect；
[0093] 步骤3)提取海底沙波特征线
[0094] 3.1)确定求导方向
[0095] 提取G(i，j)对应的最优剖面方向a(i，j)，根据步骤3.1.1)至3.1.4)确定对应G(i，j)的求导方向：
[0096] 3.1.1)若a(i，j)在区间I或区间V，则沿y＝0方向求导；
[0097] 3.1.2)若a(i，j)在区间II或区间VI，则沿y＝x方向求导；
[0098] 3.1.3)若a(i，j)在区间III或区间VII，则沿x＝0方向求导；
[0099] 3.1.4)若a(i，j)在区间IV或区间VII，则沿y＝-x方向求导；
[0100] 其中，区间I＝(0°，22.5°]∪(337.5°，360°]，区间II＝(22.5°，67.5°]，区间III＝(67.5°，112.5°]，区间IV＝(112.5°，157.5°]，区间V＝(157.5°，202.5°]，区间VI＝(202.5°，247.5°]，区间VII＝(247.5°，292.5°]，区间VIII＝(292.5°，337.5°]；
[0101] 3.2)判定极值，即判定沙波脊线和谷线特征点
[0102] 3.2.1)构建水深数据G(i，j)的八邻域窗口矩阵B，：
[0103]
[0104] 3.2.2)建立极小值矩阵 令c(i，j)＝0，Crest也称为沙波脊线矩阵；
[0105] 3.2.3)建立极大值矩阵 令t(i，j)＝0，Trough也称为沙波谷线矩阵；
[0106] 3.2.4)根据下述步骤a)至d)来判断沙波特征：
[0107] a)当G(i，j)的沿y＝0方向求导时：
[0108] i)若[G(i，j)-G(i-1，j)]＜0且[G(i，j)-G(i+1，j)]＜0，则G(i，j)为极小值，即沙波脊线特征点，c(i，j)＝1；
[0109] ii)b)若[G(i，j)-G(i-1，j)]＞0且[G(i，j)-G(i+1，j)]＞0，则G(i，j)为极大值，即沙波谷线特征点，t(i，j)＝1；
[0110] b)当G(i，j)的沿y＝x方向求导时：
[0111] i)若[G(i，j)-G(i-1，j+1)]＜0且[G(i，j)-G(i+1，j-1)]＜0，则G(i，j)为极小值，即沙波脊线特征点，c(i，j)＝1；
[0112] ii)b)若[G(i，j)-G(i-1，j+1)]＞0且[G(i，j)-G(i+1，j-1)]＞0，则G(i，j)为极大值，即沙波谷线特征点，t(i，j)＝1；
[0113] c)当G(i，j)的沿x＝0方向求导时：
[0114] i)若[G(i，j)-G(i，j-1)]＜0且[G(i，j)-G(i，j+1)]＜0，则G(i，j)为极小值，即沙波脊线特征点，c(i，j)＝1；
[0115] ii)b)若[G(i，j)-G(i，j-1)]＞0且[G(i，j)-G(i，j+1)]＞0，则G(i，j)为极大值，即沙波谷线特征点，t(i，j)＝1；
[0116] d)当G(i，j)的沿y＝-x方向求导时：
[0117] i)若[G(i，j)-G(i-1，j-1)]＜0且[G(i，j)-G(i+1，j+1)]＜0，则G(i，j)为极小值，即沙波脊线特征点，c(i，j)＝1；
[0118] ii)b)若[G(i，j)-G(i-1，j-1)]＞0且[G(i，j)-G(i+1，j+1)]＞0，则G(i，j)为极大值，即沙波谷线特征点，t(i，j)＝1；
[0119] 3.3)根据如步骤1.2)所述的数据索引Indexi，提取对应的G(i，j)，依次循环步骤3.1)至3.2)；
[0120] 3.4)输出极小值矩阵Crest和极大值矩阵Trough，也是沙波脊线矩阵Crest和沙波谷线矩阵Trough。
[0121] 实施例2
[0122] 为验证“一种基于ODP的海底沙波特征自动识别方法”的有效性和正确性，使用实施例1中之技术流程，选取一较有代表性沙波区进行实验，具体流程：
[0123] 1)测绘数字水深矩阵：在实验中使用多波束测深系统实测的沙波数据进行方法的检验。高精度水深数据采用R2Sonic2024超高分辨率多波束测深系统，经数据编辑和改正处理后获取多波束离散数据集合，采用距离反比加权方法构建沙波区的水深矩阵为图4；
[0124] 2)构建最优剖面方向矩阵：按如图2所示步骤，构建的沙波区的最优剖面方向矩阵为图5；
[0125] 3)提取海底沙波特征线：按如图3所示步骤，根据水深数据点的最优剖面方向确定求导方向，判定极值并输出沙波区脊线和谷线特征点为图6。
[0126] 4)与人工识别特征线的对比：将该方法提取的沙波区脊线和谷线与人工识别特征线之间求得相关系数分别为0.84和0.81。
[0127] 实验表明，采用本方法自动识别的海底沙波特征线与人工识别特征线的相关系数可达80％以上。
[0128] 实施例3
[0129] 为进一步验证“一种基于ODP的海底沙波特征自动识别方法”的有效性和正确性，采用水文分析方法和坡向变率分析方法对实施例2中相同沙波区进行对比实验，分别统计了在不同分类阈值条件下，水文分析方法、坡向变率分析方法提取沙波区特征点个数以及与人工识别特征线的相关系数：
[0130] 1)水文分析方法提取沙波区脊线和谷线特征点个数，以及与人工识别特征线的相关系数如图7，其中水文分析方法提取沙波区脊线和谷线与人工识别特征线的相关系数最大值分别为0.76与0.62；
[0131] 2)坡向变率分析方法提取沙波区脊线和谷线特征点个数，以及与人工识别方法的相关系数如图8，其中坡向变率分析方法提取沙波区脊线和谷线与人工识别特征线的相关系数最大值分别为0.66与0.51。
[0132] 实验表明，采用本方法与传统方法相比自动识别准确率平均提升约30％。本方法在保证识别的特征线准确性的前提下，避免了分类阈值选取过程的重复操作，地形自动化识别程度得到进一步提升，具有重要的实际应用价值。
