一种结合已有DOM和DEM的高精度自动刺像控点方法

王明明，王 佳，冯仲科，申朝永

(1. 北京林业大学精准林业北京市重点实验室，北京 100083; 2. 贵州省第三测绘院，贵州 贵阳 550004)

摄影测量作业过程中，为了控制加密，提高作业精度，一般须在测区内布设控制点。外业通过GPS测量等手段实地获取控制点坐标，内业中在单张像片上标注出控制点位置，通过内外业可得到控制点在地面测量坐标系和像平面坐标系中的坐标。在地理国情普查、重点区域变化监测等方面，为了满足现势性需求，通过航摄更新相关地理信息产品。当已具有时相1的DOM和DEM数据时，本文提出了一种利用IPS软件和其辅助软件，为时相2的航摄像片刺入控制点的方法。由于IPS软件在影像匹配阶段，利用几何自动统计方法的智能过滤器，可使匹配点提取质量最高达0.1像素[1]，因此，其刺点精度远远优于人眼识别刺点。技术流程如图1所示。

本文选择的试验区域位于贵州省安顺市平坝区。时相2航摄面积约12 km2，航摄地面分辨率优于10 cm，共拍摄一个架次7条航带156张像片，平均航高2202 m。已有此测区先期DOM和DEM，测区内最高点海拔1504 m，最低点海拔1078 m，涉及建筑物、森林、农田等多种地物。

1 软件介绍

处理过程中涉及IPS软件和自主编写的IPS控制点自动采集辅助程序。

1.1 IPS软件系统

IPS是ICAROS公司推出的快速精确处理航空影像的软件系统，可用于各种大规模制图和工程应用中，包括灾害监测、农业和林业分析、管道制图和能源审计等。IPS自动化程度高，从采集到完成所需时间快，操作干预和错误少，支持各种广泛任务要求、传感器有效载荷和处理条件。

1.2 IPS控制点自动采集辅助程序

为了辅助IPS软件完成高精度自动刺像控点工作，基于Microsoft Visual Studio 2010开发平台，采用C#语言，开发了IPS控制点自动采集辅助程序。该程序主要包括4部分：①裁切已有DOM；②生成DOM中心点坐标；③生成相机文件；④根据像点读取控制点坐标。如图2所示。

2 基础数据准备

2.1 DOM裁切

已有时相1的DOM和DEM数据,如图3所示。在ArcGIS中读取DOM左下角和右上角点位坐标，根据时相2航摄单张像片覆盖实地面积分析DOM裁切尺寸，以裁切后DOM覆盖范围略大于单张像片覆盖范围为宜。如图4所示，虚线表示待裁切DOM边界，L(XL，YL)和R(XR，YR)分别表示DOM左下角点和右上角点，a和b表示在横向和纵向裁切DOM的长度。为了将DOM规则裁切，生成如图4实线表示的格网，先计算格网左下角的起始点QS坐标(XQS，YQS)，之后再计算完全覆盖DOM所需要格网的列数CN和行数RN，生成裁切依据文件。

试验数据中，可读取到(XL,YL)=(613 725,2 924 051)，(XR,YR)=(618 200,2 929 458)。航飞单张像片表示地面范围约为0.8 km2，为了使裁切后的DOM与像片能精确匹配，此次DOM裁切规则为1500 m×1500 m，即a=1500 m,b=1500 m。下面结合试验数据，说明裁切依据文件生成过程。

(1) 根据a和b的位数，将(XL，YL)对应数位的数字全部换为0，得到QS0(XQS0，YQS0)，如本文(XQS0，YQS0)=(610 000,2 920 000)。

(2) 在(XQS0，YQS0)上分别加a和b的倍数，直至得到不大于(XL，YL)的最大值，将其赋给(XQS，YQS)。此处采用循环实现:

for (XQS=XQS0;XQS<=XL;XQS=XQS+a)

XQS=XQS-a;

for (YQS=YQS0;YQS<=YL;YQS=YQS+b)

YQS=YQS-b;

得到(XQS，YQS)=(613 000,2 923 000)。

(3) DOM右上角的坐标与起始点坐标作差，除以横向和纵向的裁切长度，之后取整再加一。根据式

(1)

可计算出裁切格网的列数CN=4，行数RN=5。

基于以上方法，IPS控制点自动采集辅助程序可生成裁切DOM所需要的TXT依据文件，此文件内容是以影像名、左下角X坐标、左下角Y坐标、右上角X坐标、右上角Y坐标的顺序表示一个图形。如613.000-2 923.000 613 000 2 923 000 614 500 2 924 500，表示影像名为613.000-2 923.000的DOM其左下角坐标为(613 000,2 923 000)，右上角坐标为(614 500,2 924 500)。之后利用裁切依据文件在Inpho中的OrthoVista模块中裁切DOM影像，得到裁切后的DOM及其对应的TFW头文件。

2.2 DOM虚拟POS生成

裁切后的DOM以影像的形式加载到IPS中作匹配处理，故需要POS数据和相机文件，POS数据一般格式为(影像名,X,Y,Z,Omega,Phi,Kappa)。DOM完成裁切后，小幅DOM以影像左下角坐标X-Y方式命名。DOM的虚拟POS生成过程中，X、Y坐标选择裁切后DOM的中心点坐标，Z坐标选择时相2航飞高度。已有DOM是标准正射影像，因此其翻滚角(Omega)、俯仰角(Phi)、航向角(Kappa)均为零。

从上文中可知裁切后DOM的左下角坐标(X左，Y左)和右上角坐标(X右，Y右)，可得到DOM中心点坐标

(2)

以DOM613.000-2 923.000为例，则可知其POS为(613.000-2 923.000,613 750,2 923 750,2202,0,0,0)。同样可得到其他裁切后DOM的POS数据。

2.3 DOM虚拟相机生成

为了与时相2像片有更为相似的匹配条件，根据时相2相机文件生成DOM虚拟相机文件。在IPS辅助软件中，生成DOM虚拟相机文件需要输入焦距f、相对航高H、正射影像长度a′和宽度b′(像素表示)、正射影像实地长度a、宽度b等参数。以上参数中，焦距参考时相2相机文件；相对航高根据平均航高和测区平均海拔作差得到；a、b已知，a′、b′可从裁切后DOM文件的属性中得到。根据以上参数，结合式(3)可先求得DOM地面采样间隔GSD，再根据式(4)求得像元大小μ。试验数据采集使用的是飞思P65+相机，其参数见表1。

(3)

表1 飞思P65+相机参数

裁切后DOM中，a=b=1500 m,a′=b′=15 000,则根据式(3)可得GSD=0.1 m，参考表1中的焦距f=45.729 1 mm，DOM像元大小μ=0.005 431 mm。总结得出DOM虚拟相机参数见表2。

表2 DOM虚拟相机参数

3 自动刺入控制点

3.1 影像匹配

基础数据准备完成后，在IPS软件中建立工程。将裁切后的时相1 DOM与时相2像片均视为影像文件，在工程中加载影像文件、相机文件、POS数据、测区最高点和最低点海拔等。在相机文件夹下，存放两个相机文件，一个是时相2拍摄相机文件，另一个是裁切后DOM的虚拟相机文件，建成工程后，修改DOM对应的相机文件为其虚拟相机文件。进入GeoViewer视图窗口，以DOM和像片所覆盖范围有重叠即有关系为准，手动添加DOM与像片之间的连接关系，完成连接关系的添加，如图5所示。

之后在Match模块中开始同名点匹配，完成后所有匹配点均列在主窗口的Points列表中。这些点既有时相1 DOM与时相2像片之间的同名点，也有仅存在于时相2像片之间的同名点。完成匹配后，在工程中的DP文件夹下有后缀为GPF和IPF的两个文件，其中GPF为地面点文件，IPB为像点二进制文件，通过“Tools-Bin to Directory”可将IPB文件转换为IPF文件。GPF文件中记录了所有匹配点在物方坐标系中的坐标，每张影像(裁切后的DOM和时相2像片)都有对应的IPF文件，记载了与此影像有关的所有同名点个数、点名及各点在此像片平面坐标系中的坐标。

3.2 制作控制点

对于同时在DOM和时相2像片上的匹配点p，图6反映了其在DOM像方平面坐标系与物方平面坐标系中的位置。根据IPF文件可知其在DOM像平面坐标系中的坐标(xp,yp),在上文中，OrthoVista软件在完成裁切后，会自动生成对应文件的TFW头文件。TFW中存储了影像在X方向和Y方向的地面采样间隔GSD、旋转系数，以及影像左上角栅格的中心点位(X0,Y0)。由式(5)可得DOM物方平面坐标系中心坐标(Xc,Yc)

(5)

由点p的像方平面坐标、DOM物方平面坐标中心点，以及采样间隔GSD,可求出p点在物方平面坐标系中的坐标

(6)

根据以上原理，在IPS辅助软件的根据像点读取控制点坐标模块中，输入GPF文件、IPF文件和TFW文件等，计算所有DOM与时相2像片匹配出的点的物方平面坐标(X,Y)，根据物方平面坐标和DEM可提出Z坐标，将(X,Y,Z)坐标一并写入GPF文件中。用新生成的GPF文件替换原有GPF文件，打开工程后可见DOM与时相2匹配点的物方坐标已被改写，而时相2像片间匹配的点坐标依然为(0,0,0)。

改变已写入物方坐标的点属性(从“tie”改为“xyz_check”)，之后通过平差，系统会给出每个点的残差大小，以残差大小为依据，结合目视检查的方法删掉有明显错误的“xyz_check”点。经过若干次平差后，删除“xyz_check”点以外的其他所有点和工程中的DOM影像。IPS软件中GPF文件和IPB(IPF)文件通过点名字段连接，通过上面的步骤可知，剩余点在GPF文件中存有其物方坐标，在IPB(IPF)文件中存有对应的像平面坐标，将“xyz_check”点属性改为“xyz_control”即完成了区域内刺像控点的工作。

3.3 后续处理

引入控制点后完成空三处理，平差之后的成果可以在IPS软件中进一步处理得到时相2的地信产品，也可直接将其导出为PATB格式，或分析其数据格式，将其转换为Inpho或其他后续软件文件格式，以期在其他软件中完成地信产品生产。在示例数据中，平差之后得到RMSPixel=1.421,RMSX_Control=0.224,RMSY_Control=0.191,RMSZ_Control=0.133。依次执行DTMe—OrthoRectifyMe—StitchMe—Seam line Editor等步骤，完成了时相2 DOM影像的生成。通过对完成后的DOM影像的分析，其结果完全符合1∶1000国家规范。

4 结 语

本文阐述了在地理信息产品更新过程中，采用计算机匹配的方法找出时相1 DOM与时相2像片之间的同名点，再通过DOM的本身属性和其对应的TFW文件，以及DEM对应点的高程数据，共同得到匹配点的物方坐标。通过平差探测与目视检查等方法排除了错误点位后，将这些点作为控制点存储，即完成了像控点高精度刺入的目的。以此为基础，进行下一步的匹配、解算、地理信息产品的生成等。

无人机数据的空间和时间分辨率远优于遥感数据，因此在需要以无人机影像作为遥感数据分析依据时，DOM同样可用高分遥感影像如WorldView、Google Earth等代替，处理所得结果可以与目标遥感影像高度匹配。在此处理中，采用计算机自动匹配的方法得到的点位精度远高于人工刺入点位，可以达到亚像元级别，同时得到的控制点更为均匀稠密地分布在测区内，为精确的空三处理提供基础保证。