一种基于地表不变特征的矿区DEM变化检测方法与流程

本发明涉及一种dem变化检测方法，具体是一种适用于因采矿活动导致的地形变化、地表破坏严重区域的基于地表不变特征的dem变化检测方法，属于矿区地形测量与监测技术领域。

背景技术：

数字高程模型(digitalelevationmodel)简称dem，是通过有限的地形高程数据实现对地面地形的数字化模拟(即地形表面形态的数字化表达)，它是用一组有序数值阵列形式表示地面高程的一种实体地面模型，dem通常作为进行地形分析与相关科学研究的基础数据模型。

矿区地表由于受到开采扰动影响，导致地表形变灾害频发，地形变化相对复杂，然而地形图更新远赶不上地物目标的实地变化速度，利用卫星遥感影像进行dem变化检测，不仅能够缩短现有的地形测量周期，还能保证救灾防灾及安全生产对数据的及时性要求。

矿区dem是进行矿区地表变形灾害监测与环境规划治理的基础数据，常被作为高程基准对影像进行正射校正，但dem数据与进行变化检测所使用的影像通常无法同时获取，如果高程信息发生变化，它将对二维变化检测的结果产生影响。

现有地形三维变化大多采用d-insar或lidar方法，前者可以实现高精度的形变检测，但在地形起伏剧烈的区域相位解缠困难，高分辨率d-insar实施中，常使用外部dem以减少相位解缠误差与地理编码误差；后者地面观测范围有限，大范围机载lidar观测目前价格昂贵，难以大面积实施，且目前利用高分辨率卫星影像建立的dem精度还不足以达到厘米级，但是我国自主研制的资源三号(zy-3)卫星几何定位精度可达到平面优于2m，高程优于3m的定位精度，高于常被使用的astergdem(advancedspacebornethermalemissionandreflectionradiometerglobaldigitalelevationmodel)和srtm(shuttleradartopographymission)数据精度，因此利用周期性卫星立体影像研究dem的变化，不仅是三维变化检测的有力补充，还能为其他相关科学研究提供更为准确的基准数据。

传统dem变化检测通常是根据地形变化前后的dem进行求差获取，先对dem配准，然后进行变化检测，没有考虑地形变化对配准算法的影响，也没有考虑水平位置变化对dem高程的影响；由于获取途径、成像方式与时间、轨道参数、影像匹配、尺度、dem插值等差异，不同时相间的dem配准并不严格，导致dem出现虚假变化；另外，矿区内自然地形与人工地形交错，突变地形与渐变地形共生，地形复杂多变，即便是布设的控制点，也会受开采扰动影响而发生变化，影响dem与变化检测精度。

dem变化检测与dem配准密不可分，目前配准方法通常以点特征配准为主，如经常使用的最近邻点算法(iterativeclosestpoint，icp)及改进算法，以点、线相似不变性或者线特征实现城区航空影像与机载激光雷达点云的配准；也有利用虚拟控制点的方法改进传统相关系数匹配算法用于机载lidar点云和pos数据共同辅助的航空影像连接点自动匹配；或者利用点云数据中建筑物的三维轮廓线构建建筑物角特征与影像上的直线特征进行配准实现lidar点云与航空影像的配准，然而现有配准大多基于最小二乘的全局匹配，并没有考虑不同时相条件下的变化信息对配准精度的影响。

技术实现要素：

针对上述问题，本发明提供一种基于地表不变特征的矿区dem变化检测方法，通过对点、线、面等地形特征基元的相似度测量，建立适合矿区地形特点的地表不变特征判断规则，不仅能够为dem的变化检测提供准确的基准信息，且能够为dem配准法方程的建立与平差解算提供初始与基准参数，继而实现dem的变化信息解算。

为实现上述目的，本基于地表不变特征的矿区dem变化检测方法具体包括以下步骤：

a.选取矿区所在位置，收集矿区地形地貌特点、地质采矿、开采、已有dem变化区域信息；

b.基于地表不变gcp特征的矿区dem数据建模；

首先针对矿区地形dem生成与变化检测需求，根据卫星观测时间获取对应矿区内ti和tj两个不同时相的立体影像数据，建立成像几何模型；

其次在分析矿区地形特点基础上，在有明显地标物的位置处布设地面控制点，并对gcp进行不变特征判断；

最后根据成像几何模型和不变gcp，利用区域网平差技术对前视、正视、后视图像进行rpc模型精化，获取矿区大范围内基础地形dem；

c.利用数字地形分析中地形形态提取的方法和图像处理中的边缘检测方法提取步骤b中两个不同时相dem地形中的点、线、面特征基元；

d.对步骤c中两个时相的地形特征基元进行相似度测量，通过层次分析方法构建步骤c中提取的地形特征基元的几何特征测量模型，进行地形特征基元相似度测量，并依此建立地表不变特征判断规则，提取变化特征与不变特征；

e.根据步骤d确定的地表不变特征，构建不同时相dem间配准模型并进行平差解算，结合步骤d中提取的变化特征，获取dem变化；

f.将配准后的dem与基准dem直接求差获取该时相的dem变化信息，该dem变化信息与步骤d中提取的变化特征一起，共同构成该矿区内dem变化信息。

作为本发明的进一步改进方案，所述的步骤b中对gcp进行不变特征判断的方法具体包括以下步骤：

首先提取每个gcp在两期影像上的像方坐标；

其次提取并计算gcp在两期dem上所对应的平面坐标、高程值、坡度和曲率参数；

最后分别计算两期参数的差异，若差异在所设定的阈值范围之内，则该gcp为不变gcp。

作为本发明的进一步改进方案，所述的步骤d具体包括以下步骤：

①按照相同的空间分辨率尺度对两期dem数据进行尺度一致重采样，利用最近邻点算法实现两期dem地形特征基元的初始匹配；

②利用局部窗口搜索方法对范围内的地形特征基元的几何特征进行相似度测量，建立点特征相似度s1、线特征相似度s2和面特征相似度s3；

③计算同名特征的相似度s；

④通过相似度最大值选取与阈值判断结合的方式实现对地形特征的匹配。

作为本发明的进一步改进方案，所述的步骤②中点、线、面地形特征基元的相似度指标分量包括高程、水平位置、曲率、坡度、长度、弯曲度、面积、边界长度、边界形状。

作为本发明的进一步改进方案，所述的步骤③的具体方法包括：

假设特征l有n个比对特征分量项l＝[l1...ln]，则任意一项特征在两期dem上的相似度为

其中，li(i＝1,...,n)，li'(i＝1,...,n)分别为两期dem上对应的第i项特征值，则该特征的相似度为

其中，pi为对应s(li)的权值。

作为本发明的进一步改进方案，所述的步骤④具体方法包括：

根据相似度s计算结果范围为[0,1]，1代表没有变化，0代表最大变化，则设定地表不变特征的取值范围为[δs，1]，δs为相似度阈值，

作为本发明的进一步改进方案，所述的步骤e具体方法包括：

假设(x0,y0,z0)^t和(x,y,z)^t为同一地区待配准的两个dem不变特征上的点，采用区域网内平移旋转缩放模型，构建方程，

(x,y,z)^t＝t+sr(x0,y0,z0)^t

式中，s为缩放系数，r为旋转系数(ω，κ），t为平移系数(tx,ty,tz)^t；

泰勒级数展开线性化后建立m个不变特征点的观测方程，

-e＝ax-l,p

式中，a为系数针，x为待解配准参数，p为权阵，l为残差向量。联合直线特征和面特征约束观测方程：

考虑稳健解算情况下，利用迭代最小二乘法联合解算式中的待解参数。

与现有技术相比，本基于地表不变特征的矿区dem变化检测方法通过对点、线、面等地形特征基元的相似度测量，建立适合矿区地形特点的地表不变特征判断规则，不仅能够为dem的变化检测提供准确的基准信息，且能够为dem配准法方程的建立构建与平差解算提供初始与基准参数，继而实现dem的变化信息解算；本基于地表不变特征的矿区dem变化检测方法引入矿区地表不变特征建立顾及dem地形结构的配准模型，将地表不变特征和特征的几何约束用于构建dem配准模型，既可以顾及地形几何特点，又可以提高配准准确性；本基于地表不变特征的矿区dem变化检测方法利用地形特征基元相似度测量方法可实现对地表不变特征判断与变化特征的同步提取，特别适用于地形复杂多变的矿区。

附图说明

图1是本发明基于地表不变特征的dem变化检测处理过程示意图；

图2是本发明地表不变gcp特征的矿区dem数据建模流程图；

图3是本发明地形特征基元的相似度测量过程示意图。

具体实施方式

对于不同时相dem而言，必须要保证配准所使用的各类特征要素是不变的，本发明通过对点、线、面等地形特征基元的相似度测量，建立适合矿区地形特点的地表不变特征判断规则，不仅能够为dem的变化检测提供准确的基准信息，且能够为dem配准法方程的建立与平差解算提供初始与基准参数，继而实现dem的变化信息解算。

下面结合附图对本发明做进一步说明。

本发明在进行矿区dem变化检测时，提前根据矿区已有影像或者dem数据，确定并提取地表不变特征，然后根据不变特征的参数，建立不同时相间的dem配准模型，获取dem变化。

如图1所示，本基于地表不变特征的矿区dem变化检测方法具体包括以下步骤：

a.选取矿区所在位置，收集矿区地形地貌特点、地质采矿、开采、已有dem变化区域信息。

b.基于地表不变gcp特征的矿区dem数据建模；

首先针对矿区地形dem生成与变化检测需求，根据卫星观测时间获取对应矿区内ti和tj两个不同时相的立体影像数据，建立成像几何模型；

其次在分析矿区地形特点基础上，在有明显地标物的位置处布设地面控制点，为避免gcp受到开采扰动影响而出现位置变动，对gcp进行不变特征判断，判断方法如下：首先提取每个gcp在两期影像上的像方坐标；其次提取并计算gcp在两期dem上所对应的平面坐标、高程值、坡度和曲率参数；最后分别计算两期参数的差异，若差异在所设定的阈值范围之内，则该gcp为不变gcp；

最后根据成像几何模型和不变gcp，利用区域网平差技术对前视、正视、后视图像进行rpc模型精化，如图2所示，获取矿区大范围内基础地形dem。

c.利用数字地形分析中地形形态提取的方法和图像处理中的边缘检测方法提取步骤b中两个不同时相dem地形中的点、线、面特征基元。

d.对步骤c中的地形特征基元进行相似度测量，通过层次分析方法构建步骤c中提取的地形特征基元的几何特征测量模型，进行地形特征基元相似度测量，并依此建立地表不变特征判断规则，提取变化特征与不变特征，具体步骤为：

①按照相同的空间分辨率尺度对两期dem数据进行尺度一致重采样，利用最近邻点算法(iterativeclosestpoint，icp)实现两期dem地形特征基元的初始匹配；

②利用局部窗口搜索方法对范围内的地形特征基元的几何特征进行相似度测量，建立点特征相似度s1、线特征相似度s2和面特征相似度s3，如图3所示，点、线、面地形特征基元的相似度指标分量包括高程、水平位置、曲率、坡度、长度、弯曲度、面积、边界长度、边界形状等；

③计算同名特征的相似度s，具体方法包括：

假设特征l有n个比对特征分量项l＝[l1...ln]，则任意一项特征在两期dem上的相似度为

其中，li(i＝1,...,n)，li'(i＝1,...,n)分别为两期dem上对应的第i项特征，则该特征的相似度为

其中，pi为对应s(li)的权值；

④通过相似度最大值选取与阈值判断结合的方式实现对地形特征的匹配，根据相似度s计算结果范围为[0,1]，1代表没有变化，0代表最大变化，则设定地表不变特征的取值范围为[δs，1]，δs为相似度阈值，

e.根据步骤d确定的地表不变特征，构建不同时相dem间配准模型并进行平差解算，结合步骤d中提取的变化特征，获取dem变化，具体方法包括：

假设(x0,y0,z0)^t和(x,y,z)^t为同一地区待配准的两个dem不变特征上的点，采用区域网内平移旋转缩放模型，构建方程，

(x,y,z)^t＝t+sr(x0,y0,z0)^t

式中，s为缩放系数，r为旋转系数(ω，κ)，t为平移系数(tx,ty,tz)^t；

泰勒级数展开线性化后建立m个不变特征点的观测方程，

-e＝ax-l,p

式中，a为系数针，x为待解配准参数，p为权阵，l为残差向量。联合直线特征和面特征约束观测方程：

考虑稳健解算情况下，利用迭代最小二乘法联合解算式中的待解参数。

f.将配准后的dem与基准dem直接求差获取该时相的dem变化信息，该dem变化信息与步骤d中提取的变化特征一起，共同构成该矿区内dem变化信息。

本基于地表不变特征的矿区dem变化检测方法基于多时相影像dem在矿区的应用，在矿区地表扰动频繁的情况下，引入矿区地表不变特征，建立基于地表不变特征的dem配准与变化检测的技术方案，通过矿区点、线、面等地形特征基元的提取与相似度测量分析，更加准确地揭示矿区dem变化时空分布规律，为矿区多源地形数据融合提供新的思路与方法。