一种线阵推扫式相机在轨光学畸变参数标定方法[发明专利]

(12)发明专利申请

(10)申请公布号(10)申请公布号 CN 103673995 A(43)申请公布日 2014.03.26

(21)申请号 201310629647.6(22)申请日 2013.11.29

(71)申请人航天恒星科技有限公司

地址100086 北京市海淀区知春路82号院(72)发明人马灵霞 郝胜勇 王剑 邹同元

丁火平 马狄 郭翠翠 张荞纪强 李宝明(74)专利代理机构中国航天科技专利中心

11009

代理人安丽(51)Int.Cl.

G01C 11/02(2006.01)G03B 43/00(2006.01)

权利要求书3页 说明书9页 附图2页权利要求书3页 说明书9页 附图2页

(54)发明名称

一种线阵推扫式相机在轨光学畸变参数标定方法(57)摘要

一种线阵推扫式相机在轨光学畸变参数标定方法，步骤为：（1）采用基于控制点的自动匹配算法采集得到密集的高精度控制点信息；（2）建立各控制点处严格成像几何模型；（3）计算各控制点处理论指向矢量与实际指向矢量，建立相机内方位元素几何检校数学模型；（4）建立理论指向矢量与图像列号对应的畸变多项式模型；（5）最小二乘法解算误差方程，迭代求解得到相机内部光学畸变标定参数。本发明通过选取大量高精度地面控制点，对卫星在轨相机内部光学系统畸变进行标定，标定结果可用于提升卫星图像无控定位精度和有控定位精度，利用内畸变标定参数进行相机内方位元素修正后，图像内部畸变基本消除。

CN 103673995 ACN 103673995 A

权 利 要 求 书

1/3页

1.一种线阵推扫式相机在轨光学畸变参数标定方法，其特征在于包括下列步骤： 1）在待检校影像I1、控制影像I2上利用控制点自动匹配算法，采集控制点信息，所述的控制点信息为T个匹配控制点对{PixI1i，PixI2i}，其中i=1,2,3,…,T）；采用格网筛选策略，使匹配控制点对在待检校影像区域内均匀分布；记录每一个控制点对{PixI1i，PixI2i}的待检校影像坐标（i,j）以及对应控制影像特征点在WGS84坐标系下的三维位置坐标（lat,lon,height）；

2）获取各控制点处的相机内方位元素以及辅助数据，所述的辅助数据包括卫星的轨道、姿态、行时；建立各控制点处的严格成像模型

其中[XS,YS,ZS]为各控制点对应的

行时时刻下卫星在协议地心坐标系中的位置，即控制点处轨道位置（PX,PY,PZ）；[XG,YG,ZG]为各控制点对应像元所对应的地面目标点在协议地心坐标系中的坐标；psiX、psiY分别为各控制点对应图像的像元主光轴单位矢量与卫星本体坐标系X轴、Y轴的夹角；u为比例因

在卫星发射前由地面测量获取；M1为各控子；M0为卫星本体坐标系相对于相机的安装矩阵，

制点对应的行时时刻下卫星至轨道坐标系旋转矩阵，由星上测量的姿态角构成；M2为各控制点对应的行时时刻下轨道至J2000.0坐标系旋转矩阵；M3为该时刻下J2000.0至WGS84坐标系旋转矩阵；M为卫星本体坐标系相对于相机的安装补偿矩阵；

3）根据步骤2）得到的严格成像模型，计算各控制点对应理论指向矢量，并获得待检校影像点对应在沿轨和垂轨方向理论探元指向角；

31）对任一组控制点，利用控制点地面点三维位置坐标（lat,lon,height），计算控制点处光轴在相机坐标系中指向矢量V，将其进行归一化 处理后作为该控制点处的理论指向矢量V理论；

其中a为地球长半轴，b为地球短半轴；

32）根据步骤31）得到的理论指向矢量，计算获得探元在沿轨和垂轨方向的理论探元指向角

2

CN 103673995 A

权 利 要 求 书

2/3页

4）建立图像列与到理想探元指向角之间的函数模型其中

S为图像列号，A0、A1、A2、A3以及B0、B1、B2、B3为函数

模型系数；对各控制点对构建参数解算误差方程

x=(ATA+E)-1(ATb+x)；其

中

E为单位矩阵；

5）采用最小二乘法解算步骤4）得到的参数解算误差方程；得到内检校参数。

2.根据权利要求1所述的一种线阵推扫式相机在轨光学畸变参数标定方法，其特征在于：步骤1）中采集控制点信息的具体方法为：

11）基于SIFT算法对待检校原始图像I1进行特征点提取，得到M个特征点PixI1i

（i=1,2,3,…,M），M为正整数；记录各特征点的SIFT特征向量；

12）基于SIFT算法对高精度控制影像I2进行特征点提取，得到N个特征点PixI2i

（i=1,2,3,…,N），N为正整数；记录各特征点处的SIFT特征向量；

13）采用欧式距离作为相似性度量准则对两幅影像的特征点进行匹配，得到T个匹配控制点对{PixI1i，PixI2i}（i=1,2,3,…,T）。

3.根据权利要求1所述的一种线阵推扫式相机在轨光学畸变参数标定方法，其特征在于：步骤2）中建立各控制点处的严格成像模型的具体方法为：

21）根据控制点对{PixI1i，PixI2i}的图像坐标（i,j），获取其对应的行时； 22）利用拉格朗日插值算法插值计算控制点对{PixI1i，PixI2i}对应行时的轨道位置

3

CN 103673995 A

权 利 要 求 书

3/3页

（PX,PY,PZ,VX,VY,VZ）；

23）对卫星下传的姿态四元数进行坐标系转换处理，利用拉格朗日算法插值计算控制点对{PixI1i，PixI2i}对应行时时刻处相机相对于轨道坐标系的三轴姿态角（Roll，Pitch，Yaw）；

24）根据步骤21）-步骤23）得到的结果，建立控制点处的严格成像模型，针对线阵推扫成像相机某一时刻获取的图像上的某一点，构建遥感影像的严格成像模型。

4

CN 103673995 A

说 明 书

一种线阵推扫式相机在轨光学畸变参数标定方法

1/9页

技术领域

本发明涉及一种线阵推扫式相机在轨光学畸变参数标定方法，属于遥感卫星图像

处理技术领域，用于遥感卫星搭载的线性推扫式相机的在轨定位精度优化提升和图像内部几何质量的提升优化。

[0001]

背景技术

随着空间技术的不断发展，卫星遥感影像已成为当前地理空间信息服务和应用的

重要数据来源，其广泛的应用前景已引起世界各国的普遍重视。与几年前高分辨率卫星数据的获取能力相比，最近几年世界各国均争相研制自主的高分辨率对地观测系统，遥感影像数据的供应即将要达到“井喷”状态。[0003] 截止2012年底，我国已经发射了环境、资源、遥感等多系列多颗光学遥感卫星，其中已发射的遥感卫星中资源一号02星定位精度为7Km，资源二号03星定位精度200m，资源02C定位精度100m。尽管国产卫星几何定位精度逐步提升，但和国外商业高分辨光学遥感卫星如Geoeye、WorldView、Pleiades等相比，国产卫星在几何定位精度方面尤其是图像的内部几何精度方面还存在很大的差距，无法满足国土资源、测绘等领域正射影像产品的制作要求，严重地制约了我国高分辨率卫星影像的应用性能。[0004] 受发射冲力，在轨环境等因素的影响，航天遥感光学相机的内部几何畸变参数在轨后发生了很大改变，地面测量的相机内部参数在用于图像的几何处理时会导致较大的内部畸变，尤其是远离主点的图像的内部畸变可达到30米量级，严重制约着其定位精度的提高，并影响后续的图像拼接镶嵌等高级应用。

[0002]

发明内容

本发明解决的技术问题是：克服现有技术的不足，提供一种线阵推扫式相机在轨光学畸变参数标定方法，用于提升相机获取影像的内部几何定位精度，达到相机内检校的目的。

[0006] 本发明的技术方案是：一种线阵推扫式相机在轨光学畸变参数标定方法，包括下列步骤：[0007] 1）在待检校影像I1、控制影像I2上利用控制点自动匹配算法，采集控制点信息，所述的控制点信息为T个匹配控制点对{PixI1i，PixI2i}，其中i=1,2,3,…,T）；采用格网筛选策略，使匹配控制点对在待检校影像区域内均匀分布；记录每一个控制点对{PixI1i，PixI2i}的待检校影像坐标（i,j）以及对应控制影像特征点在WGS84坐标系下的三维位置坐标（lat,lon,height）；

[0008] 2）获取各控制点处的相机内方位元素以及辅助数据，所述的辅助数据包括卫星的轨道、姿态、行时；建立各控制点处的严格成像模型

[0005]

5

CN 103673995 A

说 明 书

2/9页

其中[XS,YS,ZS]为各控制点对应的

行时时刻下卫星在协议地心坐标系中的位置，即控制点处轨道位置（PX,PY,PZ）；[XG,YG,ZG]为各控制点对应像元所对应的地面目标点在协议地心坐标系中的坐标；psiX、psiY分别为各控制点对应图像的像元主光轴单位矢量与卫星本体坐标系X轴、Y轴的夹角；u为比例因子；M0为卫星本体坐标系相对于相机的安装矩阵，在卫星发射前由地面测量获取；M1为各控制点对应的行时时刻下卫星至轨道坐标系旋转矩阵，由星上测量的姿态角构成；M2为各控制点对应的行时时刻下轨道至J2000.0坐标系旋转矩阵；M3为该时刻下J2000.0至WGS84坐标系旋转矩阵；M为卫星本体坐标系相对于相机的安装补偿矩阵；[0009] 3）根据步骤2）得到的严格成像模型，计算各控制点对应理论指向矢量，并获得待检校影像点对应在沿轨和垂轨方向理论探元指向角；[0010] 31）对任一组控制点，利用控制点地面点三维位置坐标（lat,lon,height），计算控制点处光轴在相机坐标系中指向矢量V，将其进行归一化处理后作为该控制点处的理论指向矢量V理论；

[0011]

[0012] 其中a为地球长半轴，b为地球短半

轴；

[0013]

32）根据步骤31）得到的理论指向矢量，计算获得探元在沿轨和垂轨方向的理论探

元指向角

[0014] 4）建立图像列与到理想探元指向角之间的函数模型其中

S为图像列号，A0、A1、A2、A3以及B0、B1、B2、B3为函数

模型系数；对各控制点对构建参数解算误差方程

6

CN 103673995 A

说 明 书

3/9页

[0015]

x=(ATA+E)-1(ATb+x)；其中

E为单位矩阵；

5）采用最小二乘法解算步骤4）得到的参数解算误差方程；得到内检校参数。

[0017] 步骤1）中采集控制点信息的具体方法为：[0018] 11）基于SIFT算法对待检校原始图像I1进行特征点提取，得到M个特征点PixI1i（i=1,2,3,…,M），M为正整数；记录各特征点的SIFT特征向量；[0019] 12）基于SIFT算法对高精度控制影像I2进行特征点提取，得到N个特征点PixI2i（i=1,2,3,…,N），N为正整数；记录各特征点处的SIFT特征向量；[0020] 13）采用欧式距离作为相似性度量准则对两幅影像的特征点进行匹配，得到T个匹配控制点对{PixI1i，PixI2i}（i=1,2,3,…,T）；[0021] 步骤2）中建立各控制点处的严格成像模型的具体方法为：[0022] 21）根据控制点对{PixI1i，PixI2i}的图像坐标（i,j），获取其对应的行时；[0023] 22）利用拉格朗日插值算法插值计算控制点对{PixI1i，PixI2i}对应行时的轨道位置（PX,PY,PZ,VX,VY,VZ）；[0024] 23）对卫星下传的姿态四元数进行坐标系转换处理，利用拉格朗日算法插值计算控制点对{PixI1i，PixI2i}对应行时时刻处相机相对于轨道坐标系的三轴姿态角（Roll，Pitch，Yaw）；[0025] 24）根据步骤21）-步骤23）得到的结果，建立控制点处的严格成像模型，针对线阵推扫成像相机某一时刻获取的图像上的某一点，构建遥感影像的严格成像模型。[0026] 本发明与现有技术相比的优点在于：[0027] （1）本发明对遥感卫星线阵推扫式光学相机严格成像模型以及辅助数据特性进

[0016]

7

CN 103673995 A

说 明 书

4/9页

行了深入分析，设计了包含大量均匀分布的影像高精度地面控制点信息获取、相机轨道、姿态、行时等辅助数据优化处理、几何外检校参数的引入、严格几何内检校数学模型建立、建立理想指向矢量对应的探元指向角与图像列号函数关系、探元指向角畸变多项式模型参数解算、最小二乘迭代解算内检校参数等步骤的高精度内参数检校流程，满足卫星图像在轨几何内检校需求，消除图像畸变。[0028] （2）建立几何内检校数学模型时，采用多项式模型拟合相机内部指向角畸变，避免了因内方位参数的相关性导致无法求解。在解算内部畸变参数时采用了基于最小二乘的迭代误差剔除设计，可以进一步剔除粗差控制点，提高了内畸变参数检校的精度。附图说明

图1为本发明的整体流程图。

[0030] 图2为基于待检校影像选取的控制点分布示意图。

[0029]

具体实施方式

[0031] 下面结合附图1、图2对本发明的具体实施方式进行进一步的详细描述：[0032] 1.在待检校影像I1、高精度控制影像I2上利用控制点自动匹配算法采集高精度控制点信息。[0033] （1）基于SIFT算法对待检校原始图像I1进行特征点提取，得到m个特征点PixI1i（i=1,2…,m），记录各特征点的SIFT特征向量；[0034] SIFT特征匹配算法表征如下：[0035] （1.1）确定特征点位置坐标与所在尺度。建立图像高斯金字塔，在金字塔尺度空间中的26个邻域中检测极值，若某点（x，y）在金字塔尺度空间本层以及上下两层的26个邻域中是最大或最小值时，定义该点是图像在该尺度L下的一个特征点。[0036] （1.2）特征点方向参数计算。在以特征点（x，y）为中心的邻域窗口内采样，用直方图统计领域像素的梯度方向。将梯度直方图的范围划分为36个方向，每10度代表一个方向。定义直方图峰值所在的方向为该特征点的方向参数。[0037] （1.3）计算SIFT特征向量描述。在特征点中心周围取8*8的窗口，将窗口切成2*2的子窗口。在每个子窗口上计算8个方向的梯度方向直方图，统计每个方向的累加值，作为该子窗口的方向信息。4个子窗口统计完成后最终生成该特征点处的32维度特征向量。[0038] （2）基于SIFT算法对高精度控制影像I2进行特征点提取，得到n个特征点PixI2i（i=1,2…,m），记录各特征点处的SIFT特征向量；[0039] （3）采用欧式距离作为相似性度量准则对两幅影像的特征点进行匹配，得到M个匹配控制点对；[0040] （3.1）对于待检校原始图像I1的任一控制点PixI1i，计算其SIFT特征向量与高精度控制影像I2上提取得到的m个特征向量之间的欧式距离。对于两个向量l1(x1,x2,…,xn)、l2(y1,y2,…,yn)，其欧式距离表征如下：

[0041]

8

CN 103673995 A

说 明 书

5/9页

[0042] （3.2）计算m个欧式距离中的最小值，最小值所在的特征点PixI2j即为待检校原

始图像I1特征点PixI1i的匹配特征点。[0043] （4）根据匹配得到的大量控制点对应在影像上的分布区域和数量，采用格网筛选策略，尽量确保待检校影像各区域都均匀分布一定数量的控制点。如图2所示为VRSS-1遥感卫星待检校影像控制点分布示意图，可以看出采集得到的控制点分布较为均匀。[0044] （5）对每一个匹配特征点对{PixI1i，PixI2j}，记录待检校原始图像I1特征点PixI1i处像点平面坐标（sample,line），以及高精度控制影像I2特征点PixI2j处三维位置坐标（lat,lon,height），作为采集到的高精度控制点信息。[0045] 2.获取相机外检校相机安装补偿矩阵，以及轨道、姿态、行时等辅助数据，建立各控制点处像点到物方点的严格成像模型。[0046] （1）从卫星下传的原始数据中解析待检校影像成像时间范围内的轨道、姿态、行时等数据；[0047] （2）对任一控制点Pointsi{sample,line,lat,lon,height}，根据控制点的图像坐标（sample,line），获取其对应的摄影时间scanTime；

[0048] 控制点对应的摄影时间可以从卫星下传的辅助数据中直接解析，第line行辅助数据对应的成像时间即为该控制点对应的摄影时间。

[0049] （3）利用拉格朗日插值算法插值计算摄影时间scanTime的卫星轨道位置

（PX,PY,PZ,VX,VY,VZ）；卫星按照一定频次下传轨道数据，因此，摄影时间scanTime对应的卫星轨道位置需要利用摄影时刻前后几组的轨道数据进行插值计算。本发明采用拉格朗日插值算法，利用摄影时间前后三组轨道数据计算摄影时间的卫星轨道位置。[0050] （3.1）从第一组轨道数据开始，判断该组轨道数据的生成时间、下一组轨道数据的生成时间与scanTime之间的关系，若scanTime大于第i组轨道数据的生成时间，同时小于第i+1组轨道数据的生成时间，则记录i为距离摄影时间最近的轨道数据序号。[0051] （3.2）利用第i-1，i,i+1组轨道数据，基于拉格朗日算法计算摄影时刻的轨道位置以及速度。拉格朗日插值算法表述如下：[0052] 对于已知y=f(x)的函数表(xi,f(xi))（i=0,1,…,n），对于在[xo,xn]范围内任一x，有：

[0053]

[0054] （4）利用拉格朗日算法插值计算摄影时间scanTime相机相对于轨道坐标系的三

轴姿态角（Roll,Pitch,Yaw）；[0055] （5）建立控制点处的严格成像模型，针对任一控制点，利用外检校安装补偿矩阵、轨道、姿态、行时等数据构建遥感影像的严格几何成像模型Model[0056] 线阵推扫相机的严格成像几何模型如下所述。

[0057]

9

CN 103673995 A[0058]

说 明 书

6/9页

其中：

[0059] [XS,YS,ZS]为该时刻下卫星在协议地心坐标系中的位置，即控制点处轨道位置（PX,PY,PZ）；

[0060] [XG,YG,ZG]为该像元对应的地面目标点在协议地心坐标系中的坐标。[0061] psiX、psiY分别为图像对应的相机像元主光轴单位矢量与卫星本体坐标系X轴、Y轴的夹角

[0062] u—比例因子。

[0063] M为卫星本体坐标系相对于相机的安装补偿矩阵，从外检校参数解算结果文件获取；

[0064] M0为卫星本体坐标系相对于相机初始安装矩阵，从配置文件获取；[0065] M1为该时刻下卫星至轨道坐标系旋转矩阵，由星上测量的姿态角构成。[0066] M2为该时刻下轨道至J2000.0坐标系旋转矩阵，由升交点赤经、轨道倾角、幅角等构成。

[0067] M3为该时刻下J2000.0至WGS84坐标系旋转矩阵，需进行岁差改正、章动改正、格林尼治恒星时改正和极移改正。

[0068] 3.建立相机内方位元素几何检校数学模型，计算各控制点处理论指向矢量。[0069] （1）建立相机内畸变检校模型，在严格成像几何模型中引入已明确的外检校误差矩阵M，如下所示

[0070]

psiX=f(s)[0072] psiY=g(s)[0073] （2）对任一组控制点，计算控制点处的理论指向矢量V理论：[0074] （2.1）利用控制点地面点位置坐标（lat,lon,height）计算地面点在协议地心坐标系中的三轴位置（XG,YG,ZG），计算公式如下：

[0071]

[0075]

[0076] （2.2）利用控制点处的三轴姿态角（Roll,Pitch，Yaw）计算由相机坐标系至轨道坐

标系的旋转矩阵M1[0077] （2.3）利用控制点处的卫星轨道位置（PX,PY,PZ,VX,VY,VZ）计算卫星在J2000.0坐标系下的轨道六根数，由此计算轨道坐标系至J2000.0坐标系的旋转矩阵M2[0078] （2.4）利用控制点处的成像时间，计算J2000.0坐标系至WGS84坐标系的旋转矩阵M3

[0079] （2.5）计算控制点处的理论指向矢量V理论，计算公式为：

[0080]

10

CN 103673995 A

说 明 书

7/9页

[0081] （3）建立控制点理论指向矢量与对应像点列号之间的函数关系。[0082] （3.1）根据指向矢量计算探元在沿轨和垂轨方向的指向角

理论指向矢量对应在沿轨和垂轨方向的指向角：[0084] psiX=atan2(X,Z)[0085] psiY=-atan2(Y,Z)[0086] （3.2）建立基于探元在沿轨和垂轨方向指向角与图像列号的之间函数关系，以三阶多项式为例：

32

[0087] f(s)=A3×S+A2×S+A1×S+A0

[0083]

g(s)=B3×S3+B2×S2+B1×S+B0

[0089] （4）根据每个控制点在不同列号处对应的探元指向角，建立误差方程；针对N个控制点，可得到方程如下：

[0088]

[0090]

建立误差方程:[0092] v=Ax-b[0093] （5）最小二乘法解算误差方程参数，根据模型计算出来的残差剔除粗差，迭代至模型残差在阈值范围内，得到精确内检校参数。[0094] （5.1）针对N个控制点建立的误差方程M1，最小二乘解算方程参数，对于误差方程Ax=b，其最小二乘求解为：

T-1T

[0095] x=(AA+E)(Ab+x)

[0091] [0096]

此时即为该方程的最小二乘解。

[0097] （5.2）将方程M1的最小二乘解作为初步内畸变参数标定结果，代入各控制点处的

严格成像几何模型，计算各控制点处的检校残差，剔除残差大的控制点。

11

CN 103673995 A

说 明 书

8/9页

[0098] （5.3）利用剔除过大残差后的M个控制点，建立误差方程M2，最小二乘迭代解算方

程参数。解算出的最小二乘解作为相机畸变参数标定结果。将标定结果应用于遥感图像几何校正算法中，对严格成像几何模型参数进行修正，可以大幅提高该相机遥感图像产品的几何定位精度和图像内部定位精度。[0099] 如表1-表3所示，表1为本发明解算出的内方位参数标定结果，表2为未使用本发明前VRSS-1卫星几何定位精度测量结果，以及其误差均值和均方根误差，表3为利用本发明进行在轨内方位元素检校后VRSS-1卫星的几何定位精度测量结果，以及其误差均值和均方根误差；从表2可以看出在采用本发明前卫星的实测无控定位精度平均值为34.268米，误差的均方根为10.378。从表3可以看出，采用本发明对相机进行内方位参数标定后实测无控定位精度平均值为30.68米，误差均方根为3.223米。误差的均方根放映了误差偏离平均值的程序，误差均方根值越小说明其内部畸变越小，从表2与表3的计算结果可以看出，采用本发明后，图像的内部畸变得到有效消除，相机的内部定位精度得到了很大提升。[0100] 表1

[0101]

PAN1相机内畸变检校参数A0A1A2A3B0B1B2B3

[0102] [0103]

检校内畸变模型系数结果0.44570751415878540.00000181869345813-0.000000000188596990673790.0000000000000037739095017-1.32292097484831930.000213726116063760.0000000009234494885509-0.00000000000000476011523982

表2

12

CN 103673995 A

说 明 书

9/9页

[0104] [0105]

表3

[0106] 本发明未详细阐述的部分属于本领域公知技术。

13

CN 103673995 A

说 明 书 附 图

1/2页

图1

14

CN 103673995 A

说 明 书 附 图

2/2页

图2

15