一种基于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；所述步骤二具体过程为：(2‑1)计算高程域G中每一点的λ尺度邻域的梯度卷积并赋值给当前点，并保存为空间尺度为λ的地形域图层公式为：$L_{\lambda }\left(\overrightarrow{X}\right)=\int {\int }_G{e}^{-\frac{|\overrightarrow{X}-\overrightarrow{X^{\prime }}{|}^2}{2{\lambda }^2}}▿G\left(\overrightarrow{X^{\prime }}\right)·\frac{\overrightarrow{X^{\prime }}-\overrightarrow{X}}{|\overrightarrow{X^{\prime }}-\overrightarrow{X}|}d\overrightarrow{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'作为新的基点，重复Step 1‑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范围内的所有点对其边界进行拟合并计算其形状参数，计算公式和参数如下所示：$\overrightarrow{C}=\frac{\int {\int }_U\overrightarrow{X}d\overrightarrow{X}}{\int {\int }_{U}d\overrightarrow{X}}$      式3$r\left(\theta \right)=r_0·(1+{\Sigma }_{i=1}^n{A}_{i}sini\theta +B_i\cos i\theta )$      式4${\mathrm{\pi r}}_0^2=\int {\int }_{X}d\overrightarrow{x};A_n=\frac{\int {\int }_X\sin n\theta d\overrightarrow{X}}{\int {\int }_{X}d\overrightarrow{X}};B_n=\frac{\int {\int }_X\cos n\theta d\overrightarrow{X}}{\int {\int }_{X}d\overrightarrow{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(θ)，对应的边界点坐标为：