一种线阵推扫式行星遥感影像摄影测量快速处理方法

耿 迅，徐 青

（1.河南大学 地理与环境学院，开封 475004；2.河南大学 河南省时空大数据产业技术研究院，郑州 450000；3.信息工程大学 地理空间信息学院，郑州 450052）

引 言

地外天体表面高精度的正射影像图（Digital Orthophoto Map，DOM）与数字高程模型（Digital Elevation Model，DEM）是深空探测任务实施与行星科学研究的重要基础数据[1-4]。使用摄影测量方法对地外天体遥感影像进行几何处理是制作行星制图产品的主要技术手段[5-7]。现有行星遥感制图产品的精度与分辨率仍然较低，不能很好地满足深空探测工程任务以及深层次行星科学研究的需要。当前，利用不同航天机构获取的多探测任务数据制作高分辨率、高精度的行星制图产品是深空探测遥感测绘的研究重点，但是融合多探测任务数据进行摄影测量处理的技术难度很大，而且行星摄影测量处理过程复杂，可用软件资源有限。

结合地外天体遥感影像的特点，有针对性地设计、优化行星摄影测量处理方法，是深空探测任务科学数据处理与应用的重要工作。深空探测器通常搭载线阵传感器获取地外天体的遥感影像用于制作行星制图产品[8-11]。然而，现有行星摄影测量处理方法在用于线阵推式遥感影像时的效果不佳，通常表现出较低的处理效率[12]，美国地质调查局（United States Geological Survey，USGS）开发的ISIS（Integrated System for Imagers and Spectrometers）软件在处理生成大数据量（例如GB级）线阵推扫式行星遥感影像的正射影像图时，通常需耗时几十分钟甚至若干小时[13]。另外，在构建连接点控制网时，现有行星摄影测量处理方法针对线阵推扫式遥感影像的适用性有限，经常出现匹配失效的情况，导致后序摄影测量处理步骤以及行星制图产品精度受到影响[14]。针对现有行星摄影测量处理方法在线阵推扫式遥感影像上的不足，本文采用快速反投影算法进行正射影像纠正，并且在近似正射影像空间匹配获取连接点，提升长条带、大数据量线阵推扫式行星遥感影像的摄影测量处理效率。

1 行星摄影测量处理技术

经过几十年的积累，美国在深空探测遥感测绘领域已经形成了一套比较成熟的行星摄影测量处理技术体系，可以将这个技术体系概括为“PDS+SPICE+ISIS”，其中PDS（Planetary Data System）是行星遥感数据系统，SPICE是深空探测任务辅助信息系统，ISIS是开源行星遥感影像处理软件。另外，近年来美国国家航空航天局（National Aeronautics and Space Administration，NASA）研发的ASP（Ames Stereo Pipeline）开源软件也逐渐在行星制图领域中得到应用。

PDS是NASA研发的行星遥感数据系统，目前作为深空探测任务科学数据的标准文件格式被世界各国广泛采用。中国“嫦娥”工程与“天问一号”火星探测工程的科学数据也采用了PDS文件格式。PDS文件除了记录影像信息，也包含有影像获取时间、分辨率、地理坐标范围等辅助信息。可以使用开源图像处理引擎GDAL对PDS文件进行读写、格式转换等操作。

行星遥感影像几何处理涉及到复杂的行星坐标系转换、各种时间系统转换、探测器位置与姿态计算、相机模型构建、光束法平差等内容。NASA使用SPICE库文件对这些深空探测任务辅助数据进行归档[15]。SPICE库文件也是行星摄影测量处理时重要的辅助数据，其中SPK文件对应于探测器位置数据，CK文件对应于探测器姿态数据，利用这两个文件可以生成影像的外方位元素信息，用于构建相机模型以及光束法平差处理。另外，为便于利用这些SPICE库文件进行各种计算，SPICE辅助信息系统也提供了相应的库函数，可以使用C语言进行调用。

美国自实施“阿波罗”登月计划以来，一直十分重视深空探测遥感测绘工作，基于USGS ISIS软件制作了月球、火星、小行星等多个星体的制图产品，在深空探测着陆区选址、任务路径规划以及行星科学研究等方面发挥了重要作用。USGS开发的ISIS行星摄影测量软件系统能够支持NASA的多个深空探测任务数据。ISIS提供了很多工具软件用于行星遥感影像的摄影测量处理，其它航天机构，如欧洲航天局（European Space Agency，ESA）、日本宇宙航空研究开发机构（Japan Aerospace Exploration Agency，JAXA）也使用ISIS进行深空探测科学数据的预处理，ISIS在深空探测遥感测绘领域得到了行业的广泛认可。

ASP是由NASA的艾姆斯研究中心（Ames Research Center）开发的立体摄影测量软件，可以支持行星遥感影像与对地观测遥感影像[16]。ASP软件在生成DEM方面的效果较好，并且自动化程度较高，提供有stereo、point2dem等多个工具软件用于自动生成DEM，stereo软件主要是匹配立体影像生成点云文件，point2dem可以将点云文件转换为常用格式的DEM文件（如GeoTiff格式），并且生成DEM交会残差图辅助分析DEM成果精度。ASP支持ISIS的cube文件格式，能够读取其相机模型信息。

但是与商业摄影测量软件相比，现有开源行星摄影测量软件在交互式编辑、线阵推扫式遥感影像处理效率等方面的功能稍弱。然而，商业摄影测量软件通常无法直接用于处理月球、火星等地外天体遥感影像，主要原因是商业摄影测量软件无法支持行星遥感影像的数据格式、相机模型、坐标基准等；另外，商业摄影测量软件的核心算法（例如影像匹配）也通常结合对地观测卫星影像的特点进行了优化，而这些优化参数并不适合于处理行星遥感影像。

2 基于快速几何纠正的行星摄影测量处理方法

本文基于线阵推扫式遥感影像快速反投影算法进行行星遥感影像的几何纠正与控制网构建，并自主研发了相应的软件模块，应用于地外天体的摄影测量处理过程。基于快速几何纠正的线阵行星遥感影像摄影测量处理流程如图1所示，具体处理步骤可以分为影像预处理、相机模型构建、连接点提取、光束法平差以及制图产品生成等几个步骤。经过几何处理生成的正射影像图与数字高程模型可以加载至ArcGIS、GlobalMapper、QGIS等地理信息系统（Geographic Information System，GIS）软件中进行科学分析。由于USGS ISIS支持多个深空探测任务的科学数据，而且在行星测绘领域得到广泛应用，尤其是ISIS内置了许多深空探测任务的相机参数，因此本文的行星摄影测量处理流程也使用USGS ISIS软件进行格式转换、辐射校正以及光束法平差等处理。

图1 基于快速几何纠正的行星遥感影像摄影测量处理流程图Fig.1 Flowchart of the photogrammetric processing of planetary remote sensing images based on fast geometric rectification

在影像预处理阶段需要导入PDS格式的原始影像并转换为USGS ISIS系统支持的cube文件格式。ISIS针对不同深空探测任务的传感器开发了专门的数据导入接口，例如lronac2isis用于导入月球侦察轨道器（Lunar Reconnaissance Orbiter，LRO）窄角相机（Narrow Angle Camera，NAC）影像，hrsc2isis用于导入火星快车（Mars Express，MEX）高分辨率立体相机（Hight Resolution Stereo Camera，HRSC）影像。在预处理阶段还需要对影像进行辐射检校，例如使用lronaccal与lronacecho对LRO NAC影像进行辐射检校。在预处理阶段的一个重要步骤是调用spiceinit获取与影像相关的各个SPICE库文件，从而计算出影像在摄影时刻的位置、姿态信息，这一步相当于构建了相机模型，以便在像点与地面点之间进行坐标转换。

2.1 线阵推扫式影像地面点反投影

构建遥感影像的相机模型是开展摄影测量处理的基础[17-18]。相机模型主要用于地面点坐标P(X,Y,Z)与像点坐标p(x,y)的相互转换。以MEX HRSC影像为例，地面点P(X,Y,Z)与像点p(x,y)的坐标转换公式为

图2 线阵推扫式影像地面点反投影示意图Fig.2 Illustration of back projection of linear pushbroom images

为了提升最佳扫描线的迭代计算效率，本文采用基于物方投影面约束的快速反投影算法进行线阵推扫式影像的地面点反投影计算[18-19]。先利用线阵推扫式影像的严密几何模型建立每一条扫描线的物方投影面（例如第i条扫描线的物方投影面为SiAiBi），然后在地面点反投影迭代计算时，利用地面点至投影面的距离作为几何约束条件辅助确定地面点对应的最佳扫描线。地面点P与某一扫描线所构成的物方投影面之间的距离 L的计算公式为

其中：A、B、C、D是物方投影面的法线方向，可由空间3点Si、Ai、Bi的三维坐标计算得出。由于基于物方投影面的迭代计算过程是简单的几何计算，并且迭代收敛速度快，因此地面点反投影过程的计算效率可以得到显著提升。

2.2 快速几何纠正

线阵推扫式遥感影像的几何纠正通常采用反解法，即对于正射影像图上的每一个像素，先计算出其对应的地面点坐标，再由地面点坐标反求出像点坐标（即地面点反投影过程），然后将该像点的灰度值进行内插处理并赋值给正射影像图。USGS ISIS系统提供了cam2map工具软件用于生成DOM。由于ISIS系统的线阵推扫式相机模型采用计算效率偏低的二分法确定最佳扫描线，导致地面点反投影计算的效率较低，因此cam2map在处理线阵推扫式行星遥感影像（例如LRO NAC、MEX HRSC）时非常耗时。

本文通过优化线阵推扫式遥感影像的相机模型提升地面点反投影迭代计算效率，并且采用多线程程序设计方法进一步提高正射影像生成的处理效率，与USGS ISIS软件相比，单线程几何纠正处理的计算效率可提升3～5倍，多线程几何纠正处理的计算效率可提升10～20倍（4线程）。

2.3 连接点控制网构建

通过spiceinit获取的影像初始位置、姿态信息（即外方位元素）的精度偏低，需要利用光束法平差进行精化。光束法平差的关键是获取均匀分布的可靠连接点，在USGS ISIS系统中称为控制网（control network）。高精度的连接点控制网是获得高质量行星摄影测量成果的前提。USGS ISIS提供了autoseed、cnetref、pointreg等一系列控制网构建工具，但是其计算效率较低，并且经常出现匹配处理失败的情况，不适合大数据量的连接点提取，而且相关匹配算法难以用于畸变较为严重的立体影像（例如MEX HRSC的前、后视图像）。

本文先利用线阵推扫式遥感影像的快速反投影算法生成近似正射影像，然后在近似正射影像上采用归一化相关系数方法匹配连接点，再利用地面点反投影算法将匹配的连接点转换至原始影像空间。由于在近似正射影像空间进行匹配，立体影像的分辨率变为一致，同时也消除了地形起伏与摄影角度不同带来的几何畸变影响，并且能够利用影像初始位置、姿态信息辅助确定匹配点的起始位置，因此能够提升连接点提取步骤的计算效率与匹配精度。

2.4 光束法平差

USGS ISIS提供了jigsaw工具用于光束法平差，jigsaw软件的计算效率较高，对深空探测任务各种传感器的支持也比较好。因此，本文将匹配的连接点转换至ISIS软件支持的PVL （Parameter Value Language）控制网文件格式进行光束法平差处理。光束法平差的一个难点是设置不同类型观测值的权值，主要涉及探测器位置与姿态观测值、控制点观测值等。权值设置可以依据探测器定轨、定姿精度进行设置，也需要一定的工程经验。光束法平差通常需要经过几次迭代，每一次平差时需要剔除连接点中的粗差，通常认为当像点观测值残差小于1个像素，平差指标Sigma0值小于0.5时可以认为获得了理想的平差结果（实际精度情况受探测数据质量影响）。若要提升绝对几何定位精度，则需要在光束法平差过程中引入控制数据。但是，现有行星表面控制数据的分辨率与精度较低，尤其是缺乏高精度的平面控制数据，实际工程应用中通常引入激光测高数据或者由激光测高数据生成的DEM作为高程控制信息。

2.5 数字高程模型与正射影像图生成

在光束法平差后，利用立体像对进行影像匹配可以获取密集的匹配点，再利用立体影像空间前方交会计算出点云的三维坐标，从而自动构建出DEM。通常使用ISIS软件进行行星遥感影像的预处理与光束法平差，再利用ASP软件生成DEM。ASP软件自动生成DEM的分辨率与原始影像分辨率基本相同，但是存在较多噪声，实际应用时通常需要将自动生成的DEM重采样为3～5倍原始影像分辨率。

通过几何纠正可以消除各种畸变的影响，生成含有精确地理信息的DOM。生成DOM时需要利用影像的外方位元素以及参考DEM，因此光束法平差的精度以及DEM的精度会直接影响DOM的成果精度。在生成正射影像图时，可以利用USGS ISIS提供的cam2map软件以及自主研发的快速几何纠正软件。

3 试验与分析

试验选取4幅LRO NAC影像与5幅MEX HRSC影像进行处理（影像基本信息见表1）。LRO NAC影像数据位于“嫦娥四号”着陆区，影像分辨率为0.88～1.18 m；MEX HRSC影像数据覆盖Mars 2020任务预选着陆区的Jezero撞击坑，影像分辨率为14～22 m。利用USGS ISIS 3.5.2软件进行格式转换、辐射校正、光束法平差等处理，使用自主研发软件进行相机模型构建、几何纠正、连接点提取等处理，使用ASP 2.6.0软件生成数字高程模型。试验软件环境为64位Ubuntu 14.04操作系统（安装在VMWare虚拟机上）；硬件环境为Intel Core i7-7 500 2.7 GHz CPU与8 GB内存。月球与火星制图分别采用1 737.4 km与3 396.19 km的正球体，地图投影采用Equirectangular投影方式。

表1 测试影像基本信息Table 1 Basic information of the test images

3.1 光束法平差精度分析

基于初始的影像位置、姿态信息以及USGS ISIS软件内置的月球、火星全球DEM数据，使用自主研发的几何纠正软件模块将试验影像处理生成近似正射影像图，在近似正射影像图上提取连接点，再转换为ISIS系统支持的PVL控制网文件格式用于光束法平差。MEX HRSC影像共提取495个连接点，LRO NAC测试影像共获取187个连接点，两组实验数据的连接点分布情况见图3～4。

图3 MEX HRSC影像连接点分布示意图Fig.3 Illustration of tie points distribution of MEX HRSC images

图4 LRO NAC影像连接点分布示意图Fig.4 Illustration of tie points distribution of LRO NAC images

使用jigsaw进行光束法平差，两组实验数据的Simga0值均小于0.3，光束法平差后影像像方残差值见表2，残差分布见图5～6。分析像方残差结果可知，LRO NAC与MEX HRSC影像光束法平差的像方残差值均小于半个像素，表明影像匹配精度较高，不同条带之间外方位元素的不一致性得到消除。连接点几何定位精度见图7～8。分析结果可知，LRO NAC影像平面定位精度为0.5～1.0 m，高程定位精度为1.0～1.5 m；MEX HRSC影像平面定位精度为4～10 m，高程方向上大部分连接点的定位精度为10～20 m，少量连接点的定位精度接近30 m。结合影像分辨率可知，LRO NAC与MEX HRSC影像均可达到平面约0.5～1个像素，高程约1～2个像素的定位精度，与对地观测卫星遥感影像的几何定位精度基本一致。

图5 LRO NAC影像光束法平差后像方残差分布图Fig.5 Image residuals of the bundle adjustment for LRO NAC images

图6 MEX HRSC影像影像光束法平差后像方残差分布图Fig.6 Image residuals of the bundle adjustment for MEX HRSC images

图7 LRO NAC影像光束法平差连接点几何定位精度Fig.7 The posterior geometric accuracy of tie points for LRO NAC images

图8 MEX HRSC影像光束法平差连接点几何定位精度Fig.8 The posterior geometric accuracy of tie points for MEX HRSC images

表2 光束法平差影像残差值Table 2 Image residuals of bundle adjustment results

3.2 数字高程模型精度分析

光束法平差之后，可以利用ASP的stereo、point2dem等工具软件自动生成DEM。LRO NAC原始影像之间几何变形较小，可以直接用于密集匹配；而MEX HRSC影像由于摄影角度差异大，影像之间几何变形较大，直接在原始影像上匹配比较困难，因此可以先纠正为正射影像，然后在正射影像上匹配生成DEM。LRO NAC测试影像选用M1303619844LE与M1303640934LE组成立体像对，而MEX HRSC测试影像选用h7289 nd2、s12、s22（即下视、前视、后视）组成立体像对。自动生成的DEM与交会残差见图9～10，LRO NAC立体影像生成的DEM重采样为5 m格网间距，MEX HRSC立体影像生成的DEM重采样为50 m格网间距。分析LRO NAC与MEX HRSC立体影像的DEM交会残差图可知，大部分区域的交会残差值为1～2个像素，表明影像匹配精度与光束法平差精度较高。受光照不足、阴影、纹理稀疏等因素影响，立体影像中也有少部分区域未匹配出同名点，例如图9（a）左下角的撞击坑（图中显示为白色无效值区域）。另外，图9 （b）中DEM交会残差图显示出有规律的条纹形状，可能是由卫星平台的高频颤振引起，本文未对颤振做进一步处理。

图9 LRO NAC立体影像自动提取DEM结果Fig.9 The DEM results of LRO NAC stereo images

图10 MEX HRSC立体影像自动提取DEM结果Fig.10 The DEM results of MEX HRSC stereo images

由于月球、火星表面通常缺乏高精度的绝对控制数据，行星摄影测量的绝对定位精度通常依赖于探测器定轨、定姿数据的精度。利用行星遥感影像生成DEM时的交会残差及连接点几何定位精度能够表明立体影像之间的内部符合精度较好，实际上是对光束法平差后相对定位精度的评价。为了提升绝对定位精度，可以使用激光测高数据获取的DEM添加高程控制，或者将激光测高足印数据与影像进行融合处理，即将高精度的激光测距信息作为观测值与影像观测值进行联合平差，进一步提升绝对定位精度。

3.3 正射影像图精度分析

利用光束法平差后的精化外方位元素以及自动生成的高分辨率DEM纠正生成DOM，结果如图11～12所示。图11（a）是M1303619844LE/RE构成的正射影像镶嵌图，影像分辨率为0.9 m，图中蓝色方框标注了图11（b）、（c）局部区域影像拼接效果（M1303619844LE与M1303640934LE重叠区域），红色五角星标注了“嫦娥四号”着陆点以及“玉兔二号”月球车在影像获取时刻的位置，具体可参见图11（d）。图12（a）是利用h7289 nd2影像制作的正射影像图，影像分辨率为14 m，影像中间是Jezero撞击坑，图中蓝色方框标注了右侧图12（b）、（c）局部区域影像拼接效果（h7289 nd2与s12重叠区域），左侧为nd2通道影像，右侧为s12通道影像。由正射影像拼图结果可知，使用光束法平差精化的外方位元素以及生成的高分辨率DEM可以显著降低正射影像图的几何拼接误差。

图11 “嫦娥四号”着陆区LRO NAC正射影像图Fig.11 The digital orthophoto map of LRO NAC images for the landing site of Chang’E-4

图12 Mars 2020预选着陆区MEX HRSC正射影像图Fig.12 The digital orthophoto map of MEX HRSC images for the candidate landing site of Mars 2020

4 结 论

本文针对线阵推扫式行星遥感影像摄影测量处理过程中的正射影像生成与控制网构建步骤进行算法优化，自主研发了快速几何纠正与连接点提取软件模块，显著提升了线阵推扫式行星遥感影像的摄影测量处理效率。但是，自主研发的软件模块目前仅能支持LRO NAC与MEX HRSC的相机模型，下一步需要结合月球、火星全球遥感制图的需求，支持更多的深空探测任务数据，并采用GPU计算、集群处理等技术进一步优化算法性能，以适用于大数据量处理的需求。

中国的深空探测将持续向着“开放、合作、共建、共享”的方向发展，而科学数据的高效处理是深空探测工程任务与行星科学研究的基础。当前，综合利用国内外多探测任务数据制作更高精度、更高分辨率的行星制图产品仍然面临诸多技术难题。针对行星遥感影像特点，自主研发一系列的摄影测量软件工具是深空探测遥感测绘中必不可少的基础性工作。