卫星遥感影像产品快视图地理校正研究

李 想，高 磊

(1.江苏星月测绘科技股份有限公司南京雨花台分公司，江苏 南京 210000；2.江苏省测绘工程院，江苏 南京 210013)

0 引 言

近十年来，国产遥感卫星发射数量不断增加，空间、时间和光谱分辨率不断提升，用户对高分辨率遥感数据的需求日益提升。目前我国已有包含高分、资源、环境减灾等多个系列在内的陆地观测卫星近30颗，形成了多尺度、多观测要素、多级别产品的业务化对地观测体系，卫星遥感影像产品分发服务机制也日趋完善[1-3]。

终端用户订购获取的卫星遥感影像产品除影像数据外，通常还包括元数据、成像几何模型参数文件、快视图等[4]。其中快视图是栅格影像经拉伸增强后的真彩色或假彩色图片，用于直接显示影像数据呈现的地表环境状态，如地理范围、云雪覆盖、水陆分布等，方便用户快速了解影像质量情况进行选择。据统计，早期的国产陆地卫星遥感影像产品中的快视图都没有提供定位信息，近年来高分系列的PMS、资源一号02d、北京二号等少数产品的快视图逐步增加了定位文件，其它如高分一号和高分六号的WFV、资源三号(01、02、03星)、高分七号、高景一号等遥感影像产品中快视图仍然只是一张普通的JPG格式图片，缺少定位信息，无法结合行政区划或研究区矢量边界范围进行针对目标区域的快速浏览和数据筛选，在实际生产和应用项目中较为不便。遥感产品的元数据文件记录了影像产品的成像时间、太阳角度、观测角度、四至坐标、像元大小或影像尺寸等基本信息[5]。本文在分析图像坐标定位文件结构的基础上，探讨了利用元数据提供的影像数据坐标信息计算定位参数的方法，并通过Python编程语言结合空间数据转换库(Geospatial Data Abstraction Library，GDAL)开发了快视图定位文件生成工具，以期为相关研究和生产作业人员提供参考。

1 研究方法

本文的研究数据源是几种主要的国产中高分辨率光学遥感影像产品，其中高分、资源和环境等系列影像来源于中国资源卫星应用中心，高景影像由北京航天世景信息技术有限公司提供。

1.1 定位文件结构分析

图像以栅格数据形式存储，栅格图像中的每个像元都具有一个对应于特定行号和列号的像素坐标。要在正确的地理位置显示图像，需要建立特定的变换关系，以将像元坐标转换为真实世界坐标，通常使用6参数构建的仿射变换关系进行转换[6]：

(1)

写为方程组的形式：

x′=Ax+By+C

(2)

y′=Dx+Ey+F

(3)

式中，x、y分别为像素坐标中的列号和行号，x′、y′为该像元对应的地理坐标；A、E分别为像素在地图坐标系X和Y方向上的尺寸，采用地图单位；B和D为旋转项；C和F分别左上角像素中心点的地图坐标。由于图像坐标系和地图坐标系的起点分别位于左上角和左下角，Y正方向相反，因此E为负值。

仿射变换信息可以存储在图像文件的文件头内部，也可以存储在单独的普通文本文件中，这个文本文件就是图像坐标定位文件，也称世界文件(world file)[7]，仿射变换的6个参数以十进制数字的形式，按A、D、B、E、C、F的顺序，依次分行存储在文件中。坐标定位文件的主文件名与图像文件相同，后缀名分两种：在图像文件的扩展名最后添加字母w，或使用图像文件扩展名中的第一个和第三个字符，最后再添加字母w，如jpg格式图像的坐标定位文件后缀名可以是jpgw或者jgw，通常后者更为常用。图1为2021年1月11日某景高分一号PMS影像的坐标定位文件内容，产品号为L1A0005387434。

图1 图像坐标定位文件示例

1.2 定位参数解算原理

图像定位校正与仿射变换参数的关系如图2所示，X-O-Y为地图坐标系，影像产品的元数据中记录了影像四角的坐标信息，以此可以计算影像在X方向的距离w和Y方向的距离h，均为地图单位：

图2 图像定位校正与仿射变换参数关系示意图

ΔX1=xur-xul

(4)

ΔY1=yur-yul

(5)

(6)

同理，

ΔX2=xur-xlr

(7)

ΔY2=yur-ylr

(8)

(9)

并计算得到：

(10)

Rx=w/col

(11)

Ry=h/row

(12)

式中，xur、yur、xul、yur分别为图像右上角和左上角的地理坐标，row、col为图像行数和列数，Rx、Ry为像元在X和Y方向上的尺寸，t为图像旋转角度，以X轴逆时针方向为正。由此可以得到：

A=Rx×cost

(13)

B=Ry×sint

(14)

D=Rx×sint

(15)

E=-Ry×cost

(16)

(17)

(18)

2 编程实现

遥感影像产品按处理的程度进行分级，其中L0级和L1级产品未经几何校正，地理坐标信息存储在成像几何模型参数文件中，影像文件本身不带地理信息，大部分国产卫星遥感影像如高分、资源系列等都按L1级产品分发；L2级产品是指经过系统几何校正的影像产品，影像文件带有地理坐标，如环境、中巴系列等[8]。为适应不同级别遥感产品快视图定位的需要，本文利用Python编程语言和GDAL转换库，从元数据解析出发，结合影像四角地理坐标和快视图图像坐标构建控制点，快速计算仿射变换参数并生成坐标定位文件，主要流程如图3所示。

图3 坐标定位文件参数计算流程

2.1 元数据解析

元数据文件通常以可扩展标记语言(EXtensible Markup Language，XML)格式随卫星影像产品下发，XML是一种用于标记电子文件使其具有结构性的标记语言，以纯文本格式存储，便于读取、记录、调试，使不同系统、不同程序之间的数据共享变得更加简单。Python自带的Minidom库可以方便地对XML文件进行解析，提取出所需要的数据[9]，主要过程如图4所示。至此即可获取左上角坐标(ulx、uly)、右上角坐标(urx、ury)、左下角坐标(llx、lly)、右下角坐标(lrx、lry)共计四对值。

图4 元数据解析主要代码

2.2 地理变换参数解算

地理变换参数涉及到坐标的缩放、旋转和平移，由于地理校正后的影像存在扭曲、变形等原因，并非常规矩形，利用角点坐标不同组合计算按理论方法计算得到的值会存在差异，存在多组已知点时，需要用最小二乘法解算最佳变换参数，计算过程较为繁琐。GDAL是一个在X/MIT许可协议下的开源栅格空间数据转换库，使用抽象数据模型来解析它所支持的数据格式，抽象数据模型包括数据集、坐标系统、仿射地理坐标转换、大地控制点等，具有强大的地理数据处理能力[10]，利用其提供的函数可以通过已知控制点计算地理变换参数，具体方法如图5所示。

图5 地理变换参数解析主要代码

3 结果评价

计算地理仿射变换参数并生成坐标定位文件后，快视图就可以在地理信息系统或遥感平台软件里进行地理定位，在真实坐标系统下显示。本研究以产品号为GF1_WFV1_E121.6_N31.2_20211126_L1A0006081870的高分一号WFV影像产品为例，对快视图地理校正的效果进行评价。从目视对比看，恢复地理信息后的快视图(左半幅)与实际影像数据(右半幅)的空间位置基本一致(图6)。由于产品快视图的分辨率在160 m左右，为影像分辨率的十分之一，无法准确定位同名地物的位置，本文以影像四角像素的外侧角点为检查点，计算校正偏差(表1)，计算得到X、Y方向上偏差的中误差分别为0.010 5、0.000 5(单位/°)，偏差主要体现在经度方向上，约6个像素，该精度已能够满足粗略定位和范围显示的需要。

图6 地理校正后的快视图与影像文件叠加对比

表1 检查点校正偏差(单位/°)

4 结 语

“十四五”期间，我国将有计划地陆续发射多颗陆地观测卫星，进一步推动遥感业务化应用向市场化、产业化发展。遥感影像产品直接面向终端用户，包括影像文件、元数据、快视图等内容，其中快视图是影像文件的缩略图，为普通的图片格式，便于用户直接浏览，了解影像文件的质量概况。本研究利用影像产品附带元数据的四角坐标信息构建控制点，通过GDAL空间数据处理模块计算坐标定位参数，生成坐标定位文件，实现对快视图的地理校正，恢复快视图的地理信息，为用户进行影像产品存档、筛选、预览、建库等工作提供有效的帮助。