沿海矿山生态环境修复中遥感影像自动配准方法的应用

项谦和，杜 娟，陈春雷，汪建光

(1.中国地质大学(武汉)，湖北 武汉 430074；2.中国煤炭地质总局浙江煤炭地质局，浙江 杭州 310021；3.嘉兴市规划设计研究院有限公司，浙江 嘉兴 314000；4.浙江省测绘科学技术研究院，浙江 杭州 310030)

沿海矿山生态环境修复是以计算机技术、光电技术、网络通信技术、空间科学、信息科学为基础，以3S技术为核心，对废弃矿地现场及周边进行实地踏勘、海洋地质测绘，收集矿山水文、工程、环境及地质等方面的相关资料，并在海洋地质环境条件和矿山边坡稳定性分析的基础上，对山体开挖、矿坑回填及生态恢复进行综合设计，并对全过程提出严格的环保和安全要求的工作。由于沿海矿山实地测量耗时耗力，且在我国海岸线蜿蜒曲折的东部海岸，尤其是基岩海岸分布区，地形陡峭较难以到达，给人工实地测绘与调查带来很大困难。随着卫星遥感技术的发展，利用3S技术对测绘困难区域进行海岸线遥感解译及提取成为可能，并逐渐发展成熟。

对同一场景使用相同或不同的传感器(成像设备)，在不同条件下(天候、照度、摄像位置和角度等)获取的两个或多个图像一般都会有差异。同一场景的多次成像的差别可以表现在:不同的分辨率、不同的灰度属性、不同的位置(平移和旋转)、不同的尺度、不同的非线性变形等。图像配准是将不同时间、不同传感器(成像设备)或不同条件下(天候、照度、摄像位置和角度等)获取的同一场景的两幅或多幅图像进行匹配、叠加或处理的过程，是图像处理领域的一个基础问题。简单来说，图像配准就是将同一场景的不同图像“对齐”或进行广义的匹配，主要目的为消除图像间存在的几何畸变[1]。

自动配准的方法可分为两大类:基于区域和基于特征的配准方法[2]。①基于区域的配准方法是将待配准图像中一块区域与参考图像中的相同尺寸的区域在统计学上进行比较，其相似度评测标准是从两块区域的标准化交叉相关系数中取最大值。也可以先通过变换将图像由时域变换到频域，再进行配准。对位移量比较大的图像，可以先校正图像的旋转，再建立两幅图像之何的映射关系。但如果图像中存在比较大的噪声和灰度差异时，这个交叉相关测量标准就变得不可靠。②基于特征的配准方法有比较高的稳健性。它有特征抽取和特征配准两个过程。该算法首先从两幅图像中提取灰度变化明显的点、线、区域、轮廓等特征形成特征集，然后在两幅图像对应的特征集中利用特征匹配算法选择存在对应关系的特征对。对于非特征像素点利用插值等方法推算出对应匹配关系，从而实现两幅图像之间的逐像素配准。基于特征配准方法的关键是如何提取出稳健的特征及如何匹配特征[3-4]。提取特征一般要考虑：①是否能够提取出可靠的特征；②提取出的特征是否能够找到可靠的同名点；③特征是否具有足够的信息来完成匹配。现有的算法中选取的特征主要有点、直线(线段)、轮廓和区域(面元)等图像的几何特征。

本文通过对沿海矿山进行海洋测绘与地质调查，利用实地勘测和遥感调查，结合专题图件岸线提取和历史资料进行分析综合，形成并建立沿海矿山生态环境修复调查成果矢量数据库，提出一种基于闭合区域特征的图像配准方法，用空间关系加以约束找到同名对象即同名区域，从而实现更稳健地找到相似区域，满足对象相似匹配的局部属性相似及全局几何关系上的一致性。

1 区域对象特征提取

1.1 光谱特征—灰度均值

由于目标区域是一个像素的集合，在作了影像辐射一致性预处理后，本文认为两幅影像同名区域之间的光谱特征具有一定相似性。这里采用灰度均值作为区域的光谱特征描述。

1.2 仿射不变矩及之间的距离

区域特征的仿射不变描述子有很多，如矩仿射不变量、形状矩阵等。由于不变矩[6-7]对两幅旋转角、比例缩放较大的影像较敏感，所以本文算法取前3个作为小面元区域特征的描述。定义如下

(1)

对于基准图像的M区域和待配准图像的N个区域，可用M×N维的不变矩距离矩阵Dmn={dij}来表示目标区域之间的相似性。Φr(i)表示基准图像中第i个目标区域的第k个不变矩，Φs(j)表示待配准图像中第j个目标区域的第k个不变矩，则基准图像第i个目标区域和待配准图像中第j个目标区域的不变矩距离为

(2)

很显然，距离越小，表明两个目标区域越相似。

1.3 形状特征

目标区域形状特征与组成目标对象的像素分布的空间位置有关，与像素的灰度值没有关系。上述特征的具体描述如下：

(1)面积：目标区域所包含的像素个数。

(3)圆度：它是度量区域形状常用的量。

(4)紧凑度：它是一个描述区域形状接近于圆形程度的量，可以用来衡量对象形状的规则程度，其定义为

(3)

(4)

(5)区域形状的一维描述：为了能够方便区域对象间形状的相似性评价，须有一个统一的描述方法，本文利用了中心—距离标记图方法，对影像的形状从矢量角度进行描述。中心—距离标记图的优点在于，它可以忽略被描述区域对象的具体尺度，是一种归一化的单纯形状的简单描述方法，同时它还有力地克服了角度—距离标记图法的缺陷(即当中心点落在区域对象的矢量多边形外或矢量多边形为非凸多边形、有内凹的现象时，同一个角度值可能对应多个多边形的边界点，从而不能生成正确的矢量标记图)。矢量标记图实际上就是一种矢量多边形的一维函数表达方法，本文采用的方法是首先将两幅图像中区域对象的多边形进行采样，获取数量相同的采样点，然后分别计算多边形中线点到采样点的距离，这样实际上就构成了一条多边形的形状描述曲线。假设有两个多边形A、B，它们外形相似但存在一定的角度旋转，利用中心—距离法对它们的形状进行表达，形状归一化曲线如图1所示。

图1 多边形及其形状曲线

如果令其中一个曲线不动，而对另一个曲线相对它进行点位移动，当移动到某点x0时，两曲线的吻合程度达到最佳吻合，也就是相似度最大，那么此最大相似度即为两多边形的形状相似度，点x0就是不动曲线起始点的同名点，它们与各自所在多边形中心点的连线相对水平方向的角度差即为两幅图像的旋转角度。多边形标记图的比较公式为

CurveShape(A,B)=

j=0,1,2,…,m

(5)

1.4 空间特征提取

对基准图像中的M个区域分别求取它们的重心点，记为Gr={(Xi,Yi),i=1,2,…,M}；对待配准图像中的N个区域分别求取重心点，记为Gs={(Xj,Yj),j=1,2,…,N}。在空间关系匹配中，每一个重心点就代表了一个区域对象。

2 特征匹配

2.1 属性特征匹配

本文将区域对象属性特征的匹配分成两个步骤：

第1步是根据对象的光谱特性、仿射不变矩、形状特征等因素综合决定，作一个同名区域的初始选择。因此可以把相似性定义为

Similarity(A,B)=Spectral(A,B)×Wspectral+

Invariment(A,B)×Winvariment+

Shape(A,B)×Wshape+

Compact(A,B)×Wcompact

(6)

交叉进行匹配：

对于CM中的每一对区域对象分别取其矢量多边形，计算它们形状曲线的最大相关系数MAX(CurveShape(Ai′,Bj′))，并计算相关系数最大时两个多边形之间角度的旋转差异值θ。若

(7)

需要注意的是，当区域目标为狭长形或矩形、正方形、圆形等对称形状时，在计算相关系数时极容易出现不止一个峰值的情况，此时获得的旋转角度可能与正式的旋转角度相差180°±σ、90°±σ等(σ为容差值)。由于这类目标通常形状曲线的相关系数较大，角度却与其他组相关系数很大的目标相差某一个常量，解决方法有两种：第1种是将此类目标单独列出，不参与特征匹配的第2步(即形状曲线的相关)，待其他组目标区域完成形状曲线的相关匹配后，与匹配目标一起进入空间关系约束环节；第2种是在统计角度直方图时根据已得到的角度峰值，判断是否需要加减一个常量。

当完成上述属性匹配后，筛选得到了一组按形状曲线相关系数从大到小排列的对象CM′={(ri,si),i=0,1,2,…,m-1}，ri表示基准图像中大区域目标，si表示待配准图像中的区域目标。这组对象在形状上已基本相似(若形状曲线的相关系数≥0.9，已经可以认定这是一对同名对象)，但是考虑一方面不同位置上的目标可能具有相似的形状，如城区中的建筑群；另一方面不同的地物在分割后可能具有相似的形状，因此本文从全局空间关系上再对这组目标作约束，以保证最后得到的匹配目标不仅具有局部属性上的相似性，还具有全局空间关系上的一致性。

2.2 空间关系约束

当前中高分辨率的遥感影像主要是Gauss-Kruger投影或UTM(universal transverse mercator)投影，它们均属于等角投影体系。在此投影体系下校正所得的遥感影像中，地物之间的角度基本保持不变，而长度和面积的变形较小。在本文中，利用这一性质对属性匹配后的结果进行最终的确认。

用已经获得的重心点代替目标区域本身，在基准图像和待配准图像中分别以各重心点构造三角网，判断两图中相应的目标点构建三角形的各个对应内角是否满足等角关系，即对应的三角形是否相似，从而对属性匹配的区域对象进行确认。在点集匹配方法中，生长匹配法是一种常用且较成熟的算法，本文就采用了这种方法对区域目标在空间关系上进行匹配。具体步骤如下：

(2)若m≥3，首先利用遥感地物的保角性质，从CM′中任取3对目标对象，分别在两幅对象中构建对应的三角形，计算三角形对应内角的角度差的和。如果该数值小于预设的阈值，那么就认为两图中的3个目标对组成的2个三角形近似相似，这3个目标对则被认为是基本匹配组。以目标点间的角度作为约束条件，对CM′中剩下的目标点进行匹配。以基本组的3个目标点对作为基础，从CM′中逐一挑选新的目标点对加入各自图像中，与3个基础点构成3对新的三角形，如果两图中这3对新加入的三角形满足相似条件，那么可认为新加入的目标点对相匹配，否则不能匹配。

3 变换模型求解

图像间的变换模型有多种，常用的有刚体变换、相似变换、仿射变换以及多项式变换等，本文根据最后得到的同名度区域的数量采用了两种变换方式。经过处理后得到的同名区域数小于3对，采用相似变换模型。比例误差可由同名对象面积比得到，得到了角度差异值和比例参数后，可以利用同名对象的质心位置计算影像间的平移误差参数。因而这时候对影像间的匹配实际上是一种粗匹配。假设变换前后的图像坐标分别为(x,y)和(X,Y)，则相似变换可表示为

(8)

式中，s为比例参数；θ为逆时针方向旋转角；Δx、Δy分别为水平、垂直方向的平移量。经过处理后得到的同名区域数大于等于3对时，可以采用多项式纠正模型。

4 试验结果与分析

对具有一定重叠区域的光学影像进行了试验，试验结果如图2—图4所示。

图2 对需要配准的两幅影像进行基于凸面模型的多尺度分割

图2是对两幅影像采用基于凸面模型的多尺度分割的结果，其中图2(a)、图2(c)为待配准影像及对其分割的结果，图2(b)、图2(d)为基准影像及对其分割的结果。图3(a)、图3(b)分别显示了在两幅影像中最后确认的同名目标区域，图4为采用仿射变换模型进行配准的结果与基准图的叠加显示。为了检查配准精度，在上述试验的影像上分别选取10对均匀分布的独立点对，计算其均方误差RMSE来评价。影像配准精度计算结果见表1。RMSE定义为

图3 显示寻找到的同名区域

表1 影像配准精度

图4 配准后结果叠加显示

(9)

式中，m为点数；xi、yi为待配准影像上点的坐标；Xi、Yi为对应参考影像上点的坐标。

5 结 语

从试验结果可看出，采用面向区域的原理进行影像配准，将分割后的目标区域作为一个对象，具有良好的稳定性且特征信息丰富。本文算法从目标对象的局部特征相似性和空间关系一致性出发[8]，快速有效地实现了遥感影像的自动配准。本文以基于凸面模型的多尺度分割为基础，对分割后的影像进行区域特征提取，结合了区域特征之间的自身相似性和它们在空间关系上的一致性原则，寻找同名区域，实现影像间的自动配准，取得了较好的结果。由于本文算法先对图像进行了分割，再利用分割图像寻找整幅图像中存在的一些明显区域，因此算法的可靠性在一定程度上依赖于分割算法的优劣。

本文从地物遥感光谱特征出发对遥感影像中的岸线进行判读解译，快速有效地实现了遥感影像的快速自动配准，有助于促进沿海矿山生态修复权属调查工作准确快速找准切入点，对人工实地测量调查难以企及的地区有较好的互补作用。沿海矿山生态修复调查工作中，应注意与宗地、海籍权属调查同步结合，更加清楚地把握范围线、岸线、土地所有权、土地权属、土地性质及海域资源使用现状，为更加准确地确定沿海矿山生态修复项目权属界址线位置提供了依据与保障。

通过生态环境修复方案的实施，利用综合矿区周边交通结构、自然地形和康养资源的现状，根据上述布局原则及康养产业空间发展战略，结合矿区自然生态优化和未来康养产业发展的需要，规划以矿区为核心、以矿区为特色化、以周边有机农业为基质、以海岸线或乡村民居为补充的空间差异战略，在此基础上逐步恢复和重建矿山生态系统，达到与周边环境的协调发展，实现高效综合利用。