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

1、Landsat 8影像像元地理坐标计算Landsat 8影像像元地理坐标计算张浩1，2，葛小青1，冯旭祥1【摘 要】摘要：在遥感卫星数据处理中，计算影像中像元对应地面点的地理坐标不仅是遥感影像数据编目的主要工作之一，同时也是各级几何校正的基础，其计算精度直接影响数据产品处理的几何精度。本文在充分了解Landsat 7卫星和SPOT5卫星像元地理坐标计算算法的基础上，对USGS提出的Landsat 8卫星像元地理坐标计算算法进行了系统性整理和研究，并比较了3种传感器间扫描方式和探测器排布的差异造成计算方法的差异，最后以一景Landsat 8的模拟数据为样本，以GLS2005地面控制点库为参考，将

2、本文介绍算法得到的坐标与控制点库坐标作了对比分析，结果验证了该算法的有效性。【期刊名称】遥感信息【年(卷),期】2013(028)005【总页数】7【关键词】Landsat 8;视线矢量;像元地理坐标【文献来源】1 引 言美国的陆地资源系列卫星(Landsat)对地球进行了长达30余年的连续观测，获取了大量的地球遥感影像数据。这些数据被广泛应用于各类资源调查、生态环境监测、城乡规划与建设、重大自然灾害监测以及农林畜牧业等众多领域。在Landsat系列卫星中，从Landsat 1到Landsat 4已经不再运行，Landsat 6发射失败，目前只有Landsat 5和Landsat 7继续在轨运

3、行。而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首次采用推扫式成像扫描方式，卫星搭

4、载OLI(Operational Land Imager)和TIRS(Thermal Infrared Sensor)两个传感器，其成像波段参数如表1所示1。在遥感卫星数据处理中，计算遥感影像中像元对应地面点的地理坐标不仅是遥感影像数据编目的主要工作之一，同时也是各级几何校正的基础，其计算精度直接影响数据产品处理的几何精度。之前发射的Landsat系列卫星中，其传感器的扫描方式均为双向扫描式，与Landsat 8的两个传感器所采用的扫描方式存在明显差异，而扫描方式又直接决定了传感器坐标系下像元成像视线(Line of Sight,LOS)矢量的计算方式。因此，USGS针对Landsat 8成像

5、特点，提出了新的像元地理坐标计算算法。本文在充分了解Landsat 7、SPOT5卫星像元地理坐标计算算法的基础上，详细比较了Landsat 7、SPOT5和Landsat 8在传感器坐标系下LOS矢量计算方法的差异，并根据Landsat 8传感器的成像模型特点，详细整理了其图像像元地理坐标的计算过程。2 像元地理坐标计算像元地理坐标计算是将像元在图像空间中的平面坐标(x,y,h)转化为该像元对应地面点大地经纬度坐标(,，h)的过程2，如图2所示。Landsat 8传感器的像元地理坐标计算过程是根据像元所在的扫描行号、波段序号、传感器芯片组(Sensor Chip Assemblies,SCA

6、)编号以及像元对应探测器在SCA内的序号，计算该像元对应地面点的经纬度坐标。其计算的完整流程图如图3所示3。下面以OLI传感器为例介绍USGS提出的Landsat 8像元地理坐标计算算法的详细步骤，该算法同样适用于TIRS传感器图像像元地理坐标的计算。2.1 计算像元成像时间像元地理坐标计算过程中需要用到像元成像时刻的卫星姿态和星历数据，而这些数据可根据像元成像时间对已有卫星姿态和星历数据插值计算获得。因此，首先要计算出像元的成像时刻t，具体计算公式为：(1)其中，Tf为该像元所在帧(所有探测器成像一次所产生的数据，多光谱波段每帧对应一行图像数据，全色波段每帧对应两行图像数据)的完成时刻，其值

7、记录在辅助数据中；Ts为探测器成像结束后的采样稳定时间；Ti为探测器曝光时间，Nl为像元所在图像的行号，n为每帧图像数据的行数，Nf为像元所在的帧号，Tsample为每行数据的采样时间。2.2 LOS矢量计算像元成像LOS矢量的计算与传感器的成像扫描方式密切相关，对不同的扫描方式，其LOS矢量计算方法不同。LOS矢量最初被定义在传感器坐标系中，其计算方式与探测器在焦平面坐标系的相对位置有关，在焦平面坐标系中定义了各波段探测器相对于光学轴的偏移量，并最终用于产生单个探测器成像的LOS矢量。目前，全球范围内陆地卫星的成像方式主要有双向扫描式(Landsat 17等)和推扫式(SPOT15/Land

8、sat 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波段分布在冷焦平面上。每个探测器对应的LOS矢量与定义在传感器坐标系下的观测角度有关4，观测角度可分为平行扫描方向的角度(along_angle)和垂直扫描方向

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

10、其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

11、)和跨轨方向观测角度(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)。2.2.3 Landsat 8卫星LOS矢量计算La

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

13、94个主探测器，对于全色波段每个SCA中含有奇偶两排探测器，每排中有988个主探测器5，如图6所示6。鉴于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_shifty=b0

14、+b1*n+b2*(1.5*n2-0.5)+b3*(n*(2.5*n2-1.5)+ifov_y*y_shiftz=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)数

15、据产品，并以SCA为单位计算出已配准控制点像元在L1T图像和标准参考图像(高分辨率数字正射影像镶嵌图或比OLI更高分辨率传感器的几何校正数据产品)中的偏差，将偏差加权最小二乘拟合，计算出勒让德系数的修正量，进而修正已有系数7。采用这种方法可以根据最新卫星数据不断地对传感器进行重新几何定标，从而有效保证了数据处理的几何精度。2.3 卫星姿态计算卫星的姿态运动是卫星绕自身质心的转动运动8，在本文介绍的算法中，使用姿态数据构造出从导航坐标系到轨道坐标系的转换矩阵。Landsat 8卫星用roll、pitch、yaw分别表示卫星姿态中的翻滚角、俯仰角和偏航角。姿态数据采集频率为50Hz，可以根据已有姿

16、态数据线性插值得到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矢量最初被定义在传感器坐标系，因此需要被转换到导航参考坐标系，然

17、后使用姿态数据将LOS转换到轨道坐标系，最后再被转换到ECEF坐标系下，具体转换过程如下所示：(8)其中，Ms2a为传感器坐标系到导航坐标系的转换矩阵，其值由卫星发射前的几何定标数据确定；Ma2o为导航坐标系到轨道坐标系的转换矩阵，其值由卫星姿态数据构造出；Mo2e为轨道坐标系到ECEF坐标系的转换矩阵，其值由卫星星历数据构造出；为LOS矢量在ECEF坐标系下的值。2.6 卫星位置调整在坐标计算过程中最终使用的是传感器的位置信息，而星历数据中记录的是卫星质心的位置信息，因此需要根据卫星质心相对于传感器的偏置调整卫星位置。首先将质心相对于传感器的位置矢量经过坐标变换转换成ECEF坐标系下的位置矢

18、量然后将该位置矢量与原卫星位置矢量相加得到卫星新位置矢量具体计算过程为：(9)(10)2.7 地面点地心经纬度计算像元成像点的坐标即为LOS与地球相交地面点的坐标，如图7所示，其运算关系如式(11)所示：(11)其中，分别为卫星和地面点在ECEF坐标系下的位置矢量，m为矢量的模，该值可根据矢量与矢量的夹角和地球半径利用三角余弦定理求得。同时由于卫星和地球之间存在相对速度，因此会带来视扰差，所以需要对LOS矢量进行校正以消除误差，该过程被首次引入到Landsat系列卫星的像元地理坐标的计算中。具体校正方法为：(12)(13)其中，e为地球自转角速度，分别为卫星和地面点在ECEF坐标系下的速度矢量

19、，vlight为光速，为调整之后的LOS矢量。最后归一化将单位矢量重新记为并带入式(11)重新计算地面点坐标。将表示为则地面点的地心纬度和地心经度可分别表示为：(14)(15)2.8 光速校正卫星高度大概为705km，地面成像点的反射光到达卫星传感器需要一定的时间，这个时间延迟会对计算地面点的坐标带来一定误差，因此需要进行校正，以消除其影响，其校正过程为：(16)(17)其中，为反射光由地面点到达传感器时间内地球转动的角度，为校正后地面点在ECEF坐标系下的位置矢量。2.9 地面点大地纬度计算成像地面点的地心纬度是指该点和地心连线与赤道面的夹角，而大地纬度是指通过该点椭球面的法线与赤道面的夹角

20、。由于二者之间不存在简单的显式几何转换关系，所以转换过程采用迭代计算的方式进行。首先根据该点地心纬度计算该点处地球半径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就是转换

21、后地面点的大地纬度，至此成像地面点的大地经纬度坐标计算完成。3 实验结果为了说明本文所介绍算法的实用性，使用USGS发布的GLS2005控制点库作为参考依据，对算法的处理结果进行分析。该库是在Landsat系列卫星影像数据的基础上提取建立的9，其中80%的控制点定位精度优于30m10。同时，该库也作为Landsat8数据几何精校正阶段的控制点来源。本文的实验数据采用USGS发布的一景Landsat 8模拟数据，该景数据在WRS-2坐标系中的Path/Row为024/038。现从GLS2005控制点库中选取均匀分布在该景数据范围内的10个点作为参考，其分布位置如图8所示。同时从该景数据中分别找出

22、每个参考点对应的扫描行号、波段号、SCA编号和像元序号，使用本文介绍的算法计算出其坐标，并与其对应参考点的坐标进行对比分析，结果如表2所示。从表2可以看出，与GLS2005控制点库中的数据相比，使用本文介绍算法计算成像地面点的经度误差和纬度误差换算成地面距离均在亚像素级，可以满足Landsat 8卫星数据预处理阶段的精度要求。4 结束语本文在充分了解Landsat 7卫星和SPOT5卫星像元地理坐标计算算法的基础上，对USGS提出的Landsat 8像元地理坐标计算算法进行了系统性的整理和研究，并比较了由于不同传感器间扫描方式的差异所带来的计算方法差异，另外在该算法中USGS还考虑到了传感器相

23、对于卫星质心的位置偏移、速度像差、光程时间延迟等，并对这些因素所带来的影响予以校正和消除，使运算流程更科学，计算精度更高，最后以一景Landsat 8的模拟数据为样本，验证了该算法的有效性。本文所介绍的算法在陆地卫星像元地理坐标计算方面具有典型性，对其他型号的陆地卫星(特别是推扫式传感器)具有重要的参考意义。参考文献：1 IRONS J R,DWYER J L,BARSI J A.The next Landsat satellite:The Landsat data continuity missionJ.Remote Sensing of Environment,2012(122):11-2

24、1.2 RIAZANOFF S,GLEYZES J P,ZOBLER D,et al.SPOT 123-4-5 Geometry HandbookR.GAEL-P135-DOC-001,2004.3 USGS.LDCM CAL/VAL ALGORITHM DESCRIPTION DOCUMENTG.2013.4 USGS EROS Data Center.Landsat 7 Image Assessment System (IAS) Geometric Alogrithm Theoretical Basis DocumentG.2011.5 REINKE N,ZAHN S,HOLLAREN D

25、,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 missio

26、n (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.

27、9 USGS Publications Warehouse.Extraction of GCP chips from GeoCover using Modified Moravec Interest Operator (MMIO) algorithmG.2006.10 USGS EROS,GLCF,NASA Headquarters.Assessment of the NASA-USGS Global Land Survey(GLS) DatabasesG.2012.doi:10.3969/j.issn.1000-3177.2013.05.010E-mail:zhanghao_gisE-mail:xxfeng【文献来源】