CN103674063B - 一种光学遥感相机在轨几何定标方法 - Google Patents

Abstract

一种光学遥感相机在轨几何定标方法，针对国产卫星平台和载荷设计特点，构建严格成像模型，采用内外方位元素联合建模的方法对影像进行几何定标处理。本方法基于严格成像模型进行几何定标处理，模型参数具有明确的物理含义；2）内外方位元素联合建模，同时实现内外方位元素的解耦；3）模型可以用于多线阵相机和多相机联合定标，具有很强的普适性。4）充分考虑成像过程中导致影像畸变的各种因素，能够实现高精度。

Description

一种光学遥感相机在轨几何定标方法

技术领域

本发明属于国产遥感卫星在轨几何定标领域，适用于线阵推扫式光学CCD遥感相机。

背景技术

随着我国航天技术的发展,国产遥感卫星空间分辨率逐步提高,遥感影像的实际定位精度成为一项重要的考察指标。线阵推扫式光学遥感相机是一种主要的光学遥感卫星载荷，对其内外方位元素进行在轨定标，能够减小影像畸变，提高影像对地定位精度。

从目前的发表文章和公开资料来看，国内遥感相机检较领域，在光学遥感相机的几何检校方面，主要研究工作集中在地面和航空摄影测量相机技术参数的标定，且大多是针对传统框幅式相机展开深入研究；而少量的对于航天线阵相机在轨几何定标的研究，则是利用一种数学模型去拟合内外方位元素的变化产生的影像畸变，无法对成像过程中各物理参数进行标定，定标精度均不能达到预定要求,与国际水平有一个量级的差距。

在国外，则已发射多颗高分辨率光学遥感卫星，如IKONOS、Quickbird、WordView-1/2、Geo-eye以及Pleiades等，并形成了相应的在轨几何定标技术方案，但是公开资料中只有简要描述，没有公开技术细节。

发明内容

本发明的技术解决问题是：通过对国产光学卫星遥感相机进行在轨定标，可以提高遥感卫星影像的几何精度，从而提高影像对地定位精度。

本发明的技术解决方案是：一种光学遥感相机在轨几何定标方法，步骤如下：

（1）卫星在选择的区域成像后得到1级图像并进行辐射校正，选择同分辨率或更高分辨率的正射影像和DEM高程图作为参考图；在经过了辐射校正的1级图像上和参考图上选取同名点，1级图像上的同名点坐标记为（m，n），参考图上的同名点的经纬度和高度值记为(Lon,Lat,h)，所述的同名点分为两部分，一部分为供计算用的控制点，一部分为验证用的检查点；

（2）根据控制点的坐标m以及CCD线阵的安装位置，确定控制点在焦平面坐标系下的坐标（x，y）；根据控制点的坐标n，计算控制点所在行的成像时刻秒计数t，进而确定卫星位置的测量坐标X_GPS(t),Y_GPS(t),Z_GPS(t)和姿态的测量角度φ(t),κ(t)；同时根据参考图上的同名点的经纬度和高度值(Lon,Lat,h)，计算同名点在地球固定坐标系下的坐标(X_obj,Y_obj,Z_obj)；

（3）利用步骤（2）的结果，结合星地接口文件，构建严格成像模型；

（4）将定标参数代入严格成像模型，得到图像定标模型，初始化定标参数；

（5）利用谱修正法对图像定标模型进行平差迭代，直至解稳定为止；

（6）利用步骤（5）解稳定后得到的定标参数确定1级图像上检查点在焦平面坐标系下的实际坐标（x’，y’），并利用传感器坐标系到地球固定坐标系的旋转矩阵，得到参考图像上对应检查点在焦平面坐标系下的坐标；

（7）将步骤（6）得到的1级图像上检查点在焦平面坐标系下的实际坐标与参考图上对应检查点在焦平面坐标系下的坐标求差，该差记为定位偏差；

（8）调整定标参数的初始值及迭代阈值，重复步骤（5）-（7），使得定位偏差最小，完成几何定标。

所述的图像定标模型为

$\left\{\begin{array}{c}x=\mathrm{\Δx}-f\frac{r_{11}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{21}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{31}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}{r_{13}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{23}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{33}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}\\ y=\mathrm{\Δy}-f\frac{r_{12}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{22}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{32}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}{r_{13}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{23}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{33}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}\end{array}\right.$

其中， $\left\{\begin{array}{c}\mathrm{\Δx}={\mathrm{\Δx}}_p+(k_1{r}^2+k_2{r}^4+k_3{r}^6)\stackrel{\‾}{x}+p_1(r^2+2{\stackrel{\‾}{x}}^2)+2p_2\stackrel{\‾}{\mathrm{xy}}+\stackrel{\‾}{y}\sin \mathrm{\θ}+\stackrel{\‾}{y}{r}^{2}b\\ \mathrm{\Δy}={\mathrm{\Δy}}_p+(k_1{r}^2+k_2{r}^4+k_3{r}^6)\stackrel{\‾}{x}+2p_1\stackrel{\‾}{\mathrm{xy}}+p_2(r^2+2{\stackrel{\‾}{y}}^2)+s\stackrel{\‾}{y}\end{array}\right.$

x,y为像元在焦平面坐标系下的坐标，f为相机主距，x_p,y_p为相机主点位置，Δx,Δy为焦平面内某一像元实际位置与理想位置之间的偏差，Δx_p,Δy_p为线阵CCD中心像点相对像主点的偏移量，b为线阵CCD的曲率系数，s为相机的比例参数，θ为线阵CCD与焦平面坐标系Y轴的夹角，k₁,k₂,k₃为径向畸变系数，p₁,p₂为偏心畸变系数，r²=(x-x_p)²+(y-y_p)²; $\stackrel{\‾}{x}=\left({x-x}_p\right);\stackrel{\‾}{y}=\left({y-y}_p\right);$ $\left[\begin{array}{c}{r}_{11},r_{12},r_{13}\\ r_{21},r_{22},r_{23}\\ r_{31},r_{31},r_{33}\end{array}\right]$ 为传感器坐标系到地球固定坐标系的旋转矩阵。

本发明与现有技术相比有益效果为：

1）基于严格成像模型进行几何定标处理，模型参数具有明确的物理含义；2）内外方位元素联合建模，同时实现内外方位元素的解耦；3）模型可以用于多线阵相机和多相机联合定标，具有很强的普适性。4）充分考虑成像过程中导致影像畸变的各种因素，能够实现高精度。

附图说明

图1为线阵推扫相机在轨成像模型示意图；

图2为线阵推扫相机焦平面安装示意图；

图3为本发明方法流程图；

图4为本发明控制点选取示意图；

图5a、5b分别为定标前控制点处的平面定位偏差和高程定位偏差；

图6a、6b分别为定标后控制点处的平面定位偏差和高程定位偏差。

具体实施方式

下面结合附图及实例对本发明做详细说明，为了更清楚的了解本发明，首先对方案中涉及的坐标系进行说明。

焦平面坐标系：以影像的左上角为坐标原点，沿着扫描线方向为Y轴，沿着卫星运行方向为X轴；

传感器坐标系：原点在线阵投影中心，Y轴沿着扫描方向，X轴沿着卫星飞行方向，Z轴按照右手法确定；

地球惯性坐标系：原点为地球质心，Z轴指向天球的北极，X轴指向春分点，Y轴按照右手法确定；

地球固定坐标系：原点为地球质心，Z轴指向地球的北极，X轴指向格林尼治子午线与地球赤道焦点，Y轴按照右手法确定。

遥感光学卫星在轨几何定标，是数据处理中的重要环节，其难点在于定标模型的构建和模型解算。本方案使用内外方位联合建模的方法，构建定标模型，利用谱修正迭代方法求解模型参数，从而实现内外方位元素的标定。

线阵推扫式遥感光学相机的成像模型如图1所示，卫星平台的姿态测量主要靠星载GPS和星敏感器等部件，Y方向为卫星飞行方向，X方向为相机线阵扫描方向（垂轨方向），O为相机投影中心，视向量OA与地面交于A，与焦平面交于a。由于卫星平台外方位元素的测量误差，以及相机主点、焦距、CCD阵列的安装偏差等内方位元素的误差，导致OA与地面实际相交于A′，与焦平面交于a′，AA′即相机对地定位误差，为了提高对地定位精度，需要减小AA′的距离。

从以上论述不难看出，为了达到高精度的定标结果，定标模型需要在严格成像模型的基础上，统筹考虑外方位元素和内方位元素的影响。

外方位元素：

采用线性***误差校正模型对卫星平台外方位元素进行建模：

卫星平台位置***误差校正模型，如（1）所示：

其中，X_sat(t),Y_sat(t),Z_sat(t)为t时刻卫星平台的实际坐标，X_GPS(t),Y_GPS(t),Z_GPS(t)为t时刻卫星平台测量坐标；φ_sat(t),κ_sat(t)为t时刻卫星平台的实际姿态，φ(t),κ(t)为t时刻卫星平台姿态的测量角度；为由于最初安置不当形成的星敏感器与陀螺等的姿态测量常差。为随时间累积的姿态飘移误差，X_off,Y_off,Z_off是卫星位置偏移参数。

内方位元素：

如图2所示为线阵推扫相机焦平面安装示意图，地物反射的光线经过镜头和光路传播***投射到焦平面上，被安装在焦平面的CCD感知。由于线阵安装和镜头参数的变化，使得实际的影像发生畸变。本方案对镜头主点、传感器主距、镜头的偏心畸变、径向畸变、CCD像元的尺寸变化、CCD的安装位置偏差（包括旋转和偏移）、CCD的弯曲这些内方位元素进行标定，模型如（2）所示：

$\left\{\begin{array}{c}\mathrm{\Δx}={\mathrm{\Δx}}_p+(k_1{r}^2+k_2{r}^4+k_3{r}^6)\stackrel{\‾}{x}+p_1(r^2+2{\stackrel{\‾}{x}}^2)+2p_2\stackrel{\‾}{\mathrm{xy}}+\stackrel{\‾}{y}\sin \mathrm{\θ}+\stackrel{\‾}{y}{r}^{2}b\\ \mathrm{\Δy}={\mathrm{\Δy}}_p+(k_1{r}^2+k_2{r}^4+k_3{r}^6)\stackrel{\‾}{x}+2p_1\stackrel{\‾}{\mathrm{xy}}+p_2(r^2+2{\stackrel{\‾}{y}}^2)+s\stackrel{\‾}{y}\end{array}\right.---\left(2\right)$

其中，(x,y)为像元在焦平面内的坐标，f为相机主距，x_p,y_p为相机主点位置，Δx,Δy为焦平面内某一像元实际位置与理想位置之间的偏差，Δx_p,Δy_p为线阵CCD中心像点相对像主点的偏移量，b为线阵CCD的曲率系数，s为相机的比例参数，θ为线阵CCD与焦平面坐标系Y轴的夹角，k₁,k₂,k₃为径向畸变系数，p₁,p₂为偏心畸变系数，r²=(x-x_p)²+(y-y_p)²;传感器主距的变化和像元尺寸的变化具有同一效应，统一用相机变化因子s表示，CCD的中心偏移与像主点的偏移有相同的相应，统一用Δx_p,Δy_p表示。

严格成像模型：

设[x,y,-f]^T为像元在传感器坐标系下的坐标，[X_sat,Y_sat,Z_sat]^T为卫星平台在地球固定坐标系下的坐标，[X_obj,Y_obj,Z_obj]^T为地面控制点在地球固定坐标系下的坐标， $r_{\mathrm{sat}}^{\mathrm{cis}}=\left[\begin{array}{c}{r}_{11},r_{12},r_{13}\\ r_{21},r_{22},r_{23}\\ r_{31},r_{31},r_{33}\end{array}\right]$ 为传感器坐标到地球惯性坐标的旋转矩阵，κ为姿态角，则

为地球惯性坐标系到地球固定坐标系的旋转矩阵，为传感器坐标系坐标到地球固定坐标系的旋转矩阵，由共线方程可知：

$\left\{\begin{array}{c}x=-f\frac{r_{11}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{21}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{31}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}{r_{13}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{23}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{33}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}\\ y=-f\frac{r_{12}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{22}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{32}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}{r_{13}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{23}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{33}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}\end{array}\right.---\left(3\right)$

如图3所示，具体步骤如下：

（1）选择一景人工地物目标丰富，道路交通发达的地区作为几何检校区，控制点的选择尽量均匀地覆盖整个图幅区域,使得定标参数更能代表整个区域的畸变特征，并且防止定标模型方程出现病态，使得参数解算稳定。卫星在该区域成像后得到1级图像，同时选择同分辨率或更高分辨率的正射影像和DEM高程图作为参考图。在经过了辐射校正的1级图像上和参考图上选取同名点，1级图像上的同名点坐标记为(m,n)，m表示列号，n表示行号，参考图上的同名点的经纬度和高度值记为(Lon,Lat,h)，Lon表示经度，Lat表示纬度，h表示高度；所述的同名点分为两部分，一部分为供计算用的控制点，一部分为验证用的检查点。

控制点选择方式如图4所示，每一个“+”表示一个控制点，选择清晰交叉点，个数k个，15<k<45。控制点在整个图幅区域内均匀分布,使得定标参数更能代表整个区域的畸变特征，并且防止定标模型方程出现病态，使得参数解算稳定。检查点的选取从统计学角度来讲一般大于10个点。

（2）数据处理，获取严格成像模型相关参数。

首先，根据列号m，计算相机坐标系下y坐标：

y=c(m-M/2) （5）

这里，M表示一行的像元个数，c表示像元尺寸。得到焦平面坐标系下的y坐标。

根据星地接口文件中说明的CCD线阵安装位置（如图2所示），得到焦平面坐标系下的x坐标。

然后，根据行号n，可以计算控制点所在行的成像时刻秒计数：

t=n·d-t₀ （6）

其中d为积分时间，t₀为该景第一行秒计数。

根据卫星秒计数可以推算出t时刻的测量坐标X_GPS(t),Y_GPS(t),Z_GPS(t)和φ(t)，κ(t)。

根据星地接口文件说明的相机安装位置,结合推算出的姿态φ(t)，κ(t)，可以计算出传感器坐标到地固系坐标的旋转矩阵r。

根据星地接口文件得到相机焦距f。

最后，根据参考图上的同名点经纬度和高度值(Lon,Lat,h)，计算得到同名点在地球固定坐标系下的坐标(X_obj,Y_obj,Z_obj)。

（3）利用步骤（2）的结果，构建严格成像模型；

上一步分别得到了焦平面坐标(x,y)、卫星测量坐标X_GPS(t),Y_GPS(t),Z_GPS(t)、焦平面坐标到地固系坐标的旋转矩阵r、参考图像上同名点在地球固定坐标系下的坐标(X_obj,Y_obj,Z_obj)。

根据卫星成像正向模型，构建严格成像模型，如公式（3）所示。

（4）将定标参数带入严格成像模型，得到图像定标模型，初始化定标参数；

当影像发生畸变时，公式（3）可表示为：

$\left\{\begin{array}{c}x=\mathrm{\&Delta;x}-f\frac{r_{11}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{21}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{31}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}{r_{13}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{23}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{33}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}\\ y=\mathrm{\&Delta;y}-f\frac{r_{12}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{22}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{32}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}{r_{13}(X_{\mathrm{obj}}-X_{\mathrm{sat}})+r_{23}(Y_{\mathrm{obj}}-Y_{\mathrm{sat}})+r_{33}(Z_{\mathrm{obj}}-Z_{\mathrm{sat}})}\end{array}\right.---\left(4\right)$

将公式（1）中当前时刻t对应的X_sat(t),Y_sat(t),Z_sat(t)作为X_sat,Y_sat,Z_sat和公式（2）一起代入公式（4）即得到光学遥感影像的定标模型，解算此模型可以得到内外方位元素的标定结果，完成定标。

（5）利用谱修正法对图像定标模型进行平差迭代，直至解稳定为止；

根据星地接口文件，得到各模型参数的初始值，由于卫星距离地面的距离远远大于地面两个控制点之间的距离，导致控制点间的视向量夹角很小，模型解算时每个控制点的误差方程出现强相关，法方程病态。这里，提出利用谱修正迭代方法，解决该问题：

X_k+1=(A^TA+I)^-1(A^TL+X^k),k=1,2,...... （7）

式中，I为单位矩阵，X^k表示第k次迭代的结果。A为误差方程系数矩阵，L为自由项。

利用上式迭代方法可以有效防止了法方程的病态问题，设定阈值ε，直到|X_k+1k-X_k|<ε为止。

（6）利用步骤（5）解稳定后得到的定标参数确定1级图像中检查点在焦平面坐标系下的实际坐标，并利用传感器坐标系到地球固定坐标系的旋转矩阵，得到参考图像上对应检查点在焦平面坐标系下的坐标；

（7）将步骤（6）得到的图像上检查点在焦平面坐标系下的坐标与参考图上检查点在焦平面坐标系下的坐标求差，该差记为定位偏差；

（8）调整定标参数的初始值及迭代阈值ε，重复步骤（5）-（7），使得定位偏差最小，完成几何定标。

实施例

利用我国国产资源02C多光谱数据对定标模型进行试验。选择2012年9月15日覆盖北京地区的一景遥感影像作为待校正影像，相机CCD像元尺寸为13μm，对应地面10m。以中国资源卫星应用中心生产的正射影像为参考影像，高程数据为SRTM3。控制点共有33个，平面精度为10米，高程精度约为15米。

拟合前畸变误差非常明显，见图5a、5b，6a、6b，箭头方向表示畸变方向，箭头长度表示畸变大小。经本方法定标后，各参数定标结果见表1,定标残差见表2。从定标结果可以看出，定标前图像x方向的畸变均为负值，三个方向定位精度均为百米量级。经过定标后，各方向的畸变残留已经不具有规律性，定标精度优于0.8像素。因此本方法较以前的定标方法达到了高精度定标的目的，且本方法可以应用到多线阵、多相机的联合定标的研究中。

表1定标结果

表2.畸变拟合结果

本发明未详细说明部分属于本领域技术人员公知常识。

Claims (1)

1.一种光学遥感相机在轨几何定标方法，其特征在于步骤如下：

(1)卫星在选择的区域成像后得到1级图像并进行辐射校正，选择同分辨率或更高分辨率的正射影像和DEM高程图作为参考图；在经过了辐射校正的1级图像上和参考图上选取同名点，1级图像上的同名点坐标记为(m，n)，参考图上的同名点的经纬度和高度值记为(Lon,Lat,h)，所述的同名点分为两部分，一部分为供计算用的控制点，一部分为验证用的检查点；

(2)根据控制点的坐标m以及CCD线阵的安装位置，确定控制点在焦平面坐标系下的坐标(x，y)；根据控制点的坐标n，计算控制点所在行的成像时刻秒计数t，进而确定卫星位置的测量坐标X_GPS(t),Y_GPS(t),Z_GPS(t)和姿态的测量角度φ(t),κ(t)；同时根据参考图上的同名点的经纬度和高度值(Lon,Lat,h)，计算同名点在地球固定坐标系下的坐标(X_obj,Y_obj,Z_obj)；

(3)利用步骤(2)的结果，结合星地接口文件，构建严格成像模型；

(4)将定标参数代入严格成像模型，得到图像定标模型，初始化定标参数；所述定标参数包括内方位元素和外方位元素，所述的图像定标模型为

$\left\{\begin{array}{c}x=\mathrm{\&Delta;}x-f\frac{r_{11}(X_{obj}-X_{sat})+r_{21}(Y_{obj}-Y_{sat})+r_{31}(Z_{obj}-Z_{sat})}{r_{13}(X_{obj}-X_{sat})+r_{23}(Y_{obj}-Y_{sat})+r_{33}(Z_{obj}-Z_{sat})}\\ y=\mathrm{\&Delta;}y-f\frac{r_{12}(X_{obj}-X_{sat})+r_{22}(Y_{obj}-Y_{sat})+r_{32}(Z_{obj}-Z_{sat})}{r_{13}(X_{obj}-X_{sat})+r_{23}(Y_{obj}-Y_{sat})+r_{33}(Z_{obj}-Z_{sat})}\end{array}\right.$

其中，

x,y为像元在焦平面坐标系下的坐标，f为相机主距，x_p,y_p为相机主点位置，Δx,Δy为焦平面内某一像元实际位置与理想位置之间的偏差，Δx_p,Δy_p为线阵CCD中心像点相对像主点的偏移量，b为线阵CCD的曲率系数，s为相机的比例参数，θ为线阵CCD与焦平面坐标系Y轴的夹角，k₁,k₂,k₃为径向畸变系数，p₁,p₂为偏心畸变系数，r²＝(x-x_p)²+(y-y_p)²；为传感器坐标系到地球固定坐标系的旋转矩阵；

(5)利用谱修正法对图像定标模型进行平差迭代，直至解稳定为止；

(6)利用步骤(5)解稳定后得到的定标参数确定1级图像上检查点在焦平面坐标系下的实际坐标(x^,，y^,)，并利用传感器坐标系到地球固定坐标系的旋转矩阵，得到参考图像上对应检查点在焦平面坐标系下的坐标；

(7)将步骤(6)得到的1级图像上检查点在焦平面坐标系下的实际坐标与参考图上对应检查点在焦平面坐标系下的坐标求差，该差记为定位偏差；

(8)调整定标参数的初始值及迭代阈值，重复步骤(5)-(7)，使得定位偏差最小，完成几何定标。