MapGIS 影像镶嵌融合控制点获取的一种有效方法

赵卫东

（河北省地质矿产勘查开发局第一地质大队，河北 邯郸）

1 镶嵌融合控制点及其gcp 文件

Mapgis 中影像的镶嵌融合功能非常强大[1]，只要有足够的控制点或其文件（*.gcp）[2]，可以随心所欲的将msi 影像进行平移、旋转、缩放，甚至局部伸缩、揉皱。Mapgis 影像文件（.msi 文件）镶嵌融合控制点文件[3]，可用记事本打开，其格式见图1。

图1 gcp 文件格式

该文件第一行为表头，各列间以空格与制表符两个字符分隔，下面为控制点信息，以英文豆点与制表符两个字符分隔。经试验，表头及下面控制点信息各列之间有英文豆点分隔即可。gcp 文件最关键的信息即是控制点的图像坐标、理论坐标，相关的坐标系称为像素坐标系、图形坐标系，其意义如下。

（1） 图像坐标，在镶嵌融合控制点采集窗口左窗口状态栏显示称为“图像坐标”，即msi 影像像素的坐标，以像素为单位。像素坐标系以影像左上角为像素坐标系的坐标原点，向右、向下为正向。msi 文件的像素信息与msi 影像转换前的原光栅图像相同，将鼠标悬停于光栅图像可出现像素列数、行数等信息。

（2） 理论坐标，为目标图形坐标。图形坐标是在镶嵌融合控制点采集窗口中右侧窗口的Mapgis 的图形坐标，亦即mapis 点、线、面文件编辑窗口显示的坐标，以mm 为单位。图形坐标系横向x 轴向右为正，纵向y 轴向上为正。当地图图面比例尺为1:1 000，且没有旋转、平移，则图形坐标（以mm 为单位）的数值即为公里网坐标数值（以m 为单位）。

所谓控制点，即是同一点具备这两种坐标，将影像的像素坐标进行校正与理论坐标对应。

以上，明确了像素坐标系与图形坐标系，以及控制点文件*.gcp 格式，现在需要获取控制点信息。

2 获取控制点

设一光栅地图其四角为经纬度坐标，图中有已知坐标的公里网格，本文研究推荐的获取控制点步骤如下。

2.1 矢量化控制点

将光栅图像初转换的msi 文件添加至Mapgis 点、线、面文件编辑窗口，用合适的子图（如十字丝形）标记控制点：第1，标记msi 影像右上角点；第2，标记右下角点；第3，从左下按逆时针顺序标记图框四角经纬度角点（此顺序只是建议）；第4，标记足够的公里网格交点，可随机性标记。

2.2 控制点赋属性与导出

第5，利用Mapgis 二次开发程序section 中“1 辅助工具”→“导入导出功能”→“点位置转属性”，让每个控制点子图具备图形坐标属性；第6，增加控制点文件“像素x”、“像素y”、“经度”、“纬度”、“理论Y”、“理论X”、“备注”等属性；第7，将msi 影像的右上、右下角点根据原光栅图像的像素行列信息赋予像素坐标，将“第3”步的四角赋予经度纬度坐标属性。

第8，section 中“1 辅助工具”→“导入导出功能”→“导出属性数据(EXCEL)”，将刚标记的控制点信息导出至excel。

2.3 图框四角经纬度坐标正算为公里网格坐标

在导出的Excel 表格中将图框四角经纬度换算为公里网坐标，注意换算时与原光栅地图投影参数一致。这步需要经纬度与公里网坐标正反算的软件，如Mapgis、二次开发的Excel、CoordTools_7.0.0.exe、arcgis等，资源丰富。

2.4 计算所有控制点的像素坐标

设msi 影像有m 列、n 行像素，初始未校正的msi影像右上角P 的图形坐标为（xP，yP），E 为任一点，其图形坐标为(xE，yE），像素坐标（xF，yF)，显然，初转换的msi 影像像素坐标系相对图形坐标系只有平移、缩放，没有旋转，则：

利用式（1）可计算各控制点的像素坐标系坐标。

2.5 计算各控制点的理论坐标

控制点的理论坐标即是将原msi 影像的控制点的图形坐标，计算至目标图形坐标（理论坐标）。根据任两个已知理论公里网坐标的控制点，计算各控制点的公里网坐标，然后按公里网格间距取整，即得各控制点的理论公里网坐标。图框四角点的理论公里网坐标，即2.3 节正算的结果，不能取整，若有其它特殊的控制点，也需要像图框四角一样，单独赋予其理论公里网坐标。

总之，本研究正是基于地图图面上有规律的公里网或经纬网，利用其规律获取其交点的理论公里网坐标或经纬度坐标，根据需要的比例尺缩放至理论坐标（目标图形坐标），与计算而得的各控制点的图像坐标，形成控制点gcp 文件，用于影像校正。

前面2.1～2.4 节的步骤均易实现，本小节问题较为复杂，这是本研究要解决的关键问题之一。抽象化后此问题即：平面坐标系有一组点形成“刚体”，其中两点对齐这组点外另两点，求这组点中其它点对齐后的坐标，这种变换称为平面两点对齐变换，下面阐述。

3 平面两点对齐变换

设A（xA，yA）、B（xB，yB）、C（xC，yC）、D（xD，yD）为平面坐标系xoy 中不重合的四点，E（xE，yE）为另外任意一点，现令A、C、E 三点组成“刚体”，做变换：A 点对齐B 点、C 点对齐D 点，则E 点将对齐F（xF，yF）点，现在已知A、B、C、D、E 五点的坐标，求F 点坐标。引用矢量的内积、外积定义与性质，及矢量旋转矩阵[4]，经推导，有

式中：

式（2）即可完成平面二点对齐变换，其中ABCDE五点的坐标均为已知，可求得F 点坐标。

虽然，式（2）在Excel 中编辑公式可以实现，但若应用Excel 内置的VBA[5]编辑自定义函数，更加方便快捷。

4 计算控制点的理论公里网格坐标

利用式（2），可求得对应F 点坐标（xF，yF），即该控制点的公里网坐标，但不是最终理论公里网坐标，需要将其按公里网格间距取整。理论公里网坐标一般为整数，甚至是10、100、1 000、…的倍数，比如，公里网格间距为500 m，按100 m 取整，即可获得正确的理论公里网坐标。

各控制点的xF、yF 坐标还可以继续变换，如缩放（可控制msi 影像校正后的比例尺）、换带计算、不同坐标系转换（简单的平移）等，这些在Excel 中容易实现。这样事情做好一次后，在Excel 中可以复制，极方便形成gcp 校正文件，这才是本研究倡导的这种获取gcp文件方法的方便快捷之处。

5 实例gcp 文件形成

图2 为获取msi 影像校正gcp 文件的实例截图（图中设置了小数位数显示，实际小数未显示全）。

图2 计算像素坐标系坐标、理论公里网坐标实例

该msi 影像有m=12 153 列、n=9 697 行像素，msi影像右上角的图形坐标（12 153，9 697），第1 行为表头，注意：大写的X、Y 表示南北向、东西向坐标，后面有114、117 表示该公里网坐标的投影中央子午线经度，小写的x、y 表示东西向、南北坐标，经度、纬度以“d.ffmmm”格式表示，如114.123 456 表示114°12′34.56″。第2 行为要校正的地图名称。第3、4 行为msi 影像右上、右下角点坐标信息。第5～8 行为地图图框四角点坐标信息，已知其经纬度坐标，公里网坐标为正算获得。第9 行及之下为公里网交点的控制点，图中公式栏中显示单元格H9:I9 的数组公式编辑状态，当“ctrl+shift+enter”后，公式编辑栏显示见图3。

图3 公式编辑栏中显示的数组公式

表明自定义函数xFyF_from_AtoB_CtoD_EtoF 存在于工作簿《Hong2024.xlsm》中。在其它工作簿中引用该自定义函数时，必须同时将《Hong2024.xlsm》打开，否则找不到该函数。

计算像素坐标所引用的公式显示见图4、图5。

图4 D5 单元计算像素x 坐标

图5 E5 单元格计算像素y 坐标

上面已经获取了形成gcp 文件的要素，要形成gcp 文件，只需将这些要素提取至gcp 文件中。可在图2 中的后面列如PQRS 四列中直接重现前面像素坐标与理论公里网坐标的四列内容（令该单元“=”另外的单元格即可），在O 列形成以1 起始的序号列，将这五列复制到已有的gcp 文件中，代替gcp 文件中原有内容，再将这五列分隔符（tab 制表符）替换为英文豆点，保存。也可编辑公式直接形成这五列内容，代替某gcp文件中原有内容，形成新的gcp 文件。

6 结论

本研究倡导的这种获取gcp 文件的方法有如下特点：

（1） 矢量化控制点比传统方法简化了很多切换窗口、缩放窗口、平移窗口、点击鼠标的步骤，也不必制作参照文件，控制点采集不考虑顺序，能有效避免人工错误，大大改善了控制点采集体验，当有大量的msi 文件需要校正时，更能提高效率。

（2） 形成一次获取gcp 文件工作表后，注意公式中对单元格的相对引用、绝对引用，在Excel 中可以复制整个实例，只粘贴替换图形坐标两列与msi 影像的像素列数、行数，处理好特殊情况如图框四角点的经纬度坐标，即可自动计算，极方便形成gcp 校正文件，这是本研究倡导的这种获取gcp 文件方法的方便快捷之处。

（3） 理论公里网坐标还可以再次转换，如：将1954 北京坐标系坐标转换为CGCS2000 坐标，或者六度带换算为三度带等，获得最终的理论公里网坐标，可轻易实现msi 影像的坐标变换、比例缩放。

（4） 若将图形坐标与理论坐标组成坐标对，可形成误差校正文件，用于Mpagis 图形文件的误差校正。