CN103927543B - 一种基于dem的月表撞击坑自动识别和边界提取方法 - Google Patents

Abstract

本发明公开了一种基于DEM的月表撞击坑自动识别和边界提取方法，包括以下步骤：步骤一：对月表DEM高程域进行流域分析，进行流域划分生成；步骤二：对DEM高程域进行空间变换，将处理后的月表高程域变为人工地形域，来识别不同尺度的凹陷区域；步骤三：使用邻域扩充算法对凹陷区域进行边缘扩展，得到完整的撞击坑区域；步骤四：使用傅里叶级数拟合撞击坑区域的边界，提取准确的撞击坑形态。解决了由于撞击坑重叠覆盖造成的复杂月表难以全部识别的问题，保证识别结果的完整性；并对撞击坑对象的边界使用傅里叶级数进行拟合，提取出撞击坑的真实形状参数，保证提取结果的准确性。

Description

一种基于DEM的月表撞击坑自动识别和边界提取方法

技术领域

本发明涉及计算机图形学和数字地形分析领域，具体涉及一种基于DEM的月表撞击坑自动识别和边界提取方法。

背景技术

环形构造，特别是撞击坑，是月球表面最显著的地形特征。对撞击坑的分布特征、形态结构和数量的分析研究有着非常重要科学意义和应用价值。

早期撞击坑矢量化主要采取人工识别的方法，现今多采用自动或半自动识别方法，主要包括基于形态拟合、基于机器学习以及基于地学信息的方法：一，撞击坑外形多为圆形和椭圆型，基于形态拟合主要使用霍夫变换、二次曲线拟合以及模板匹配等方法来识别撞击坑形态。二，机器学习是研究计算机模拟或实现人类的学习活动，该方法被引入到撞击坑的自动识别过程中。比如，利用支持向量机构建模式分类器的方法，利用面向对象方法和基于神经网络等方法对月球影像数据进行撞击坑识别研究。三，基于地学信息的方法主要通过计算地形曲率和梯度来识别撞击坑。比如，坡度变化值高的地区是撞击坑的边缘，通过构建撞击坑3D模型，模拟地形剖面分析识别撞击坑。

以上方法虽然能够根据月表影像和地形数据自动识别出撞击坑并获得其形状参数，但是其实现过程和结果仍然存在以下几个问题：

(1)不能提取撞击坑真实边界，而是用圆或椭圆近似替代。由于空间风化影响和后续撞击事件的作用，撞击坑保留下来的边界是不完整、不对称的，其形状与圆和椭圆相比有较大偏差。

(2)很少考虑撞击坑之间的重叠关系，对复杂月表区域识别效果不好。

(3)大部分使用遥感影像数据来识别和提取撞击坑。相比于影像数据，DEM是月表地形地貌的直接反应，能够真实再现撞击坑的三维空间结构，因此更加适合撞击坑识别和参数定量分析。

发明内容

为解决现有技术存在的不足，本发明公开了一种基于DEM的月表撞击坑自动识别和边界提取方法，该方法以月表DEM数据为数据源，对于月表撞击坑重叠区域，能够依次识别出不同尺度的撞击坑，保证识别结果的完整性；对撞击坑区域进行边界拟合，能够提取出撞击坑的真实形状参数，保证提取结果的准确性。月表的数字高程模型简称DEM。

为实现上述目的，本发明的具体方案如下：

一种基于DEM的月表撞击坑自动识别和边界提取方法，包括以下步骤：

步骤一：对月表高程域Z的数据进行处理，生成流域F；

步骤二：对高程域Z进行空间变换，即使用空间变换理论将处理后的月表高程域Z变为空间尺度为λ的人工地形域L_λ，依此来识别不同尺度的凹陷区域；

步骤三：使用邻域扩充算法对凹陷区域进行边缘扩展，得到完整的撞击坑区域；

步骤四：使用傅里叶级数拟合撞击坑区域的边界，提取准确的撞击坑形态。

所述步骤一的具体过程为：

(1-1)，使用高斯滤波对高程域Z进行误差去除，消除DEM噪音点，生成平滑的高程域G；

(1-2)，遍历G中每一点p₁的八邻域，计算其对应八个方向的梯度，最大梯度正方向为水流方向，如果该点是其邻域的最低点，那么水流方向为其本身；

(1-3)，归属提取，遍历G中的每一点p₁，根据其水流方向追踪其流向的下一点，循环追踪直到流向终点即最低点，获取其终点位置；

(1-4)，划分流域，根据G中的每一点p₁对应的终点位置，进行流域划分，具有相同终点的点划分为同一个流域，生成流域F。

所述p₁的八邻域为p₁点的上、下、左、右、左上、左下、右上和右下八个相邻点。

所述步骤二具体过程为：

(2-1)计算高程域G中每一点的λ尺度邻域 $\left(\stackrel{\→}{X^{\′}}\right)$ 的梯度卷积并赋值给当前点，并保存为空间尺度为λ的地形域图层公式为：

$L_{\mathrm{\λ}}\left(\stackrel{\→}{X}\right)=\∫{\∫}_G{e}^{-\frac{|\stackrel{\→}{X}-\stackrel{\→}{X^{\′}}{|}^2}{2{\mathrm{\λ}}^2}}\▿G\left(\stackrel{\→}{X^{\′}}\right)\·\frac{\stackrel{\→}{X^{\′}}-\stackrel{\→}{X}}{|\stackrel{\→}{X^{\′}}-\stackrel{\→}{X}|}d\stackrel{\→}{X^{\′}}$ 式1

其中λ尺度邻域是以计算点为中心的λ×λ矩形区域，是高程域函数；

(2-2)凹陷提取。

所述凹陷提取的具体过程为：

a.四方向二次求导，对地形域图层上的每一点计算其二次导数；

b.二值化，遍历每一点，如果该点四方向二次求导结果都大于0，则判断这一点为凹陷区域点，根据是否凹陷将G二值化，生成一张标示凹陷区和非凹陷区的图像；

c.凹陷计数，使用八邻域种子填充算法，从每个凹陷点开始，查找其周围是否存在上凹点，并将其归为同一个凹陷区，生成凹陷域D；

d.凹陷裁剪。

所述凹陷裁剪的具体过程为：将凹陷域D与流域F进行比较，对空间位置不在同一流域的凹陷进行裁剪删除；假设d为D中的一个凹陷，p₂为d中的一点，其在流域F中属于流域f，f的最低点为p₂'，如果D(p₂')不等于D(p₂)，则将该点从凹陷中剔除，生成裁剪后的凹陷域TD。

所述四方向二次求导，对地形域图层上的每一点计算其二次导数的具体过程为：

假设地形域中每一点为单位正方形，对地形域图层的点L(i,j)的一次求导公式如下式2所示：

根据以上公式再对四方向求导，即得到点L(i,j)的四方向二次导数：Dx(Dx(L))，Dy(Dy(L))，Dxy(Dxy(L))和Dyx(Dyx(L))。

所述四方向包括水平X，垂直Y，对角线XY和YX。

所述步骤三中，包括一个邻域扩展凹陷区域边缘得到撞击坑区域的算法，边缘扩展的算法具体分为四步：

Step1，从凹陷每一基点a开始，如果扩展能力a(s)>Ts(设定的扩展能力阈值)，则有能力扩展，继续下一步，否则，跳出；

Step2，对于基点a八个方向即邻域上的点a_i，判断a_i全部满足条件，则a_i可以作为扩展点a'；

Step3，找到a'点后，对a'点进行扩展参数赋值；

Step4，将新增的每个扩展点a'作为新的基点，重复step1-3，直到没有扩展点产生。

所述Step2中满足条件的条件包括以下内容：

(2-11)其他原始凹陷TD；

(2-22)扩展高度d(z)>0，保证向高地势即向原凹陷上方扩展；

(2-33)点a_i与点a连线的向量的方向，与a_i点的梯度方向基本相反，即满足公式：两者之间的角度在[138.6,221.4]之间，这样保证a'大约延梯度最高的方向扩展；

(2-44)如果a_i之前扩展过，即a_i点可能满足多个凹陷的扩展条件，则保证离当前凹陷基点a的距离更近，即D_is>a_i(d)，其中D_is为之前扩展时a_i与其扩展基点之间的距离；

(2-55)将a_i赋值给a'，即a_i点可以被选择为扩展点，a'＝a_i；

所述Step3中具体过程为：

(3-11)首先对a'的扩展能力a'(s)进行计算，如果其扩展梯度a'(g)较小，即坡度较缓，则降低其扩展能力，如果其扩展梯度大，即沿着陡坡扩展的话，则增大；

当a'(g)<0.5*a(g)，则a'(s)＝a'(s)*a'(g)/a(g)；

当a'(g)<0.3*Mg，则a'(s)＝a'(s)*a'(g)/Mg；

当a'(g)>a(g)，则a'(s)＝1；

(3-22)赋值a'的扩展距离，即扩展点a'与a点之间的距离，D_is＝a(d)；

(3-33)比较凹陷中原有最大梯度Mg与新加入的扩展点a'梯度的大小，并将其中的大值赋给Mg，即Mg＝Max(Mg,a'(g))；

其中，边缘扩展算法中参数和含义如下：

原始凹陷：TD；

扩展凹陷：U，经过扩展后包含边缘部分的完整撞击坑区域；

基点：a，原始凹陷中的点，可以扩展；

扩展点：a'，以a点为基点8邻域搜索到可以扩展的点；

扩展能力：a(s)，描述点a的扩展能力，初始值为1，如果a(s)小于算法设定的扩展能力阀值Ts＝10^-λ/5，则该点不能向外扩展；同理，a'(s)描述a'点的扩展能力。

扩展距离：a(d)，描述点a和扩展点a'之间的距离，即a(d)＝|a'-a|，同理，a_i(d)描述a点和其邻域点a_i之间的距离；

扩展高度：d(z)，描述点a和a'的地形域高度差，即d(z)＝L_λ(a')-L_λ(a)；

扩展梯度：a(g)，描述a和a'之间的方向梯度，即a(g)＝d(z)/a(d)，同理，a_i(g)描述a与其邻域点a_i之间的方向梯度；a'(g)描述扩展点a'与a之间的方向梯度；

凹陷最大梯度：凹陷中所有点计算其梯度，其中最大的梯度赋值给Mg。

所述Step4中，撞击坑的边界用以为中心极点的极坐标形式表示，使用傅里叶级数根据撞击坑对象U范围内的所有点对其边界进行拟合并计算其形状参数，计算公式和参数如下所示：

$\stackrel{\&RightArrow;}{C}=\frac{\&Integral;{\&Integral;}_U\stackrel{\&RightArrow;}{X}d\stackrel{\&RightArrow;}{X}}{\&Integral;{\&Integral;}_{U}d\stackrel{\&RightArrow;}{X}}$ 式3

$r\left(\mathrm{\&theta;}\right)=r_0\&CenterDot;(1+{\mathrm{\&Sigma;}}_{i=1}^n{A}_i\sin i\mathrm{\&theta;}+B_i\cos i\mathrm{\&theta;})$ 式4

${\mathrm{\&pi;r}}_0^2=\&Integral;{\&Integral;}_{X}d\stackrel{\&RightArrow;}{x};A_n=\frac{\&Integral;{\&Integral;}_X\sin n\mathrm{\&theta;}d\stackrel{\&RightArrow;}{X}}{\&Integral;{\&Integral;}_{X}d\stackrel{\&RightArrow;}{X}};B_n=\frac{\&Integral;{\&Integral;}_X\cos n\mathrm{\&theta;}d\stackrel{\&RightArrow;}{X}}{\&Integral;{\&Integral;}_{X}d\stackrel{\&RightArrow;}{X}}$ 式5

其中，为撞击坑中心点位置，和分别为其中心点的x和y坐标；n为撞击坑对象U中点的个数，U＝{P₁P₂...P_n}，P_i(y)和P_i(x)分别为P_i点的x和y坐标；A_i和B_i为傅里叶级数的系数；r₀为近似边界圆的半径；θ和r(θ)分别为极坐标中的极角和极半径。对于每个U来说，每个像素的默认为单位正方形，其面积即为U中点的个数sum，所以每一点P_i对应的极角根据公式(5)可以求出A_i和B_i，直到A_n和B_n，然后根据公式(4)计算边界上的极径r(θ)，对应的边界点坐标为：

本发明的有益效果：

本发明涉及一种基于月表数字高程模型(DEM)的撞击坑自动识别和边界提取的方法。本发明使用空间变换对不同尺度的撞击坑进行分类识别，解决了由于撞击坑重叠覆盖造成的复杂月表难以全部识别的问题，保证识别结果的完整性；并对撞击坑对象的边界使用傅里叶级数进行拟合，提取出撞击坑的真实形状参数，保证提取结果的准确性，能够为后续撞击机理的研究提供准确位置和形状参数信息。本发明面向月球表面撞击坑的检测和边界提取，也可应用于其他行星和小行星表面。

附图说明

图1(a)原始月表DEM图；

图1(b)空间尺度λ＝5对应空间变换的地形域结果示意图；

图1(c)空间尺度λ＝10对应空间变换的地形域结果示意图；

图1(d)空间尺度λ＝20对应空间变换的地形域结果示意图；

图2(a)空间尺度λ＝20时凹陷检测后得到的撞击坑示意图；

图2(b)空间尺度λ＝20时边缘扩展后得到的撞击坑示意图；

图3(a)空间尺度λ＝5识别的撞击坑边界拟合结果图；

图3(b)空间尺度λ＝10识别的撞击坑边界拟合结果图；

图3(c)空间尺度λ＝20识别的撞击坑边界拟合结果图；

图3(d)撞击坑边界拟合最终汇总结果图。

具体实施方式：

下面结合附图对本发明进行详细说明：

如图1所示，原始月表DEM图，空间范围为左上角经纬度[73°,-5°]，右下角[83°,5°]。本发明对于存在嵌套重叠的撞击坑，根据其空间尺度进行分类，首先识别出尺寸最小的撞击坑，然后再依次识别出尺寸较大的撞击坑，直到所有撞击坑都被识别出来。为了将能够对不同尺度的撞击坑进行识别，使用空间变换理论将原来的高程域G变为空间尺度为λ的人工地形域L_λ，用于提取空间尺度为λ的撞击坑。识别时，设置λ循环逐步增长，使其从最小尺度逐步增长到最大尺度，每一尺度识别的撞击坑汇总为最终的对象集合。

本发明的基于DEM的月表撞击坑自动识别和边界提取方法，输入数据为月表高程域Z，包括以下步骤：

步骤一：流域生成

a.高斯滤波，为消除DEM噪音点，使用高斯滤波对高程域Z进行误差去除，生成平滑的高程域G；

b.流向计算，遍历G中每一点p₁的8邻域，计算其对应八个方向的梯度，最大梯度正方向为水流方向，如果该点是其邻域的最低点，那么水流方向为其本身；

c.归属提取，遍历G中的每一点p₁，根据其水流方向追踪其流向的下一点，循环追踪直到流向终点(即最低点)，获取终点位置赋值给p₁；

d.划分流域，根据G中的每一点p₁对应的终点位置，进行流域划分，具有相同终点的点划分为同一个流域，生成流域F。

步骤二：地形域转换

设置空间尺度λ＝5，使用式(1)对高程域G中的每一点进行计算，生成人工地形域用于提取空间尺度为λ的凹陷。

步骤三：凹陷提取

a.四方向二次求导，对地形域上的每一点计算其二次导数；

b.二值化，遍历每一点，如果该点四方向二次求导结果都大于0，则判断这一点为凹陷区域点，根据是否凹陷将G二值化，生成一张标示凹陷区和非凹陷区的图像；

c.凹陷计数，使用四边界或八边界邻域种子填充算法，从每个凹陷点开始，查找其周围是否存在凹陷点，并将其归为同一个凹陷区，并对该凹陷区进行标识，生成凹陷域D；

d.凹陷裁剪，所述凹陷裁剪的具体过程为：将凹陷域D与流域F进行比较，对空间位置不在同一流域的凹陷进行裁剪删除；假设d为D中的一个凹陷，p₂为d中的一点，其在流域F中属于流域f，f的最低点为p₂'，如果D(p₂')！＝D(p₂)，则将该点从凹陷中剔除，生成裁剪后的凹陷域TD。

步骤四：边缘扩展

使用邻域扩展算法自动对凹陷域TD中每一个对象进行边缘扩展，生成对应的完整撞击坑对象U。

步骤五：边界拟合

对每一个撞击坑对象U使用公式(4)和公式(5)进行边界拟合，生成对应的撞击坑形状边界。

步骤六：循环重做

设置提取凹陷的空间尺度λ₁＝λ*2，λ＝λ₁重复步骤二到步骤六，直到λ>20。

Claims (1)

1.一种基于DEM的月表撞击坑自动识别和边界提取方法，其特征是，包括以下步骤：

步骤一：对月表高程域Z的数据进行处理，生成流域F；

步骤二：对高程域Z进行空间变换，即使用空间变换理论将处理后的月表高程域Z变为空间尺度为λ的人工地形域L_λ，依此来识别不同尺度的凹陷；

步骤三：使用邻域扩充算法对凹陷进行边缘扩展，得到完整的撞击坑区域；

步骤四：使用傅里叶级数拟合撞击坑区域的边界，提取准确的撞击坑形态；

所述步骤一的具体过程为：

(1-1)，使用高斯滤波对高程域Z进行误差去除，消除DEM噪音点，生成平滑的高程域G；

(1-2)，遍历G中每一点p₁的八邻域，计算其对应八个方向的梯度，最大梯度正方向为水流方向，如果该点是其邻域的最低点，那么水流方向为其本身；

(1-3)，归属提取，遍历G中的每一点p₁，根据其水流方向追踪其流向的下一点，循环追踪直到流向终点即最低点，获取其终点位置；

(1-4)，划分流域，根据G中的每一点p₁对应的终点位置，进行流域划分，具有相同终点的点划分为同一个流域，生成流域F；

所述步骤二具体过程为：

(2-1)计算高程域G中每一点的λ尺度邻域的梯度卷积并赋值给当前点，并保存为空间尺度为λ的地形域图层公式为：

$L_{\mathrm{\&lambda;}}\left(\stackrel{\&RightArrow;}{X}\right)=\&Integral;{\&Integral;}_G{e}^{-\frac{|\stackrel{\&RightArrow;}{X}-\stackrel{\&RightArrow;}{X^{\&prime;}}{|}^2}{2{\mathrm{\&lambda;}}^2}}\&dtri;G\left(\stackrel{\&RightArrow;}{X^{\&prime;}}\right)\&CenterDot;\frac{\stackrel{\&RightArrow;}{X^{\&prime;}}-\stackrel{\&RightArrow;}{X}}{|\stackrel{\&RightArrow;}{X^{\&prime;}}-\stackrel{\&RightArrow;}{X}|}d\stackrel{\&RightArrow;}{X^{\&prime;}}$ 式1

其中λ尺度邻域是以计算点为中心的λ×λ矩形区域；

(2-2)凹陷提取；

所述凹陷提取的具体过程为：

a.四方向二次求导，对地形域图层上的每一点计算其二次导数；

b.二值化，遍历每一点，如果该点四方向二次求导结果都大于0，则判断这一点为凹陷区域点，根据是否凹陷将G二值化，生成一张标示凹陷区和非凹陷区的图像；

c.凹陷计数，使用八邻域种子填充算法，从每个凹陷点开始，查找其周围是否存在凹陷点，并将其归为同一个凹陷区，生成凹陷域D；

d.凹陷裁剪；

所述凹陷裁剪的具体过程为：将凹陷域D与流域F进行比较，对空间位置不在同一流域的凹陷进行裁剪删除；假设d为D中的一个凹陷，p₂为d中的一点，其在流域F中属于流域f，f的最低点为p₂'，如果D(p₂')不等于D(p₂)，则将该点从凹陷中剔除，生成裁剪后的凹陷域TD；

所述四方向二次求导，对地形域图层上的每一点计算其二次导数的具体过程为：

假设地形域中每一点为单位正方形，对地形域图层的点L(i,j)的一次求导公式如下式2所示：

式2

根据以上公式再对四方向求导，即得到点L(i,j)的四方向二次导数：Dx(Dx(L))，Dy(Dy(L))，Dxy(Dxy(L))和Dyx(Dyx(L))；

所述步骤三中，包括一个邻域扩展凹陷区域边缘得到撞击坑区域的算法，边缘扩展的算法具体分为四步：

Step1，从凹陷每一基点a开始，如果扩展能力a(s)>Ts，其中Ts为设定的扩展能力阈值，则有能力扩展，继续下一步，否则，跳出；

Step2，对于基点a八个方向即邻域上的点a_i，判断a_i全部满足条件，则a_i可以作为扩展点a'；

Step3，找到a点后，对a'点进行扩展参数赋值；

Step4，将新增的每个扩展点a'作为新的基点，重复Step1-3，直到没有扩展点产生；

所述Step2中满足条件的条件包括以下内容：

(2-11)

(2-22)扩展高度d(z)>0，保证向高地势即向原凹陷上方扩展；

(2-33)点a_i与点a连线的向量的方向，与a_i点的梯度方向基本相反，即满足公式：两者之间的角度在[138.6,221.4]之间，这样保证a'大约延梯度最高的方向扩展；

(2-44)如果a_i之前扩展过，即a_i点可能满足多个凹陷的扩展条件，则保证离当前凹陷基点a的距离更近，即D_is>a_i(d)，其中D_is为之前扩展时a_i与其扩展基点之间的距离；

(2-55)将a_i赋值给a'，即a_i点可以被选择为扩展点，a'＝a_i；

扩展距离：a(d)，描述点a和扩展点a'之间的距离，即a(d)＝|a'-a|，同理，a_i(d)描述a点和其邻域点a_i之间的距离；

扩展高度：d(z)，描述点a和a'的地形域高度差，即d(z)＝L_λ(a')-L_λ(a)；

所述Step3中具体过程为：

(3-11)首先对a'的扩展能力a'(s)进行计算，如果其扩展梯度a'(g)较小，即坡度较缓，则降低其扩展能力，如果其扩展梯度大，即沿着陡坡扩展的话，则增大；

当a'(g)<0.5*a(g)，则a'(s)＝a'(s)*a'(g)/a(g)；

当a'(g)<0.3*Mg，则a'(s)＝a'(s)*a'(g)/Mg；

当a'(g)>a(g)，则a'(s)＝1；

(3-22)赋值a'的扩展距离，即扩展点a'与a点之间的距离，D_is＝a(d)；

(3-33)比较凹陷中原有最大梯度Mg与新加入的扩展点a'梯度的大小，并将其中的大值赋给Mg，即Mg＝Max(Mg，a'(g))；

扩展梯度：a(g)，描述a和a'之间的方向梯度，即a(g)＝d(z)/a(d)，同理，a_i(g)描述a与其邻域点a_i之间的方向梯度；a'(g)描述扩展点a'与a之间的方向梯度；

所述Step4中，撞击坑的边界用以为中心极点的极坐标形式表示，使用傅里叶级数根据撞击坑对象U范围内的所有点对其边界进行拟合并计算其形状参数，计算公式和参数如下所示：

$\stackrel{\&RightArrow;}{C}=\frac{\&Integral;{\&Integral;}_U\stackrel{\&RightArrow;}{X}d\stackrel{\&RightArrow;}{X}}{\&Integral;{\&Integral;}_{U}d\stackrel{\&RightArrow;}{X}}$ 式3

$r\left(\mathrm{\&theta;}\right)=r_0\&CenterDot;(1+{\mathrm{\&Sigma;}}_{i=1}^n{A}_{i}sini\mathrm{\&theta;}+B_i\cos i\mathrm{\&theta;})$ 式4

${\mathrm{\&pi;r}}_0^2=\&Integral;{\&Integral;}_{X}d\stackrel{\&RightArrow;}{x};A_n=\frac{\&Integral;{\&Integral;}_X\sin n\mathrm{\&theta;}d\stackrel{\&RightArrow;}{X}}{\&Integral;{\&Integral;}_{X}d\stackrel{\&RightArrow;}{X}};B_n=\frac{\&Integral;{\&Integral;}_X\cos n\mathrm{\&theta;}d\stackrel{\&RightArrow;}{X}}{\&Integral;{\&Integral;}_{X}d\stackrel{\&RightArrow;}{X}}$ 式5

其中，为撞击坑中心点位置，和分别为其中心点的x和y坐标；n为撞击坑对象U中点的个数，U＝{P₁P₂...P_n}，P_i(y)和P_i(x)分别为P_i点的y和x坐标；A_i和B_i为傅里叶级数的系数；r₀为近似边界圆的半径；θ和r(θ)分别为极坐标中的极角和极半径；对于每个U来说，每个像素的默认为单位正方形，其面积即为U中点的个数sum，所以每一点P_i对应的极角根据公式(5)可以求出A_i和B_i，直到A_n和B_n，然后根据公式(4)计算边界上的极径r(θ)，对应的边界点坐标为：