CN101852623B - 一种卫星光学遥感相机内方元素在轨检校方法 - Google Patents

Abstract

本发明涉及一种卫星光学遥感相机内方元素在轨检校方法，该方法根据通常的光学镜头畸变模型建立了3阶内方元素模型，并基于相对角度误差受外方元素误差影响小的原理，从视向量之间的相对角度关系中提取了内方元素的模型参数，并得到拟合残差，通过调整相机焦距f，使得拟合残差最小，完成在轨检校，本发明方法实现了内外方元素解耦，并通过在轨检校光学遥感相机的几何参数，提高了遥感卫星图像的内部几何精度，减小了光学畸变误差。

Description

一种卫星光学遥感相机内方元素在轨检校方法

技术领域

本发明属于遥感卫星在轨几何检校领域，特别是涉及一种卫星光学遥感相机内方元素在轨检校方法。

背景技术

随着国产遥感卫星空间分辨率的提高，遥感影像的内部几何精度成为一项重要的考察指标。线阵推扫式光学遥感相机是一种主要的光学遥感卫星载荷，其内方位元素检校能够提取光学畸变模型参数，减小畸变，提高影像内部几何精度。

目前国内卫星遥感领域，在光学遥感相机的内方元素几何检校方面，除了在三线阵立体成像相机的在轨动态检校方面有一定的研究外，其他还有对于面阵相机地面检校的研究，其中三线阵相机动态检校的重点是研究三线阵相机的主点和焦距及三个相机之间的几何关系；面阵相机地面检校则要借助地面测试设备如转台和光源发生器，并校正的是面阵模型，但是对于单线阵推扫式相机在轨检校方面则缺少相应的研究。

发明内容

本发明的目的在于克服现有技术的上述不足，提供一种卫星光学遥感相机内方元素在轨检校方法，该方法通过在轨检校光学遥感相机的几何参数，提高了遥感卫星图像的内部几何精度，减小了光学畸变误差。

本发明的上述目的是通过如下技术方案予以实现的：

一种卫星光学遥感相机内方元素在轨检校方法，包括如下步骤：

建立OXYZ三维坐标系，其中相机在坐标系中的位置为O₁，地球中心在坐标系中的位置为O₂，坐标系中X向为线阵推扫方向，Y向为垂轨方向，Z向为遥感相机视轴，设地面某控制点P，向量OP与地球相交于P点，OP与焦平面相交于P₂点，由于光学部件设计与加工缺陷的存在，视向量O₁P在焦平面的实际成像位置为P₁点，设P₂点的Y坐标为y₂，P₁点的Y坐标为y₁，则根据径向畸变模型公式得出y₂与y₁的关系式：

y₂＝c₀+c₁y₁+c₂y₁ ²+c₃y₁ ³

计算模型参数c₀、c₁、c₂和c₃，将计算得到的模型参数c₀、c₁、c₂和c₃代入上式，得到P₂点修正后的Y坐标y₂′，进一步得到拟合残差Δy₂＝y′₂-y₂，调整相机焦距f，使得拟合残差Δy₂最小，完成在轨检校，其中计算模型参数c₀、c₁、c₂和c₃的步骤如下：

(1)选择几何检校区，卫星在几何检校区成像后得到1级图像，同时选择与1级图像具有相同分辨率或更高分辨率的正射影像和DEM高程图作为参考图，在1级图像上选取k个控制点，同时在所述参考图上选取k个同名点，记录1级图像上k个控制点的坐标(m，n)，m表示列号，n表示行号，记录参考图上k个同名点的经纬度坐标和高度值(Lon，Lat，h)，Lon表示经度，Lat表示纬度，h表示高度，其中k为正整数；

(2)根据列号m，计算OXYZ三维坐标系下控制点的y坐标，

y＝c(m-M/2)，其中M表示一行的像元个数，c表示像元尺寸，

得到控制点线阵向量w₀(0，y，f)，进一步归一化得到控制点线阵向量w_i；

(3)根据行号n，计算控制点所在行的成像时刻秒计数t，

t＝n·d-t₀，其中d为积分时间，t₀为1级图像第一行秒计数，

根据秒计数t，利用卫星下传的GPS数据，计算成像时刻卫星在地固系的位置坐标向量P_t，并假设卫星和相机的坐标重合，根据参考图上同名点的经纬度坐标和高度值(Lon，Lat，h)，计算得到地固系下控制点的坐标向量P_j，控制点坐标向量P_j减去卫星位置坐标向量P_t得到控制点视向量v₀，进一步归一化得到控制点视向量v_i；

(4)在k个控制点的线阵向量w_i中选择最接近相机视轴的向量w_ij，计算向量w_ij与视轴的夹角α：

y₀表示离相机视轴最近的控制点的y坐标；

分别以最接近相机视轴的向量w_ij、v_ij为基准，计算其它k-1个向量与所述向量w_ij、v_ij的夹角，得到k-1个控制点线阵向量间夹角Aw_i和k-1个控制点视向量间夹角Av_i，每个夹角都加上一个夹角α，得到Aw’_i和Av’_l；

本步骤中以最接近相机视轴的向量w_ij、v_ij为基准，并假设该向量的畸变可以忽略，得到其它向量与基准向量w_ij、v_ij之间的夹角后，两组夹角数据的差异反映了光学几何畸变的特征，这是本方法中能够解耦卫星外方位元素影响的关键。

(5)计算线阵坐标Bw_i和Bv_i，计算公式为：

Bw_i＝f*tan(Aw′_i)

Bv_i＝f*tan(Av′_i)，f为相机焦距；

线阵坐标Bw_i和Bv_i分别对应公式y₂＝c₀+c₁y₁+c₂y₁ ²+c₃y₁ ³中的y₁和y₂，将Bw_i和Bv_i代入公式y₂＝c₀+c₁y₁+c₂y₁ ²+c₃y₁ ³，用最小二乘法拟合得到模型参数c₀、c₁、c₂和c₃。

步骤(5)中的公式是将两组夹角数据转换为线阵Y坐标方向的数值，使得到的畸变模型系数单位和Y坐标单位一致，方便应用该模型；对于直接利用角度进行几何校正的程序，则可以直接进行角度畸变计算c₀、c₁、c₂和c₃，而模型中的y₁和y₂也对应于角度值。

在上述卫星光学遥感相机内方元素在轨检校方法中，步骤(1)中几何检校区选择人工地物目标丰富，道路交通发达的地区。

在上述卫星光学遥感相机内方元素在轨检校方法中，步骤(1)中在1级图像上选取沿Y方向均匀分布，沿X方向较窄区域中的k个控制点，并且20＜k＜40。

本发明与现有技术相比具有如下优点：

(1)本发明在轨检校方法根据通常的光学镜头畸变模型建立了3阶内方元素模型，并基于相对角度误差受外方元素误差影响小的原理，从视向量之间的相对角度关系中提取了内方元素的模型参数，实现了内外方元素解耦；

(2)本发明在轨检校方法通过对光学高阶畸变建模，并通过最小二乘法拟合得到模型参数，从而得到修正后的线阵坐标，并进一步得到拟合残差，通过调整焦距，使得拟合残差最小，完成几何检校，本发明方法能够达到较高的检校精度；

(3)本发明在轨检校方法采用3阶模型，具有相当的柔韧性，对于不同的焦距误差均能予以吸收，并通过模型参数平衡能够达到总体模型的高精度；

(4)本发明在选择控制点时，选取沿像元排列方向均匀分布的控制点，并且控制点的选择区域为沿X推扫方向的较窄区域，并且越窄约好，目的是使得时间引起的姿态误差达到最小。

附图说明

图1为本发明单线阵推扫相机在轨成像模型；

图2为本发明内方元素几何检校流程图；

图3为本发明内方元素几何检校过程中控制点选取示意图；

图4为本发明实施例中相机畸变拟合前后散点图；

图5为本发明实施例中相机畸变拟合后残差图。

具体实时方式

下面结合附图通过具体实施例对本发明进行进一步详细的描述：

如图1所示为本发明单线阵推扫相机在轨成像模型，建立OXYZ三维坐标系，其中航天相机在坐标系中的位置为O₁，地球中心在坐标系中的位置为O₂，图1中左图表示线阵推扫的情况，右图表示相机视向量和地球相交的情况。相机焦距为f，X为线阵推扫方向，Y为垂轨方向，Z为相机视轴。设地面某控制点P，向量O₁P与地球相交于P点，与焦平面XY相交于P₂点，由于光学部件设计和加工的缺陷，光学相机总是存在一定的几何畸变，使得视向量O₁P在焦平面的实际成像位置为P₁点。

通常光学畸变为径向畸变，提出建立一个3阶多项式畸变模型来表达该径向畸变。

r′＝c₀+c₁r+c₂r²+c₃r³        (1)

其中，r为畸变后实际成像距离即OP₁，r′为共线径向距离OP₂。

地球中心为O₂，卫星和地球的距离为O₁O₂。

向量O₁P₂在OXYZ相机坐标系下构成视向量w(0，y₂，f)，P₂点的Y坐标为y₂。向量O₁P₁在OXYZ相机坐标系下构成视向量v(0，y₁，f)，P₁点的Y坐标为y₁。

根据公式(1)的径向畸变模型，y₁和y₂的关系为：

y₂＝c₀+c₁y₁+c₂y₁ ²+c₃y₁ ³            (2)

计算模型参数c₀、c₁、c₂和c₃，将计算得到的模型参数c₀、c₁、c₂和c₃代入公式(2)，得到P₂点修正后的Y坐标y₂′，进一步得到拟合残差Δy₂＝y′₂-y₂，调整相机焦距f，使得拟合残差Δy₂最小，完成在轨检校。

其中计算模型参数c₀、c₁、c₂和c₃的步骤如下：

步骤一：控制点数据采集。

选择一景人工地物目标丰富，道路交通发达的地区作为几何检校区，卫星在该区域成像后得到1级图像，同时选择与1级图像同分辨率或更高分辨率的正射影像和DEM高程图作为参考图。在经过了辐射校正的1级图像上选取控制点，同时在参考图上选取同名点。记录1级图像上控制点坐标(m，n)，m表示列号，n表示行号，记录参考图上的经纬度坐标和高度值(Lon，Lat，h)，Lon表示经度，Lat表示纬度，h表示高度。

如图3所示为本发明内方元素几何检校过程中控制点选取示意图，每一个“+”表示一个控制点，选择清晰交叉点，个数k个，20＜k＜40。控制点沿像元排列方向选择，均匀分布；控制点的选择区域为沿X推扫方向的较窄区域，并且越窄约好，使得时间引起的姿态误差达到最小。

步骤二：数据处理

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

y＝c(m-M/2)        (3)

这里，M表示一行的像元个数，c表示像元尺寸。得到控制点线阵向量w₀(0，y，f)，进一步归一化得到控制点线阵向量w_i。如图2所示为本发明内方元素几何检校流程图。

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

t＝n·d-t₀            (4)

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

根据卫星秒计数可以推算出成像时刻的卫星位置在地固系下的向量P_t，并假设卫星质心和相机位置O₁重合。

根据地面控制点的经纬度和高度值(Lon，Lat，h)，计算得到地固系下控制点坐标向量P_j。

控制点坐标向量P_j减去卫星位置向量P_t得到控制点视向量v₀，并对该向量归一化得到控制点视向量v_i。

步骤三：角度计算

上述步骤二分别得到了归一化线阵向量w_i和归一化视向量v_i，对于k个控制点，有线阵向量w_i1～k和视向量v_i1～k。

首先，在一系列线阵向量w_i1～k中选择最接近相机视轴的向量w_ij，计算其与视轴的夹角a，该夹角计算公式为：

$a=\arctan \left(\frac{y_0}{f}\right)---\left(5\right)$

其中y₀表示离相机视轴最近的控制点的y坐标，

分别以最接近相机视轴的的向量w_ij、v_ij为基准，计算其他k-1个向量与所述向量w_ij、v_ij的夹角，得到k-1个控制点线阵向量间夹角Aw_i和k-1个控制点视向量间夹角Av_i，每个夹角都加上一个常值角度α，得到新的Aw’_i和Av’_l。

根据三角关系计算在OXYZ相机坐标系下的各个y坐标，即线阵坐标，得到Bw_i和Bv_i，计算公式为：

Bw_i＝f*tan(Aw′_i)

Bv_i＝f*tan(Av′_i)        (6)

线阵坐标Bw_i和Bv_i分别对应公式y₂＝c₀+c₁y₁+c₂y₁ ²+c₃y₁ ³中的y₁和y₂，将Bw_i和Bv_i代入y₂＝c₀+c₁y₁+c₂y₁ ²+c₃y₁ ³，用最小二乘法拟合得到模型参数c₀、c₁、c₂和c₃。

基于模型参数c₀、c₁、c₂和c₃，可以得到修正后的P₂点的Y坐标y₂′，进一步得到拟合残差Δy₂＝y′₂-y₂。调整相机焦距f，使得拟合残差Δy₂最小，完成在轨检校。

下面列举一个具体的实施例

以我国环境减灾卫星1B的CCD2数据为参考，抽取2009年8月29日位于安徽和江苏北部地区的一景遥感影像，进行试验计算。相机CCD像元尺寸0.065mm，对应地面30米。影像参考图选择landsat的ETM正射影像，高程数据为STRM90。控制点共有21个，平面精度约为50米，高程精度约为5米。

如图4所示为本发明实施例相机畸变拟合前后散点图，由图4可知拟合前畸变误差非常明显，其中‘☆’表示控制点畸变误差坐标，‘+’表示控制点误差拟合修正坐标，本实施例分别选择了5个焦距值，得到5组结果，见表1为畸变拟合结果图，图4的相机畸变拟合前后散点图为选取表1中序号1数据得到的结果图。

表1

序号	焦距/mm	C0	C1	C2	C3	残差/mm
							1.	140.8	0.00095899	0.99776	1.7785E-5	2.0774E-5	0.015083
2.	140.9	0.0023983	0.99847	1.7874E-5	2.0789E-5	0.015094
							3.	141.0	0.0038375	0.99918	1.7964E-5	2.0804E-5	0.015104
4.	141.1	0.0052768	0.99989	1.8053E-5	2.0819E-5	0.015115
							5.	141.2	0.006716	1.0006	1.8142E-5	2.0833E-5	0.015126

从结果可以看出，3阶畸变模型得到的残差大约在0.015mm，等于2.3个像元，残差分布见图5，说明残差已经没有明显规律。考虑参考点精度也大约在2个像元左右，因此说明本发明在轨检校方法达到了预期精度，同时，本发明方法采用的3阶模型，具有相当的柔韧性，对于不同的焦距差异，均能予以吸收，图5相机畸变拟合后各个控制点的残差图。

以上所述，仅为本发明最佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

Claims (3)

1.一种卫星光学遥感相机内方元素在轨检校方法，其特征在于包括如下步骤：

建立OXYZ三维坐标系，其中相机在坐标系中的位置为O₁，地球中心在坐标系中的位置为O₂，坐标系中X向为线阵推扫方向，Y向为垂轨方向，Z向为遥感相机视轴，设地面某控制点P，向量O₁P与地球相交于P点，O₁P与焦平面相交于P₂点，由于光学部件设计与加工缺陷的存在，视向量O₁P在焦平面的实际成像位置为P₁点，设P₂点的Y坐标为y₂，P₁点的Y坐标为y₁，则根据径向畸变模型公式得出y₂与y₁的关系式：

y₂＝c₀+c₁y₁+c₂y₁ ²+c₃y₁ ³

计算模型参数c₀、c₁、c₂和c₃，将计算得到的模型参数c₀、c₁、c₂和c₃代入上式，得到P₂点修正后的Y坐标y₂′，进一步得到拟合残差Δy₂＝y′₂-y₂，调整相机焦距f，使得拟合残差Δy₂最小，完成在轨检校，其中计算模型参数c₀、c₁、c₂和c₃的步骤如下：

(1)选择几何检校区，卫星在几何检校区成像后得到1级图像，同时选择与1级图像具有相同分辨率或更高分辨率的正射影像和DEM高程图作为参考图，在1级图像上选取k个控制点，同时在所述参考图上选取k个同名点，记录1级图像上k个控制点的坐标(m，n)，m表示列号，n表示行号，记录参考图上k个同名点的经纬度坐标和高度值(Lon，Lat，h)，Lon表示经度，Lat表示纬度，h表示高度，其中k为正整数，且20＜k＜40；

(2)根据列号m，计算OXYZ三维坐标系下控制点的y坐标，

y＝c(m-M/2)，其中M表示一行的像元个数，c表示像元尺寸，

得到控制点线阵向量w₀(0，y，f)，进一步归一化得到控制点线阵向量w_i；

(3)根据行号n，计算控制点所在行的成像时刻秒计数t， 

t＝n·d-t₀，其中d为积分时间，t₀为1级图像第一行秒计数，

根据秒计数t，利用卫星下传的GPS数据，计算成像时刻卫星在地固系的位置坐标向量P_t，并假设卫星和相机的坐标重合，根据参考图上同名点的经纬度坐标和高度值(Lon，Lat，h)，计算得到地固系下控制点的坐标向量P_j，控制点坐标向量P_j减去卫星位置坐标向量P_t得到控制点视向量v₀，进一步归一化得到控制点视向量v_i；

(4)在k个控制点的线阵向量w_i中选择最接近相机视轴的向量w_ij，计算向量w_ij与视轴的夹角a：

y₀表示离相机视轴最近的控制点的y坐标；

分别以最接近相机视轴的向量w_ij、v_ij为基准，计算其它k-1个向量与所述向量w_ij、v_ij的夹角，得到k-1个控制点线阵向量间夹角Aw_i和k-1个控制点视向量间夹角Av_i，每个夹角都加上一个夹角a，得到Aw’_i和Av’_i；

(5)计算线阵坐标Bw_i和Bv_i，计算公式为：

Bw_i＝f*tan(Aw′_i)

Bv_i＝f*tan(Av′_i)，f为相机焦距；

线阵坐标Bw_i和Bv_i分别对应公式y₂＝c₀+c₁y₁+c₂y₁ ²+c₃y₁ ³中的y₁和y₂，将Bw_i和Bv_i代入公式y₂＝c₀+c₁y₁+c₂y₁ ²+c₃y₁ ³，用最小二乘法拟合得到模型参数c₀、c₁、c₂和c₃。

2.根据权利要求1所述的一种卫星光学遥感相机内方元素在轨检校方法，其特征在于：所述步骤(1)中几何检校区选择人工地物目标丰富，道路交通发达的地区。

3.根据权利要求1所述的一种卫星光学遥感相机内方元素在轨检校方法，其特征在于：所述步骤(1)中在1级图像上选取沿Y方向均匀分布，沿X方向较窄区域中的k个控制点，并且20＜k＜40。 

Images