Landsat 8影像像元地理坐标计算

张浩;葛小青;冯旭祥

【摘 要】在遥感卫星数据处理中,计算影像中像元对应地面点的地理坐标不仅是遥感影像数据编目的主要工作之一,同时也是各级几何校正的基础,其计算精度直接影响数据产品处理的几何精度.本文在充分了解Landsat 7卫星和SPOT5卫星像元地理坐标计算算法的基础上,对USGS提出的Landsat 8卫星像元地理坐标计算算法进行了系统性整理和研究,并比较了3种传感器间扫描方式和探测器排布的差异造成计算方法的差异,最后以一景Landsat 8的模拟数据为样本,以GLS2005地面控制点库为参考,将本文介绍算法得到的坐标与控制点库坐标作了对比分析,结果验证了该算法的有效性.

【期刊名称】《遥感信息》 【年(卷),期】2013(028)005 【总页数】7页(P52-58)

【关键词】Landsat 8;视线矢量;像元地理坐标 【作 者】张浩;葛小青;冯旭祥

【作者单位】中国科学院遥感与数字地球研究所,北京100094;中国科学院大学,北京100094;中国科学院遥感与数字地球研究所,北京100094;中国科学院遥感与数字地球研究所,北京100094 【正文语种】中 文 【中图分类】TP79

1 引 言

美国的陆地资源系列卫星(Landsat)对地球进行了长达30余年的连续观测，获取了大量的地球遥感影像数据。这些数据被广泛应用于各类资源调查、生态环境监测、城乡规划与建设、重大自然灾害监测以及农林畜牧业等众多领域。在Landsat系列卫星中，从Landsat 1到Landsat 4已经不再运行，Landsat 6发射失败，目前只有Landsat 5和Landsat 7继续在轨运行。而Landsat 5已发射近30年，其设备严重老化，各项性能不断下降，由于电子元件的急剧退化，已于2011年11月停止成像。Landsat 7在2003年扫描行矫正器(Scan Line Corrector,SLC)发生严重故障后，目前只能得到缺损的图像数据。鉴于此，美国地质调查局(USGS)于2013年2月11日成功发射了Landsat 8卫星，以延续Landsat系列卫星的对地观测任务。Landsat 8卫星的整体外观如图1所示，卫星轨道高度为705km，轨道周期为98.2分钟，重访周期为16天。与之前Landsat系列卫星不同的是，Landsat 8首次采用推扫式成像扫描方式，卫星搭载OLI(Operational Land Imager)和TIRS(Thermal Infrared Sensor)两个传感器，其成像波段参数如表1所示[1]。

图1 Landsat 8整体外观图表1 Landsat 8卫星传感器波段参数

波段号波段名称波长范围(单位:μm)地面分辨率(单位:m)波段范围传感器Band 1海岸/气溶胶(Coastal/Aerosal)0.433～0.45330可见光/近红外OLIBand 2蓝(Blue)0.450～0.51530可见光/近红外OLIBand 3绿(Green)0.525～0.60030可见光/近红外OLIBand 4红(Red)0.630～0.68030可见光/近红外OLIBand 5近红外(NIR)0.845～0.88530可见光/近红外OLIBand 6短波红外1(SWIR-1)1.560～1.66030短波红外OLIBand 7短波红外2(SWIR-2)2.100～2.30030短波红外OLIBand 8全色(PAN)0.500～0.68015可见光/近红外OLIBand 9卷云

(Cirrus)1.360～1.39030短波红外OLIBand 10长波红外1(LWIR-1)10.30～11.30100热红外TIRSBand 11长波红外2(LWIR-2)11.50～12.50100热红外TIRS

在遥感卫星数据处理中，计算遥感影像中像元对应地面点的地理坐标不仅是遥感影像数据编目的主要工作之一，同时也是各级几何校正的基础，其计算精度直接影响数据产品处理的几何精度。之前发射的Landsat系列卫星中，其传感器的扫描方式均为双向扫描式，与Landsat 8的两个传感器所采用的扫描方式存在明显差异，而扫描方式又直接决定了传感器坐标系下像元成像视线(Line of Sight,LOS)矢量的计算方式。因此，USGS针对Landsat 8成像特点，提出了新的像元地理坐标计算算法。本文在充分了解Landsat 7、SPOT5卫星像元地理坐标计算算法的基础上，详细比较了Landsat 7、SPOT5和Landsat 8在传感器坐标系下LOS矢量计算方法的差异，并根据Landsat 8传感器的成像模型特点，详细整理了其图像像元地理坐标的计算过程。 2 像元地理坐标计算

像元地理坐标计算是将像元在图像空间中的平面坐标(x,y,h)转化为该像元对应地面点大地经纬度坐标(λ,φ，h)的过程[2]，如图2所示。 图2 像元地理坐标计算示意图

Landsat 8传感器的像元地理坐标计算过程是根据像元所在的扫描行号、波段序号、传感器芯片组(Sensor Chip Assemblies,SCA)编号以及像元对应探测器在SCA内的序号，计算该像元对应地面点的经纬度坐标。其计算的完整流程图如图3所示[3]。

下面以OLI传感器为例介绍USGS提出的Landsat 8像元地理坐标计算算法的详细步骤，该算法同样适用于TIRS传感器图像像元地理坐标的计算。 图3 像元地理坐标计算流程图

2.1 计算像元成像时间

像元地理坐标计算过程中需要用到像元成像时刻的卫星姿态和星历数据，而这些数据可根据像元成像时间对已有卫星姿态和星历数据插值计算获得。因此，首先要计算出像元的成像时刻t，具体计算公式为： (1)

其中，Tf为该像元所在帧(所有探测器成像一次所产生的数据，多光谱波段每帧对应一行图像数据，全色波段每帧对应两行图像数据)的完成时刻，其值记录在辅助数据中；Ts为探测器成像结束后的采样稳定时间；Ti为探测器曝光时间，Nl为像元所在图像的行号，n为每帧图像数据的行数，Nf为像元所在的帧号，Tsample为每行数据的采样时间。 2.2 LOS矢量计算

像元成像LOS矢量的计算与传感器的成像扫描方式密切相关，对不同的扫描方式，其LOS矢量计算方法不同。LOS矢量最初被定义在传感器坐标系中，其计算方式与探测器在焦平面坐标系的相对位置有关，在焦平面坐标系中定义了各波段探测器相对于光学轴的偏移量，并最终用于产生单个探测器成像的LOS矢量。目前，全球范围内陆地卫星的成像方式主要有双向扫描式(Landsat 1～7等)和推扫式(SPOT1～5/Landsat 8等)两种。本小节详细介绍了Landsat 7、Spot5和Landsat 8所搭载传感器像元成像LOS矢量的计算方法，并详细比较了不同计算方法之间的差异。

2.2.1 Landsat 7卫星LOS矢量计算

Landsat 7的ETM+(Enhanced Thematic Mapper Plus)传感器是扫描镜双向摆动扫描方式的典型代表，其探测器在焦平面坐标系中的位置如图4所示，其中1、2、3、4、8波段分布在主焦平面上，5、6、7波段分布在冷焦平面上。

图4 Landsat 7探测器在焦平面分布示意图

每个探测器对应的LOS矢量与定义在传感器坐标系下的观测角度有关[4]，观测角度可分为平行扫描方向的角度(along_angle)和垂直扫描方向的角度

(across_angle)。由于ETM+传感器的扫描镜是双向摆动扫描，因此平行扫描方向角度又分为正向平行扫描方向角度和反向平行扫描方向角度，同理垂直扫描方向角度也分为正向垂直扫描方向角度和反向垂直扫描方向角度。通过探测器在焦平面坐标系中的偏移量便可以计算出其对应的观测角度：正反向平行扫描方向角度可由像元在扫描行内的时间、整个扫描行的扫描时间、扫描镜的起止角度、像元所在波段中心在平行扫描方向的偏移和该波段奇数探测器的偏移计算得出；正反垂直扫描方向的角度则由探测器所在扫描行中的位置及该扫描行对应的SLC角度决定。由于ETM+中每个探测器在传感器坐标系中的观测角度随着扫描镜的旋转而变化，因此其LOS矢量是一组扫描时间的函数，LOS矢量值随着扫描时间做周期性变化。 根据平行扫描方向角度和垂直扫描方向角度便可以计算出该探测器在传感器坐标系下的矢量：

x=sin(across_angle)*cos(along_angle)y=sin(along_angle)z=cos(across_angle)*cos(along_angle) (2)

最后归一化记单位矢量为 (3)

2.2.2 SPOT5卫星LOS矢量计算

法国SPOT5卫星的HRG(High Resolution Geometric)传感器是推扫式成像方式的典型代表，其每个探测器对应的LOS矢量由其定义在传感器坐标系中的沿轨方向观测角度(ΨX)和跨轨方向观测角度(ΨY)来决定，见图5。由于HRG传感器具备

侧摆能力，因此在计算探测器观测角度时应考虑区域选择镜(Strip-Selection Mirror,SSM)的旋转位置，所以其LOS矢量值依赖于SSM的旋转角度，但在SSM旋转角度不变的情况下，同一探测器成像LOS矢量在传感器坐标系下的值也是固定的。辅助数据中记录了每一个探测器的观测角度，其矢量表示为： x=-tan[(ΨY)p]y=tan[(ΨX)p]z=-1 (4)

其中，(ΨX)p为第P个探测器沿轨方向的观测角度，(ΨY)p为第P个探测器跨轨方向的观测角度。最后归一化记单位矢量为参见式(3)。 图5 传感器坐标系下的视线矢量 2.2.3 Landsat 8卫星LOS矢量计算

Landsat 8的OLI传感器同样采用推扫式成像方式，但不同于SPOT5的HRG传感器将同一波段的所有的探测器在焦平面上依次排开，OLI的焦平面上分两行排列着14组SCA，奇偶交错依次排开，序号相邻的两个SCA间在跨轨方向上有一定的重叠。每个SCA上含有9个成像波段和3个定标波段，每个波段的探测器又可分为主探测器和辅探测器，主探测器个数与辅探测器个数相同，正常情况下主辅探测器相间交错成像，每个探测器的偏移量均记录在校正参数文件(Calibration Parameter File,CPF)中，卫星地面控制系统可根据实际情况在每对主辅探测器之间选择成像的探测器。对于一个多光谱波段每个SCA中含有494个主探测器，对于全色波段每个SCA中含有奇偶两排探测器，每排中有988个主探测器[5]，如图6所示[6]。

图6 OLI探测器在焦平面上的分布示意图

鉴于OLI各波段探测器的分布特点和卫星发射之后几何定标的方便，OLI传感器采用三阶勒让德多项式(Legendre Polynomials)来计算探测器成像LOS矢量在传感器坐标系中的值，参见式(5)。计算LOS矢量时以SCA为单位，每个SCA对应两

组勒让德多项式系数，分别为沿轨方向系数和跨轨方向系数。这两组系数最初在卫星发射之前通过实验室几何定标获得，并记录在CPF中。 x=a0+a1*n+a2*(1.5*n2-0.5)+ a3*(n*(2.5*n2-1.5))+ifov_x*x_shift y=b0+b1*n+b2*(1.5*n2-0.5)+ b3*(n*(2.5*n2-1.5))+ifov_y*y_shift z=1 (5)

其中，ai(i=0,1,2,3)为沿轨方向的勒让德多项式系数，bi(i=0,1,2,3)为跨轨方向的勒让德多项式系数，n为该探测器序号在整个SCA中的归一化值，取值范围为[-1,1]，ifov_x、x_shift为该探测器沿轨方向的瞬时视场角和偏移量，ifov_y、y_shift为该探测器跨轨方向的瞬时视场角和偏移量。最后归一化记单位矢量为参见式(3)。

卫星发射之后，传感器采集数据由Landsat 8卫星数据处理系统经过精确地形校正生成L1T(Level 1 Terrain)数据产品，并以SCA为单位计算出已配准控制点像元在L1T图像和标准参考图像(高分辨率数字正射影像镶嵌图或比OLI更高分辨率传感器的几何校正数据产品)中的偏差，将偏差加权最小二乘拟合，计算出勒让德系数的修正量，进而修正已有系数[7]。采用这种方法可以根据最新卫星数据不断地对传感器进行重新几何定标，从而有效保证了数据处理的几何精度。 2.3 卫星姿态计算

卫星的姿态运动是卫星绕自身质心的转动运动[8]，在本文介绍的算法中，使用姿态数据构造出从导航坐标系到轨道坐标系的转换矩阵。Landsat 8卫星用roll、pitch、yaw分别表示卫星姿态中的翻滚角、俯仰角和偏航角。姿态数据采集频率为50Hz，可以根据已有姿态数据线性插值得到t时刻的姿态数据。

2.4 卫星位置和速度计算

卫星星历数据中记录了卫星的位置和速度信息，Landsat 8的星历数据定义在地心地固(Earth Center Earth Fixed,ECEF)坐标系下。星历数据的采集频率为1Hz，使用3个t时刻之前的星历数据和一个t时刻之后的星历数据，采用拉格朗日插值得到t时刻的星历数据，卫星位置和速度数据的拉格朗日插值公式见式(6)和式(7)： (6) (7)

其中，为卫星位置矢量，为卫星速度矢量，ti为联系位置和速度的统一时刻。 2.5 LOS矢量坐标系转换

由于Landsat8的LOS矢量最初被定义在传感器坐标系，因此需要被转换到导航参考坐标系，然后使用姿态数据将LOS转换到轨道坐标系，最后再被转换到ECEF坐标系下，具体转换过程如下所示： (8)

其中，Ms2a为传感器坐标系到导航坐标系的转换矩阵，其值由卫星发射前的几何定标数据确定；Ma2o为导航坐标系到轨道坐标系的转换矩阵，其值由卫星姿态数据构造出；Mo2e为轨道坐标系到ECEF坐标系的转换矩阵，其值由卫星星历数据构造出；为LOS矢量在ECEF坐标系下的值。 2.6 卫星位置调整

在坐标计算过程中最终使用的是传感器的位置信息，而星历数据中记录的是卫星质心的位置信息，因此需要根据卫星质心相对于传感器的偏置调整卫星位置。首先将质心相对于传感器的位置矢量经过坐标变换转换成ECEF坐标系下的位置矢量然后

将该位置矢量与原卫星位置矢量相加得到卫星新位置矢量具体计算过程为： (9) (10)

2.7 地面点地心经纬度计算

像元成像点的坐标即为LOS与地球相交地面点的坐标，如图7所示，其运算关系如式(11)所示：

图7 LOS矢量与地面相交示意图 (11)

其中，分别为卫星和地面点在ECEF坐标系下的位置矢量，m为矢量的模，该值可根据矢量与矢量的夹角和地球半径利用三角余弦定理求得。

同时由于卫星和地球之间存在相对速度，因此会带来视扰差，所以需要对LOS矢量进行校正以消除误差，该过程被首次引入到Landsat系列卫星的像元地理坐标的计算中。具体校正方法为： (12) (13)

其中，Ωe为地球自转角速度，分别为卫星和地面点在ECEF坐标系下的速度矢量，vlight为光速，为调整之后的LOS矢量。

最后归一化将单位矢量重新记为并带入式(11)重新计算地面点坐标。将表示为则地面点的地心纬度φ和地心经度φ可分别表示为：

(14) (15) 2.8 光速校正

卫星高度大概为705km，地面成像点的反射光到达卫星传感器需要一定的时间，这个时间延迟会对计算地面点的坐标带来一定误差，因此需要进行校正，以消除其影响，其校正过程为： (16) (17)

其中，Δα为反射光由地面点到达传感器时间内地球转动的角度，为校正后地面点在ECEF坐标系下的位置矢量。 2.9 地面点大地纬度计算

成像地面点的地心纬度是指该点和地心连线与赤道面的夹角，而大地纬度是指通过该点椭球面的法线与赤道面的夹角。由于二者之间不存在简单的显式几何转换关系，所以转换过程采用迭代计算的方式进行。 首先根据该点地心纬度计算该点处地球半径R： (18)

其中，a为地球所在椭球体的长半轴长度，e为地球所在椭球体的偏心率，φp为该点的大地纬度(初始值采用其地心纬度φp0)。 其次计算卫星的纬度φS，并求φS与φp的差值Δφ：

(19) Δφ=φS-φp (20)

接着计算卫星高度H： (21)

其中，r为地面点位置矢量的模。检验本次计算的卫星高度H与上次计算结果的差值，若大于规定的阈值(一般取0.001m)，则： (22)

将新计算的φp代入式(18)重新迭代计算，直到卫星高度H与上次计算结果的差值小于规定的阈值或迭代次数超过20次则结束迭代，此时的φp就是转换后地面点的大地纬度，至此成像地面点的大地经纬度坐标计算完成。 3 实验结果

为了说明本文所介绍算法的实用性，使用USGS发布的GLS2005控制点库作为参考依据，对算法的处理结果进行分析。该库是在Landsat系列卫星影像数据的基础上提取建立的[9]，其中80%的控制点定位精度优于30m[10]。同时，该库也作为Landsat8数据几何精校正阶段的控制点来源。本文的实验数据采用USGS发布的一景Landsat 8模拟数据，该景数据在WRS-2坐标系中的Path/Row为024/038。现从GLS2005控制点库中选取均匀分布在该景数据范围内的10个点作为参考，其分布位置如图8所示。同时从该景数据中分别找出每个参考点对应的扫描行号、波段号、SCA编号和像元序号，使用本文介绍的算法计算出其坐标，并与其对应参考点的坐标进行对比分析，结果如表2所示。

图8 样本点位置分布示意图

从表2可以看出，与GLS2005控制点库中的数据相比，使用本文介绍算法计算成像地面点的经度误差和纬度误差换算成地面距离均在亚像素级，可以满足Landsat 8卫星数据预处理阶段的精度要求。

表2 本文计算方法与GCP库比较结果序号GCP点库坐标本文算法计算坐标相对误差经度纬度经度纬度经度差纬度差地面距离(单位:m)1-93°05′53.27″32°23′08.74″-93°05′52.99″32°23′08.64″0.28″-0.10″8.02-92°44′51.42″30°58′36.22″-92°44′51.44″30°58′35.90″-0.02″-0.32″9.73-93°28′29.16″31°30′47.28″-93°28′29.11″31°30′46.63″0.05″-0.65″20.04-93°28′23.94″32°20′36.61″-93°28′23.86″32°20′36.72″0.08″0.11″4.05-91°59′38.73″31°45′05.12″-91°59′39.16″31°45′04.36″-0.43″-0.76″24.86-92°16′52.53″32°06′17.07″-92°16′52.63″32°06′17.02″-0.10″-0.05″3.07-93°49′58.45″31°21′14.57″-93°49′58.46″31°21′15.19″-0.01″0.62″19.28-93°45′31.04″32°04′23.16″-93°45′31.21″32°04′23.30″-0.17″0.14″6.29-93°12′12.86″31°27′41.81″-93°12′12.98″31°27′41.74″-0.12″-0.07″3.910-92°50′12.76″32°13′32.69″-92°50′12.53″32°13′32.50″0.23″-0.19″8.3 4 结束语

本文在充分了解Landsat 7卫星和SPOT5卫星像元地理坐标计算算法的基础上，对USGS提出的Landsat 8像元地理坐标计算算法进行了系统性的整理和研究，并比较了由于不同传感器间扫描方式的差异所带来的计算方法差异，另外在该算法中USGS还考虑到了传感器相对于卫星质心的位置偏移、速度像差、光程时间延迟等，并对这些因素所带来的影响予以校正和消除，使运算流程更科学，计算精度更高，最后以一景Landsat 8的模拟数据为样本，验证了该算法的有效性。本文所介绍的算法在陆地卫星像元地理坐标计算方面具有典型性，对其他型号的陆地卫

星(特别是推扫式传感器)具有重要的参考意义。 参考文献：

[1] IRONS J R,DWYER J L,BARSI J A.The next Landsat satellite:The Landsat data continuity mission[J].Remote Sensing of Environment,2012(122):11-21. [2] RIAZANOFF S,GLEYZES J P,ZOBLER D,et al.SPOT 123-4-5 Geometry Handbook[R].GAEL-P135-DOC-001,2004.

[3] USGS.LDCM CAL/VAL ALGORITHM DESCRIPTION DOCUMENT[G].2013. [4] USGS EROS Data Center.Landsat 7 Image Assessment System (IAS) Geometric Alogrithm Theoretical Basis Document[G].2011.

[5] REINKE N,ZAHN S,HOLLAREN D,et al.USGS Landsat data continuity mission (LDCM) level 0 reformatted archive (L0RA) data format control book (DFCB)[G].Department of the Interior U.S.Geological Survey,2011. [6] SCHOTT J R,RAQUENO R V,RAQUENO N,et al.A synthetic sensor/image simulation tool to support the Landsat data continuity mission (LDCM)[G].Digital Imaging and Remote Sensing Lab,2012.

[7] BECKMANN T,ENGEBRETSON C,ROCKVAM T,et al.USGS Landsat data continuity mission (LDCM) image assessment subsystem (IAS) Software design document (SDD)[G].Department of the Interior U.S.Geological Survey,2012.

[8] 章仁为.卫星轨道姿态动力学与控制[M].北京:北京航空航天大学出版社,1998. [9] USGS Publications Warehouse.Extraction of GCP chips from GeoCover using Modified Moravec Interest Operator (MMIO) algorithm[G].2006. [10] USGS EROS,GLCF,NASA Headquarters.Assessment of the NASA-USGS

Global Land Survey(GLS) Databases[G].2012.