一种山区卫星影像自动配准到地理底图上的方法

摘要

本发明公开了一种山区卫星影像自动配准到地理底图上的方法，它涉及遥感图像预处理中的几何精纠正技术领域，其具体步骤为：(一)、山区卫星影像地性线提取；(二)、DEM（Digital？Elevation？Model,数字高程模型）地性线提取；(三)、同名地物点库构建；(四)、几何精纠正运算；(五)、纠正结果误差估计。它实现了计算机自动纠正山地区域卫星影像的技术，能在保证控制点数量和质量满足应用需求的前提下，解决仅使用DEM地性线纠正山区卫星影像时的人工采集控制点主观性强、耗时耗力的不足。

说明书

技术领域：

本发明涉及一种山区卫星影像自动配准到地理底图上的方法，属于 遥感图像预处理中的几何精纠正技术领域。

背景技术：

作为遥感数字图像预处理不可缺少的环节，几何精纠正，特别是将 遥感图像快速、精确配准到地理底图上，是实现遥感数据普及化、大众 化应用，以及几乎所有的专业化遥感应用的基础。中低分辨率的山区卫 星影像，因小的道路、水系与周围地物混合，使用传统方式进行配准存 在采集控制点控制点较难找、数量不足、匹配精度不高的问题。从相对 稳定的DEM中提取地性线的几何精纠正方法，可达到增加山地区域的控 制点数量和控制点精度的目的。但是，若利用DEM或DEM生成的立体图 形来人工寻找或采集控制点，如山脊、沟谷、山峰、凹地等，受人为主 观判断影响，误差较大，且不容易实现计算机批处理。因此，针对山区 遥感影像中常规控制点难找，而基于DEM地性线手工采集控制点费时又 费力等问题，发明一种将山区卫星遥感图像自动配准到地理底图上方法， 对于满足遥感影像在山地区域的绝大多数行业应用需求具有重要的实践 意义。

发明内容：

针对上述问题，本发明要解决的技术问题是提供一种山区卫星影像 自动配准到地理底图上的方法。

本发明的一种山区卫星影像自动配准到地理底图上的方法，其具体 步骤为：一、山区卫星影像地性线提取；二、DEM地性线提取；三、同名 地物点库构建；四、几何精纠正运算；五、纠正结果误差估计。

作为优选，所述的山区卫星影像地性线提取步骤为：

1、利用图像格式转换工具，将彩色的原始遥感图像转为灰度图像， 记为g；

2、读取图像g，使用Sobel边缘算子对图像进行水平和垂直方向的 滤波，然后求取模值，得到梯度模值图像，记为bgm；

3、读取图像g，设定半径为R的圆形结构元素，记为disk(R)，进行 形态学开操作、图像腐蚀、图像重建、形态学关操作、图像膨胀、图像 重建、求反，得到图像，记为g₂；

4、标记g₂中局部极大值，得到极大值图像，记为fgm，将原图中极 大值处像素值设为255；

5、设置参数为A的结构元素，记为ones(A,A)，分别对极大值图像fgm 进行关操作、腐蚀、开操作，标记出前景，得到前景图像，记为fgm₂；

6、将图像g₂使用合适的阈值转化为二值图像，记为BW；

7、计算BW的距离D，进行分水岭变换，可使用Matlab工具箱中 的工具watershed函数；

8、修改梯度模值图像bgm，使梯度模值图像在标记的前景对象和背 景对象中有最小值，得到梯度模值图像bgm₂，然后对bgm₂进行分水岭变 换，得到分割遥感图像后的图像数字矩阵，记为L；

9、将矩阵L，导出为灰度图像，附加坐标信息，转换为带有地理坐 标信息的图像文件，记为IMG；

10、将IMG栅格矢量化，得到遥感图像的地性线矢量文件，记为 IMG_R；

11、将提取的矢量结果IMG_R与遥感影像叠加显示，目视判读地性 线是否符合应用需求，若不符合则返回步骤3和步骤5进行调整。

作为优选，所述的DEM地性线提取步骤为：

1、将DEM中每个像素点的海拔数值的高低理解为该像素点得到降水 后的流向，统计出每个像素点上的水流方向；

2、计算出每个像素点的累计汇水量；

3、预定义阈值T，当某点上的汇水量大于T时，则该像素点标记为 沟谷线上的点；

4、将DEM进行反转运算，把海拔高的区域变为海拔低的区域，海拔 低的区域变为海拔高的地区，根据步骤1-步骤3提取出山脊线上的点；

5、将选出的特征区域矢量化、线状要素细化、去掉细小断线，得到 沟谷线、山脊线，即DEM地性线，记为DEM_R。

作为优选，所述的同名地物点库构建的步骤为：

1、建立IMG_R和DEM_R拓扑关系，提取两矢量线的结点和节点；

2、根据结点属性表和节点属性表，标识出节点属性表中为结点的节 点记录，遍历这些记录，将结点所在节点的邻节点号加入节点属性表， 再附上结点坐标(u_x,u_y)、邻节点的坐标(a_x,a_y)、(b_x,b_y),可得到用于同 名地物点匹配的几何特征点库；

3、对所有交叉点进行循环遍历，去除位于X轴或者Y轴上的结点， 根据同名地物具有几何形状相似、空间位置偏差仅限几个像元大小、几 何偏移具有同向性的特点，筛选符合条件的结点对；

4、将同名地物点的坐标提取出来，建立同名地物点库；

5、整理为地理信息系统软件(如ArcGIS软件)可识别的几何精纠 正控制点txt文件。

作为优选，所述步骤3筛选符合条件的结点对的具体步骤包括：

A、设A_m是遥感图像上待纠正的结点，B_n是DEM上的参考结点；A_m 有k个相邻节点，其中与A_m成夹角的邻节点为A_m-1、A_m+1；B_n有L个 相邻节点，其中与B_n成夹角的邻节点为B_n-1、B_n+1。结点A_m和B_n的欧 式几何空间距离Length_A_mB_n≤D_max，其中D_max的大小由遥感影 像的空间分辨率和几何粗纠正的偏移误差决定；

B、点B_n相对于A_m的方位角为λ，系统级校正后的遥感图像偏移地理 底图的角度为φ，则λ≤φ；

C、δ为同名地物点偏移的角度最大值。设结点A_m在 二维平面中的投影为(u_x,u_y)；a为A_m与邻点A_m-1构成的向量，其在二 维平面中的坐标表达式为a＝{(a_x-u_x),(a_y-u_y)}；b为A_m与邻点A_m-1构成 的向量,其在二维平面中的坐标表达式为a＝{(b_x-b_x),(b_y-b_y)},则：

${\theta }_{A_m}=arccos\left[\frac{(a_x-u_x)\times (b_x-u_x)+(a_y-u_y)\times (b_y-u_y)}{\sqrt{{(a_x-u_x)}^2+{(a_y-u_y)}^2\times \sqrt{{(b_x-u_x)}^2+{(b_y-u_y)}^2}}}\right]---\left(1\right)$

同理可计算出参考结点与相邻节点的夹角的值；

D、设α_A、β_A分别为向量与X轴、Y轴的夹角，与之匹配的 与X轴、Y轴的夹角为α_B、β_B，它们之间的角度偏差最大为ξ， 则α_A-α_B≤ξ、β_A-β_B≤ξ；同理，构成夹角的另一向量与X轴、Y 轴的夹角也必须满足以上条件。以α_A、β_A的计算为例，

${\alpha }_A=arccos\left[\frac{(a_x-u_x)}{\sqrt{{(a_x-u_x)}^2+(a_y-u_y)2}}\right]---\left(2\right)$

${\beta }_A=arccos\left[\frac{(a_y-u_y)}{\sqrt{{(a_x-u_x)}^2+(a_y-u_y)2}}\right]---\left(3\right)$

作为优选，所述的几何精纠正运算步骤为：

1、几何精纠正转换模型构建；将地理信息系统软件可识别的几何精 纠正控制点txt文件导入几何精纠正模块，用RubberSheeting的模型 建立待纠正卫星图像和地理底图栅格像元点之间的坐标对应关系，并对 卫星图像矩阵进行变换；

2、图像重采样；使用双向线性内插法确定校正后图像上每个像素点 在原卫星图像中的亮度值。

作为优选，所述的纠正结果误差估计步骤为：

1、建立M×M像元大小的矢量网格，然后随机选取n个网格，在这 些网格中随机选取1个或几个沟谷线控制点，测量其距离真实同名地物 点的平均距离，建立用于估计总体误差的误差样本x₁、x₂、x₃…x_n，当 n＜30时，样本量满足小样本检验的条件；

2、使用误差样本的均值衡量纠正后的卫星影像存在的几何总体误 差μ，其计算公式为：

$\bar{x}=\frac{x_1+x_2+x_3+...x_n}{n}---\left(4\right)$

上式置信水平为1-α，置信下限为置信上限为 $\bar{x}+t_{\alpha /2}\times s/\sqrt{n},$其中s为样本标准差，$s=\sqrt{\frac{{\Sigma }_{i-1}^n(x_i-\bar{x})}{n-1}};$

3、当几何总体误差控制在1-2个原卫星图像像元大小的范围之内时， 匹配结果被认为达到应用需求；反之则返回修改同名地物点匹配的参数 值。

本发明的有益效果为：它实现了计算机自动纠正山地区域卫星影像 的技术，能在保证控制点数量和质量满足应用需求的前提下，解决仅使 用DEM地性线纠正山区卫星影像时的人工采集控制点主观性强、耗时耗 力的不足。

附图说明：

为了易于说明，本发明由下述的具体实施及附图作以详细描述。

图1为本发明中一种山区卫星影像自动配准到地理底图上的方法流 程图；

图2为本发明中地性线结点与邻节点示意图；

图3为实施例中基于DEM提取的地性线示意图；

图4为实施例中基于卫星影像提取的地性线示意图；

图5为实施例中DEM沟谷线与卫星影像叠加显示图(自动配准前)；

图6为实施例中DEM沟谷线与卫星影像叠加显示图(自动配准后)。

具体实施方式：

为使本发明的目的、技术方案和优点更加清楚明了，下面通过附图 中示出的具体实施例来描述本发明。但是应该理解，这些描述只是示例 性的，而并非要限制本发明的范围。此外，在以下说明中，省略了对公 知结构和技术的描述，以避免不必要地混淆本发明的概念。

参看图1，本具体实施方式采用以下技术方案：其具体步骤为：一、 山区卫星影像地性线提取；二、DEM地性线提取；三、同名地物点库构建； 四、几何精纠正运算；五、纠正结果误差估计。

作为优选，所述的山区卫星影像地性线提取步骤为：

1、利用图像格式转换工具，将彩色的原始遥感图像转为灰度图像， 记为g；

2、读取图像g，使用Sobel边缘算子对图像进行水平和垂直方向的 滤波，然后求取模值，得到梯度模值图像，记为bgm；

3、读取图像g，设定半径为R的圆形结构元素，记为disk(R)，进行 形态学开操作、图像腐蚀、图像重建、形态学关操作、图像膨胀、图像 重建、求反，得到图像，记为g₂；

4、标记g₂中局部极大值，得到极大值图像，记为fgm，将原图中极 大值处像素值设为255；

5、设置参数为A的结构元素，记为ones(A,A)，分别对极大值图像fgm 进行关操作、腐蚀、开操作，标记出前景，得到前景图像，记为fgm₂；

6、将图像g₂使用合适的阈值转化为二值图像，记为BW；

7、计算BW的距离D，进行分水岭变换，可使用Matlab工具箱中 的工具watershed函数；

8、修改梯度模值图像bgm，使梯度模值图像在标记的前景对象和背 景对象中有最小值，得到梯度模值图像bgm₂，然后对bgm₂进行分水岭变 换，得到分割遥感图像后的图像数字矩阵，记为L；

9、将矩阵L，导出为灰度图像，附加坐标信息，转换为带有地理坐 标信息的图像文件，记为IMG；

10、将IMG栅格矢量化，得到遥感图像的地性线矢量文件，记为 IMG_R；

11、将提取的矢量结果IMG_R与遥感影像叠加显示，目视判读地性 线是否符合应用需求，若不符合则返回步骤3和步骤5进行调整。

上述步骤3和步骤5中合适的结构元素参数R或A值是提取遥感图 像地性线的关键，可通过多次试验或根据经验定义它们的大小。

进一步的，所述的DEM地性线提取步骤为：

1、将DEM中每个像素点的海拔数值的高低理解为该像素点得到降水 后的流向，统计出每个像素点上的水流方向；

2、计算出每个像素点的累计汇水量；

3、预定义阈值T，当某点上的汇水量大于T时，则该像素点标记为 沟谷线上的点；

4、将DEM进行反转运算，把海拔高的区域变为海拔低的区域，海拔 低的区域变为海拔高的地区，根据步骤1-步骤3提取出山脊线上的点；

5、将选出的特征区域矢量化、线状要素细化、去掉细小断线，得到 沟谷线、山脊线，即DEM地性线，记为DEM_R。

上述步骤4中阈值T是提取DEM地性线的关键，它与DEM的比例尺及 其栅格点大小有关，可通过多次试验或根据经验定义它们的大小。

进一步的，所述的同名地物点库构建的步骤为：

1、建立IMG_R和DEM_R拓扑关系，提取两矢量线的结点和节点；

2、根据结点属性表和节点属性表，标识出节点属性表中为结点的节 点记录，遍历这些记录，将结点所在节点的邻节点号加入节点属性表， 再附上结点坐标(u_x,u_y)、邻节点的坐标(a_x,a_y)、(b_x,b_y),可得到用于同 名地物点匹配的几何特征点库；

3、对所有交叉点进行循环遍历，去除位于X轴或者Y轴上的结点， 根据同名地物具有几何形状相似、空间位置偏差仅限几个像元大小、几 何偏移具有同向性的特点，筛选符合条件的结点对；

4、将同名地物点的坐标提取出来，建立同名地物点库；

5、整理为地理信息系统软件(如ArcGIS软件)可识别的几何精纠 正控制点txt文件。

进一步的，所述步骤3筛选符合条件的结点对的具体步骤包括：

如图2所示，设A_m是遥感图像上待纠正的结点，B_n是DEM上的参 考结点；A_m有k个相邻节点，其中与A_m成夹角的邻节点为A_m-1、A_m+1； B_n有L个相邻节点，其中与B_n成夹角的邻节点为B_n-1、B_n+1。结点A_m和B_n匹配为同名地物点的条件可表达为：

A、设A_m是遥感图像上待纠正的结点，B_n是DEM上的参考结点；A_m有k个相邻节点，其中与A_m成夹角的邻节点为A_m-1、A_m+1；B_n有L个 相邻节点，其中与B_n成夹角的邻节点为B_n-1、B_n+1。结点A_m和B_n的欧 式几何空间距离Length_A_mB_n≤D_max，其中D_max的大小由遥感影 像的空间分辨率和几何粗纠正的偏移误差决定；

B、点B_n相对于A_m的方位角为λ，系统级校正后的遥感图像偏移地理 底图的角度为φ，则λ≤φ；

C、δ为同名地物点偏移的角度最大值。设结点A_m在 二维平面中的投影为(u_x,u_y)；a为A_m与邻点A_m-1构成的向量，其在二 维平面中的坐标表达式为a＝{(a_x-u_x),(a_y-u_y)}；b为A_m与邻点A_m-1构成 的向量,其在二维平面中的坐标表达式为a＝{(b_x-b_x),(b_y-b_y)},则：

${\theta }_{A_m}=arccos\left[\frac{(a_x-u_x)\times (b_x-u_x)+(a_y-u_y)\times (b_y-u_y)}{\sqrt{{(a_x-u_x)}^2+{(a_y-u_y)}^2\times \sqrt{{(b_x-u_x)}^2+{(b_y-u_y)}^2}}}\right]---\left(1\right)$

同理可计算出参考结点与相邻节点的夹角的值；

D、设α_A、β_A分别为向量与X轴、Y轴的夹角，与之匹配的 与X轴、Y轴的夹角为α_B、β_B，它们之间的角度偏差最大为ξ， 则α_A-α_B≤ξ、β_A-β_B≤ξ；同理，构成夹角的另一向量与X轴、Y 轴的夹角也必须满足以上条件。以α_A、β_A的计算为例，

${\alpha }_A=arccos\left[\frac{(a_x-u_x)}{\sqrt{{(a_x-u_x)}^2+(a_y-u_y)2}}\right]---\left(2\right)$

${\beta }_A=arccos\left[\frac{(a_y-u_y)}{\sqrt{{(a_x-u_x)}^2+(a_y-u_y)2}}\right]---\left(3\right)$

上述步骤A-D的特征在于：D_max、λ、δ、ξ因几何精纠正的遥 感影像载荷参数和系统级几何纠正的精度不同而不同，若写为计算机可 批处理的程序，需根据经验或多次实验预定义。

进一步的，所述的几何精纠正运算步骤为：

1、几何精纠正转换模型构建；将地理信息系统软件可识别的几何精 纠正控制点txt文件导入几何精纠正模块，用RubberSheeting的模型 建立待纠正卫星图像和地理底图栅格像元点之间的坐标对应关系，并对 卫星图像矩阵进行变换；

2、图像重采样；使用双向线性内插法确定校正后图像上每个像素点 在原卫星图像中的亮度值。

进一步的，所述的纠正结果误差估计步骤为：

1、建立M×M像元大小的矢量网格，然后随机选取n个网格，在这 些网格中随机选取1个或几个沟谷线控制点，测量其距离真实同名地物 点的平均距离，建立用于估计总体误差的误差样本x₁、x₂、x₃…x_n，当 n＜30时，样本量满足小样本检验的条件；

2、使用误差样本的均值衡量纠正后的卫星影像存在的几何总体误 差μ，其计算公式为：

$\bar{x}=\frac{x_1+x_2+x_3+...x_n}{n}---\left(4\right)$

上式置信水平为1-α，置信下限为置信上限为 $\bar{x}+t_{\alpha /2}\times s/\sqrt{n},$其中s为样本标准差，$s=\sqrt{\frac{{\Sigma }_{i-1}^n(x_i-\bar{x})}{n-1}};$

3、当几何总体误差控制在1-2个原卫星图像像元大小的范围之内时， 匹配结果被认为达到应用需求；反之则返回修改同名地物点匹配的参数 值的准确度。

本具体实施方式提出了一种山区卫星影像自动配准到地理底图的方 法。该方法以山区DEM地性线和遥感图像地性线为空间参照数据，根据 水文分析算法提取DEM地性线；使用基于标记的分水岭算法提取遥感影 像地性线；根据地性线结点空间关系，计算机自动采集几何特征点(拐 点或交叉点)；依据一定半径的空间范围内同名地物点几何特征相似、空 间偏移同向性原则，计算机自动搜寻和匹配同名地物点。

实施例：

以2006年11月18日获取的Landsat5TM影像为待匹配的卫星影 像，其卫星条带号为130，轨道号为043，是系统级纠正处理后的数据产 品。以1:50000地形图上的等高线地图层生成的DEM为地理底图。

取结构元素的参数R＝2，A＝1，得到卫星图像的地性线IMG_R， 共有10731条遥感影像地性线，如图3所示。

将阈值T预定义为10，提取DEM地性线IMG_R，共提取80862条沟 谷线，如图4所示。同名地物点空间搜索的距离D_max＝500m，遥感 图像偏移地理底图的角度为φ＝[0,90]∪[270,360]，同名地物点偏移的允许 角度δ＝15°，地性线与X轴、Y轴夹角的角度允许偏差ξ＝15°。实验区 总面积约164km²，计算机自动从遥感影像中快速采集待用特征点7516个， DEM中快速采集到待用特征点83813个。经同名地物点自动匹配获得了可 用于几何精纠正的控制点对多达1322个。

利用ArcMAP软件中的Georeference模块进行纠正处理，得到配准 后的卫星图像，如图5-6所示。待纠正点和参考点之间总均方差为 0.0001m。

建立间隔为500m的格网，并对其进行编码。随机选取无云覆盖的27 个格网，在每个格网中设置一个误差检验样本点对。共得到27个误差样 本，如表1所示。样本的平均误差标准差s＝18.53。取置信度 为90％时，自由度为27的t分布的双侧分位数表得t_0.1/2(26)＝1.71，样本 的允许误差为60.8，误差样本的置信上限为49.31，置信下限为37.15。 该结果表明，总体误差μ＝43.23m，误差区间估计值为(37.15m，49.31m)。 平均几何误差不超过2个像元，以90％的可信度可判读，总体纠正误差在 1-2个像元大小之间。

表1几何精纠正误差样本列表(单位：m)

以上显示和描述了本发明的基本原理和主要特征和本发明的优点。 本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施 例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范 围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要 求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其 等效物界定。