天宫二号卫星宽波段成像光谱仪地理定位方法

杨涛，何玉青，胡秀清

(1.北京理工大学 光电学院，北京 100081；2.国家卫星气象中心，北京 100081)

0 引言

对地观测卫星遥感图像的地理定位是建立星载遥感仪器观测图像的每个像元与地基坐标系中位置关系的过程。传统的遥感图像地理定位方法分为参数法和非参数法[1]。非参数法主要是通过地面点与遥感数据来建立对应关系，因此选取合适的地面控制点是该方法的关键。但是对于运行中的遥感卫星来说，非参数法需要大量的地面点来进行辅助，并且在云覆盖较多的情况下不易分辨地面点，因此不适合用于业务运行中的卫星。该方法主要用于早期星上测量手段缺乏或者星上测量数据丢失的情况，如QuickBrid、印度资源卫星、天绘一号等[2-4]。此外早期的星载高分辨率成像仪器常采用非参数定位方法，发展出了直接线性变换(direct lineaer transform，DLT)[5]、有理多项式系数(rational polynomial coefficients，RPC)[6]等一系列的模型。参数法是根据遥感仪器观测几何及其空间位置和指向，建立观测像元与地面观测位置之间的模型，其主要是利用精确测定的卫星位置、速度、姿态以及遥感仪器的扫描几何和时序等参数来计算出每个像元的经度、纬度和观测角度。目前绝大部分卫星都是使用参数法进行地理定位，如NOAA/AVHRR3[7]、Terra/MODIS[8]、Landsat-7/TM、NPP/VIIRS[9]，其不同之处主要体现在各个遥感仪器内部光学几何结构以及各个卫星轨道计算方法。

天宫二号(tiangong No.2，TG-2)是我国首个真正意义上的空间实验室，于2016年9月15日发射升空，其中包含了空间冷原子钟、三维成像微波高度计和多角度宽波段成像光谱仪等14个精密载荷。宽波段成像光谱仪是其中的有效载荷之一，其主要用于观测获取海洋、大气和陆地的各类参数，用于观测海洋水色海温以及优化天气预报[10]。其对天气预报进行优化主要是利用观测数据可以改善气溶胶反演数据，而这些科学和行业应用需要高质量的地理数据进行辅助，如经纬度和观测角。因此对卫星定位数据的质量和精度提出了较高要求。

宽波段成像光谱仪采用的是多电荷耦合元件(charge-coupled device，CCD)拼接直视推扫的扫描方法。本文采用参数法的原理，针对其特殊的扫描特点，分析了宽波段成像光谱仪的工作原理与几何模型，对整个仪器的地理定位方法进行分析与实验，并对定位结果进行了初步评估。

1 宽波段成像光谱仪的工作原理与几何模型

宽波段成像光谱仪根据工作波段分为三大模块：可见光近红外模块、短波红外模块和热红外模块。宽波段成像光谱仪总视场为42°，在考虑到可选用的探测器不能满足条件的情况下，各模块采用了单相机小视场合成总视场的方案：可见光近红外采用3个子相机；短波和热红外模块各采用2个子相机。同一波段的各个相机在同一平面内扇形展开，以可见光近红外波段的相机为例，单个相机的视场略大于14°，3个相机两两成14°的夹角，合成一个42°的总视场。可见光近红外波段相机的焦平面是由一个1 024×512(空间维×光谱维)面阵CCD构成，3个相机拼接成一幅完整图像。因为相机的视场比夹角略大，因此3个相机成像时会出现重叠区域。TG-2轨道高度约为400 km，宽波段成像光谱仪可见光近红外模块的地面分辨率为100 m，瞬时视场为250 μrad。相机之间存在着32个暗电平。在去处暗电平之后，2个相机之间的重叠区域观测角约为0.4°。根据瞬时视场计算可得重叠区域约为28个像元。

宽波段成像光谱仪安装在卫星的对地面，开口朝向地球，3个模块分层安置，同时进行推扫扫描，扫描方向就是卫星飞行方向，借助于卫星绕地球运行，获取地球的二维景象。宽波段成像光谱仪可见光近红外相机观测地球景象的原理图如图1所示。

图1 宽波段成像光谱仪可见光近红外相机观测地球景象原理图

宽波段成像光谱仪单次采样可以获取4个数据：可见光近红外相机获取的15个通道的图像数据、短波红外相机获取的2个通道的图像数据、热红外相机获取的2个通道的图像数据以及相同时间段内平台探测设备获取的辅助参数数据。图像信息数据包含的是各通道的图像数据以及该数据起始时间、截止时间以及一个总计数时间。辅助参数数据包含定位所需的参数，即卫星轨道GPS位置与速度、卫星姿态角度信息、与GPS姿态数据相对应的协调世界时(UTC)以及其他平台参数，如温度、船下点、轨道高度等。

3个模块的成像原理大致相同，不同之处只是相机数量、地面空间分辨率和各个相机之间安装角度。短波和热红外模块因为采用2个子相机，每个相机的光学视场为22°，两两夹角为21°。短波红外相机的CCD为800×2，地面分辨率为200 m；热红外相机的CCD尺寸为400×2，地面分辨率为400 m。以下所述的定位方法均是以3个可见光近红外波段相机为例，一个可见光近红外数据即为可见光近红外相机获取的图像数据以及与其匹配的平台辅助参数数据。本次开机共获取了600多个可见光近红外数据，其中一半数据为海洋观测区域，另一半数据为陆地区域，以中国为主。

2 宽波段成像光谱仪地理定位方法及流程

本节对宽波段成像光谱仪地理定位方法的整个流程都进行了详细的设计，主要包括读取原始数据并对数据做预处理、通过处理后的数据建立焦平面坐标系视向量与复合转换矩阵、利用视向量计算地理定位信息、通过调整安装矩阵的方法对定位结果进行了修正。定位方法具体流程图如图2所示。

图2 地理定位方法流程图

2.1 原始数据读取与预处理

传统的参数法定位大致需要2组参数，一组是与卫星相关的，主要是卫星轨道位置、速度以及运行时的姿态等；另一组与仪器相关，主要是仪器内部安装几何、探测器焦距像元尺寸等。其中仪器相关的数据在仪器研制过程中就已经进行了设计和测量，可以直接获取。而卫星轨道相关数据，传统上是通过卫星两行报中的描述卫星运行的6个参数，利用卫星轨道模型来计算得到卫星实时位置信息。TG-2自身携带了GPS接收器以及星敏感器，可以实时测量得到卫星的位置以及姿态。因此对于该仪器的地理定位方法中不需要计算轨道参数，直接通过读取L0级数据中的位置和姿态参数进行地理定位。

宽波段成像光谱仪的辅助信息和图像信息分别存在2种不同格式的L0级数据中。因为2种数据的L0级文件格式不同，时间的采样频率与表示形式也均不相同，所以2种数据的时间不是一一对应的关系，而进行卫星定位时需要每帧图像精确的观测时间。所以需要进行预处理，即对输出的数据时间进行统一，保证每一帧图像数据都有辅助数据对应。

仪器主要的时间问题是由于图像信息采样频率、辅助信息采样频率、平台辅助参数的更新频率以及时间表示不同造成的。其中平台辅助参数的更新频率最小，其时间是UTC时间。辅助信息采样频率是平台辅助参数的更新频率的一倍，可见光图像信息采样频率又是辅助信息采样频率的40倍，二者的时间均是以开机时刻作为起始时间零的计数时间。

由于辅助数据包与图像数据包采用的记录时间不同，所以先要进行时间的统一，将辅助数据包采用的计数时间和图像数据包采用的UTC时间都转换为国际原子时。在统一时间单位后，通过时间来对辅助参数数据进行插值，实现辅助数据与图像数据一一对应，插值的流程如下：

①利用图像信息起始时间与截止时间获得当前数据包的总观测时间；

②根据采样频率比以及总时间得到每帧图像的具体精确时间；

③利用每帧图像的具体时间和平台辅助信息的时间进行插值，获得每帧图像对应的具体辅助信息。

2.2 焦平面坐标系视向量及复合转换矩阵构建

地理定位方法的基本原理就是通过仪器焦平面引出对地观测的向量，通过该向量与地球模型的交点来得到XYZ坐标，再经过坐标系转换最终得到经纬度等地理信息。首先要建立焦平面模型，来获取成像光谱仪焦平面坐标系下的视向量。

虽然宽波段成像光谱仪的可见光波段由3个相机组成，但是在地理定位的过程中并不相互影响，因此本文把3个相机作为3个独立的个体，分别建立焦平面模型。单个宽波段成像光谱仪的可见光相机是由一个1 024×512(空间维×光谱维)面阵CCD扫描成像，本文把扫描方向也就是飞行方向作为x轴正方向，穿轨方向为y轴方向，y轴正方向与图像的数据存储方向一致，如图3所示。

图3 成像光谱仪焦平面坐标系

由此可以得到焦平面上每个像元的x、y坐标以及焦距f，可以建立焦平面模型，也就是焦平面视向量ufoc。

(1)

焦平面视向量是根据焦平面坐标系建立，需要将其转换到测地坐标系下，再结合地球模型才可以得到视向量与地球交点，利用交点的坐标信息通过坐标系转换可以获取所需的经纬度信息。其中涉及到望远系统、仪器、卫星本体、轨道、地心惯性(earth centered inertial，ECI)、地心旋转(earth centered rotate，ECR)和大地测量8个坐标系以及一系列的坐标转换矩阵，如图4所示，把这一系列坐标转换矩阵组合可以得到一个复合转换矩阵。

图4 多个坐标系转换图

由焦平面坐标系转换到望远系统坐标系的转换矩阵Ttel/foc主要与仪器的光学成像系统有关，一般由仪器研制方在发射前提前测量得到，通常使用单位矩阵代替。但是对于宽波段成像光谱仪来说，由于其是多CCD拼接成像，对应的焦平面有3个，这3个焦平面模型与望远系统坐标系不尽相同。其中垂直扫描的相机可以使用单位矩阵，而两侧的相机与中间的相机成14°角，根据焦平面坐标系，其扫描方向为x方向，因此两侧的相机等于绕x轴旋转正负14°，绕x轴旋转的公式为：

(2)

式中：θ为旋转的角度。由望远系统坐标系转换到仪器坐标系的转换矩阵Tinst/tel主要与仪器扫描机理有关，该仪器成像方式为推扫成像，没有扫描镜调整方向，因此可以使用单位矩阵来表示这个旋转矩阵。

(3)

仪器坐标系与卫星本体的坐标系之间的关系是由仪器的安装矩阵确定，其是由仪器研制方在场地测试中计算得到的，用于调整仪器安装时带来的误差。可以先将转换矩阵Tsat/inst表示为单位矩阵，在卫星在轨运行之后需要继续进行调整，以修正定位误差[11]。

卫星的姿态决定了卫星坐标系与轨道坐标系之间的关系，该姿态由卫星上的星敏感器提供，表示形式为欧拉角，即滚动(roll)角、俯仰(pitch)角与偏航(yaw)角3个姿态角。其在直角坐标系下对应的就是绕X轴、Y轴与Z轴旋转，而在卫星地理中一般按照ZXY的顺序计算转换矩阵[12]，因此该转换矩阵Torb/sat可以表示为：

(4)

式中：y为偏航角；p为俯仰角；r为滚动角。

轨道坐标系与地心惯性系之间的关系是由卫星的瞬时位置与速度计算得到的，可以通过在ECI坐标系中形成轨道坐标系轴来构建从轨道转换为ECI的旋转矩阵Teci/orb，轨道坐标系的原点为卫星质心，Z轴方向对应的是卫星星下点方向，Y轴方向是卫星角动量方向的负方向，X轴方向由右手准则确定，即Z轴方向与Y轴方向的叉积。卫星星下点方向就是卫星在ECI坐标系下的坐标归一化后的方向，而卫星角动量方向的负方向是垂直于卫星星下点方向与卫星速度方向的所构成的平面，公式如下：

(5)

式中：p表示卫星在ECI坐标系下的瞬时位置；v表示卫星在ECI坐标系下的瞬时速度。

从ECI坐标系到ECR坐标系的变换主要是由于时间不同、地球旋转带来的。此外还要考虑岁差、章动与极移这些变化的影响。这些矩阵可以通过国际天文联合会提供的岁差章动极移模型得到。

将这些转换矩阵构成一个复合转换矩阵T：

T=Tecr/eciTeci/orbTorb/satTsat/instTinst/telTtel/foc

(6)

通过复合转换矩阵与焦平面视向量，就可以得到在ECR坐标系下的视向量模型uecr

uecr=T×ufoc

(7)

2.3 获取地理坐标

X=p+d·u

(8)

式中：p为卫星在ECR坐标系下的位置矢量；d为观测点到卫星的距离，其是由通过地球椭圆模型重新缩放的视向量与卫星位置矢量构成的：

(9)

式中：u′与p′就是根据地球椭圆模型长短轴重新缩放的视向量与卫星位置矢量；a为长轴，b为短轴：

(10)

卫星运行观察地球的几何关系如图5所示。

图5 观测交点几何模型

最后将此位置矢量X(x1,x2,x3)换算为大地测量坐标(lat，lon，h)，此外地形和大气折反射也会对遥感图像地理定位精度造成影响，利用全球数字高程模型(digital elevation model，DEM)进行修正，本文地理定位采用来自美国国家影像与制图局和世界数字化图的30″分辨率DEM数据[13]。利用DEM找到对应区域高度信息，通过迭代的方法利用高度信息对ECR坐标系下视向量进行不断修正直到高度信息相符，把修正后的视向量转换到大地坐标系下得到修正后的经纬度和高度信息。

(11)

2.4 定位修正

因为3个相机的定位数据是单独计算的，所以定位数据也会出现重叠区域。把3个相机获取的定位数据拼接在一起时，需要考虑到重叠区域，拼接方法和图像的拼接方法一致。2个相机之间有28个像元的重叠。考虑到边缘像元分辨率下降的问题，于是在剔除重叠区域的时候，剔除各个相机更边缘的区域，保留了靠近中心的区域。即两边的相机各自剔除14个像元，而中心相机两边各剔除14个像元。

利用拼接后的图像与定位数据，主要选取具有特征的地面点，如海陆边界线和湖泊，进行投影输出。投影结果与高精度海陆分界线模板进行对比，可以发现存在定位结果误差，效果如图6所示。通过调整安装矩阵可以减少误差提高定位精度。安装矩阵对图像的修正效果体现在图像上是卫星图像绕坐标轴旋转后的结果，坐标轴方向与焦平面坐标系的正方向一致，绕X轴旋转会使图像左右偏移，绕Y轴旋转会使图像前后偏移，绕Z轴旋转会使图像产生一定的旋转偏移。相应的我们可以通过计算图像中偏移像元数量来对这几个偏移角度进行迭代调整，直到图像定位误差达到最小。

图6 TG-2宽波段成像光谱仪未修正定位结果局部投影图

可见光近红外波段探测器的瞬时视场角ωIPOV为0.014 3°。本文挑选了50个带有明显地标特征的图像数据进行统计，主要统计数据中心区域海岸线模板与图像中海岸线的像素偏差。统计结果为在星下点附近左右方向平均偏移约13.41个像元，上下方向平均偏移约6.53个像元，整幅图像两侧边缘平均有个5.26像元的转动，把像元偏移量带入角度旋转公式：

(12)

式中：Δny为左右平移量；Δnx为上下平移量；Δnz为旋转平移量；L为每行的采样点3 016个。最终计算得到绕X轴旋转角φ为0.192°、绕Y轴旋转角μ为0.093°、绕Z轴旋转角ψ为0.099°。

把相应的角度带入公式(4)可以计算出修正后的安装矩阵，修正前后的定位改进效果可以参看图7。

图7 TG-2宽波段成像光谱仪局部放大图像定位结果图

在定位精度提高方面，除了对安装矩阵进行调整的修正方法外，还可以采用基于视向量的修正方法[14]，即通过统计像元在扫面和穿轨方向上的误差来建立像元与误差的模型，利用模型修正视向量，从而减少定位误差提高定位精度。本文也利用这一方法对第一种方法采用的50个数据进行了修正，其中，修正模型为视向量误差随影像列号成线性变化的一次线性模型[15]。通过分别统计2种方法的修正结果，我们发现修正安装矩阵后平均误差像元数为2.07，采用视向量修正模型后平均误差像元数为1.99，进而每幅定位图像与海陆模板的误差仍处于2个像元左右，由此这2种方法的修正效果差异不明显。与此同时，由于修正安装矩阵的方法更容易实现，因而本文在后续的分析中仍然采用了对安装矩阵进行调整的修正方法。

3 定位实验与结果

3.1 实验数据

宽波段光谱成像仪在TG-2入轨不久就开机获取数据，此次开机运行持续了6个多月的时间。主要包含有2016年9月、10月、12月以及2017年1月、2月、3月的数据。12月份的数据最多每天都有数据，其他月份仅有十几天有数据，其中每天的数据也不是整天所有轨道的数据都存在，以几个或几十个可见光近红外数据为主。

在这600多个可见光近红外数据的基础上，本文做了进一步的筛选。因为数据中存在了整个图像都是海洋的这种数据，这类数据对定位结果的判断没有帮助，因此剔除了这类数据。剩下的共301个可见光近红外数据以中国陆地区域为主或者是包含海岸线区域，把这些数据作为有效数据。因为一个数据包含的图像数据的垂直像元数量很多，大约有为50 000多个。为了显示与存储的方便，本文以4 800列为单位，将一个数据分解成了若干份。即一份实验用的图像数据为3 016×4 800，如图8所示。

图8 TG-2可见光近红外相机原始图像

3.2 实验结果与误差分析

本文从6个月的300多个有效可见光近红外数据中每个月都随机挑选了几个可见光近红外数据进行地理定位，利用修正后的定位结果与高精度海陆分界线模板比较，结果如图7所示，图中白色线条为高精度海陆分界线模板的投影。从局部放大的图像中可以看到海陆分界线与图像基本契合，二者之间存在几个像元的误差。

为了检验定位算法和修正矩阵的适用性，本文对所有的300多个有效可见光近红外数据都进行地理定位，并且进行投影输出。与平台辅助信息提供的船下点经纬度比较，误差在0.01°以内。本文把图7中对应数据星下点的经纬度与平台辐射信息提供的经纬度进行了对比，如表1所示。

此外本文还把海陆模板叠加在投影图像上，统计了二者相差的像元数量，考虑到海陆模板的精度问题，每幅定位图像与海陆模板的误差在2个像元左右，如表2所示。此外本文对宽波段光谱成像仪的其他2个波段的图像也进行了地理定位，由于短波红外波段相机与热红外波段的相机空间分辨率都低于可见光波段的相机，定位精度有所提升，热红外波段相机误差在1个像元左右。

表1 实验定位信息与平台辅助信息比较表

表2 定位结果误差统计表

定位结果中主要误差是由卫星观测时间的偏差、卫星姿态测量精度、卫星位置测量精度、发射后仪器安装偏差等造成的。其中安装偏差是误差的主要来源，其是由于安装误差和卫星发射过程中的剧烈震动，以及发射后的环境的影响，导致仪器与卫星之间的位置关系发生改变，即使采用地面测试得到安装矩阵依旧无法完全消除定位误差，需要通过对图像进行统计分析，调整安装矩阵来进行修正。另外，卫星姿态测量和卫星位置测量精度还存在随机误差，它会随着时间的变化而变化，造成图像出现一定的平移与旋转，这种变化不易进行模拟。

4 结束语

天宫二号宽波段成像光谱仪为多CCD拼接直视推扫模式。本文详细设计了基于该仪器扫描模式的参数型定位方法，其中包括数据读取与预处理、建立焦平面视向量与坐标系转换矩阵、地理参数计算以及定位修正这几个过程。本文通过卫星的空间位置姿态参数，建立了仪器焦平面与地理空间位置的对应模型。通过大量实际观测的定位结果建立旋转角度与变形数量的关系，修正安装矩阵，提高模型定位精度。利用该模型实现了对TG-2卫星宽波段成像光谱仪图像的地理定位。利用图像的地理位置进行投影后，与高精度海陆模板进行比较，定位精度约为2个像元。该研究成果为后期的TG-2宽波段成像光谱仪的图像应用，以及海洋、大气和陆地地球物理参数的定量反演奠定了坚实的基础。