CN101839713B - 一种基于带时间因子偏置矩阵的卫星影像系统误差改正方法 - Google Patents

Abstract

基于带时间因子偏置矩阵的卫星影像系统误差改正方法，将系统误差改正问题归结为求解相机绕卫星本体的两个偏置角，并针对图像几何定位误差在某段在轨测试期间的线性变化规律，用带时间因子的一次多项式表示偏置角，并基于选取的样本图像和自动提取的若干控制点平差解算该一次多项式系数，从而得到该段时间内任意卫星影像的偏置矩阵，代入共线条件方程式即实现了对几何定位系统误差的有效补偿和改正。在轨测试结果表明，该方法简单易行，能有效消除随时间变化的系统误差，保持系统几何纠正产品精度的长期稳定性。

Description

一种基于带时间因子偏置矩阵的卫星影像系统误差改正方法

技术领域

本发明属于测绘科学与技术领域，涉及一种基于带时间因子偏置矩阵的卫星影像系统误差改正方法。

背景技术

近年来，我国相继发射了多颗高分辨率光学遥感卫星，为科学试验、国土资源普查、农作物估产和防灾减灾等领域提供了大量的数据，在我国国民经济的发展中发挥着重要作用。全色CCD相机是高分辨率卫星的重要成像载荷，一般与星体刚性固联，采用线阵推扫方式成像。在卫星入轨后，对该卫星的地面应用系统进行相关性能测试，从而为有效载荷参数的调整和成像数据预处理算法的改进提供依据，并对影像产品质量及应用潜力进行综合评价。

卫星遥感影像的系统几何纠正，是指由卫星轨道星历、姿态以及相机传感器参数等建立直接对地定位模型，然后基于该模型将影像纳入到一定的地图投影坐标系的过程。导致卫星影像对地定位误差的因素有很多，包括CCD位置偏差、相机焦距测量误差、相机安装角偏差、卫星位置测量误差、卫星速度测量误差、卫星时间校准误差、卫星姿态角测量误差以及各项测量数据的量化误差等；卫星上天之后，其地面测量值也可能发生变化。这些因素都可能导致遥感影像对地定位模型存在系统性误差，影响系统级几何纠正产品的定位精度。因此，有效改正或补偿卫星影像的系统误差，是遥感卫星地面预处理系统的关键技术。

在建立了系统误差补偿模型之后，需要利用地面控制点通过平差解算得到系统误差补偿参数。但由于各观测值系统误差之间存在强相关性，即复共线性，如相机侧摆角与相机安装角之间、星时误差与平台俯仰角之间等，联合求解会导致法方程严重病态，甚至奇异，使得求解精度不高甚至无法求解。目前，系统误差在轨测试分析方法以基于偏置矩阵的安置角常差检校和基于广义偏移矩阵的综合误差改正方法为主。理论上讲，要精确计算偏置矩阵，需要使用尽可能精确的姿态和星历数据。资源卫星应用中心最早利用地面控制点和对应的影像坐标计算了资源一号卫星相机的偏置矩阵，有效补偿了由于CCD安装误差等引起的卫星遥感影像的系统误差，提高了对地观测精度；袁修孝和张过等进而分析了三个正交方向上偏置角对目标定位精度的影响，提出了对三个偏置角分步解算的新方法。事实上，当卫星遥感影像存在一个明显的方向一致的系统误差时，其原因可能是由于传感器的安装误差，或者是姿态数据或者星历数据的系统误差，或者是地面处理模型的系统误差等。要从总的误差中分离、提取出各个因素误差，目前还有相当的难度。基于广义偏移矩阵的系统误差改正方法将相机安装误差、姿态星历数据误差、地面处理模型误差等所综合导致的成像偏差，归结为一个旋转矩阵来校正，国产资源卫星的几何在轨测试证明该方法能有效提高图像预处理几何定位精度。

近来，在国产卫星的几何在轨测试过程中，发现了一种新的卫星影像直接对地定位误差，它在一定的时间范围内随时间缓慢线性变化，具有系统误差的规律性变化的特点。由于常规基于偏置矩阵或广义偏置矩阵的系统误差改正模型是与时间无关的，故不能解决目前在轨测试遇到的新问题。因此，在理论分析的基础上，提出了一种基于带时间因子偏置矩阵的卫星影像系统误差改正方法。尽管导致国产卫星直接定位误差变化规律的因素很复杂，但是试验表明，通过建立这样一种带时间因子的偏置矩阵误差改正模型，能有效地改善遥感卫星影像直接对地定位的精度，保持系统几何纠正产品精度的长期稳定性。

发明内容

本发明所要解决的问题是：实现一种基于带时间因子偏置矩阵的卫星影像系统误差改正方法。

本发明提供的技术方案是：基于带时间因子偏置矩阵的卫星影像系统误差改正方法，将系统误差改正问题归结为求解相机绕卫星本体的两个偏置角，并针对图像几何定位误差在某段在轨测试期间的线性变化规律，用带时间因子的一次多项式表示偏置角，并基于选取的样本图像和自动提取的若干控制点平差解算该一次多项式系数，从而得到该段时间内任意卫星影像的偏置矩阵，代入共线条件方程式即实现了对几何定位系统误差的有效补偿和改正。该技术方案具体包括以下步骤：

一、样本图像选取。

从连续M个月在轨获取的卫星影像中，抽取若干景作为样本图像，这里，M取不大于6的正整数。样本图像获取的日期各不相同，且均匀覆盖该时间段；令样本图像个数为S，S≥1；计算每一景样本图像的相对成像时刻t_i，其中i＝1，2，3，...，S，t_i为从该段时间起始日期起算的天数(如果t₁为该段时间起始日期的第1天，则t₁=1)；

二、控制点提取。

以控制点影像或地形图数据为参考，提取控制点并得到其在样本图像上的像点坐标；每景样本图像上的控制点数为N_i(i＝1，2，3，...，S)，其中N_i＞0，控制点总数

三、基于样本图像和控制点计算偏置角随时间变化的一次多项式系数。

1.将每一景样本图像摄影光线的对地定位偏差分解为沿轨道和垂直于轨道两个方向，令

和

代表摄影光线定位误差分解到这两个方向后对应的偏置角度，根据星载线阵推扫影像的严格成像几何模型，对第j(j＝1，2，...，N)个控制点列共线条件方程式，如公式(1)：

$\left[\begin{array}{c}{X}_j\\ Y_j\\ Z_j\end{array}\right]=\left[\begin{array}{c}{\mathrm{Xs}}_j\\ {\mathrm{Ys}}_j\\ {\mathrm{Zs}}_j\end{array}\right]+m_j{R}_{C2\mathrm{Tj}}{R}_{\mathrm{GFj}}{R}_{\mathrm{FBj}}{R}_{\mathrm{vi}}{R}_{\mathrm{BS}}\left[\begin{array}{c}{x}_j\\ y_j\\ -f\end{array}\right]---\left(1\right)$

式中，

(X_j，Y_j，Z_j)是控制点的地心直角坐标；

(Xs_j，Ys_j，Zs_j)是控制点成像时刻相机投影中心的地心直角坐标；

(x_j，y_j)是控制点对应像点的相机焦平面坐标；

f是相机的焦距；

m_j是比例尺系数；

R_BS是相机坐标系到卫星本体坐标系的旋转矩阵；

R_Vi是偏置矩阵；

R_FBj是本体坐标系到轨道坐标系的旋转矩阵，即姿态矩阵；

R_GFj是轨道坐标系到地心惯性系的旋转矩阵；

R_C2Tj是从地心惯性系到地固系的旋转矩阵；

i是控制点所在样本图像的序号。

用带时间因子的一次多项式表示构造R_vi的两个偏置角

和

如公式(2)所示；

式中，t_i表示第i景样本图像的相对成像时刻；

为一次多项式系数。

2.对系数

赋初值0，即R_vi初始为单位矩阵；

3.将

的值代入公式(3)；

$\left[\begin{array}{c}{\stackrel{\‾}{X}}_j\\ {\stackrel{\‾}{Y}}_j\\ {\stackrel{\‾}{Z}}_j\end{array}\right]=\left[\begin{array}{c}{\mathrm{Xs}}_j\\ {\mathrm{Ys}}_j\\ {\mathrm{Zs}}_j\end{array}\right]+m_j{R}_{C2\mathrm{Tj}}{R}_{\mathrm{GFj}}{R}_{\mathrm{FBj}}{R}_v{R}_{\mathrm{BS}}\left[\begin{array}{c}{x}_j\\ y_j\\ -f\end{array}\right]---\left(3\right)$

式中，

是由

的迭代值计算得到的控制点地心直角坐标。

4.基于公式(3)，对第j(j＝1，2，...，N)个控制点列误差方程式(4)；

A_jQ＝L_j    (4)

$L_j=\left[\begin{array}{c}{X}_j-{\stackrel{\‾}{X}}_j\\ Y_j-{\stackrel{\‾}{Y}}_j\\ Z_j-{\stackrel{\‾}{Z}}_j\end{array}\right]$

式中，A_j是误差方程式的系数矩阵；

Q是由待求系数迭代累加值构成的向量；

L_j是误差向量，反映了控制点坐标计算值与真实值之间的差异。

5.基于公式(5)计算法方程系数矩阵；

$A^{T}A=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_j^T{A}_j---\left(5\right)$

$A^{T}L=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_j^T{L}_j$

6.基于公式(6)，平差计算待求系数的迭代累加值向量Q；

Q＝(A^TA)^-1(A^TL)    (6)

7.计算

的累加值，返回执行步骤3进行迭代计算，当

均小于阈值10^-6时，迭代停止。

四、系统误差改正：

当系统误差在某段时间内呈近似线性变化时，基于求解的一次多项式系数

对于该段时间内的某一景卫星影像，根据其相对成像时刻构造偏置矩阵，

代入共线条件方程式之后，能有效地改进影像的对地定位精度。

本发明方法简单易行，能有效消除或补偿遥感卫星影像随时间变化的系统误差，保证几何定位精度的长期稳定性。

在构成传统偏置矩阵的偏置角中添加带时间因子的一次项，代入共线方程严格成像几何模型后，从而建立了卫星影像随时间变化的系统误差改正模型。在系统误差呈近似线性变化的每一段时间内，基于地面控制点列误差方程式，最小二乘平差解算该偏置矩阵的一次多项式系数。在轨测试的大量实验证明该方法简单易行，能有效消除随时间变化的系统误差，保证产品几何定位精度的长期稳定性。

附图说明

图1是相机偏置角与其直接对地定位误差的关系示意图。

图2是系统误差改正前相机沿CCD线阵方向几何定位误差随时间的变化趋势图。

图3是系统误差改正前相机沿轨道方向几何定位误差随时间的变化趋势图。

图4是系统误差改正后相机沿CCD线阵方向几何定位误差随时间的变化趋势图。

图5是系统误差改正后相机沿轨道方向几何定位误差随时间的变化趋势图。

具体实施方式

下面结合附图对本发明做进一步详细描述。

为了保证卫星影像的直接对地定位精度在规定范围之内，必须建立相应的系统误差改正模型。严格意义上讲，偏置矩阵是由于传感器与卫星平台之间的安装偏差所导致的传感器坐标系与卫星平台坐标系之间的正交旋转矩阵，用来校正传感器与卫星平台坐标系之间不重合而导致的成像偏差，偏置角即构成该旋转矩阵的三个转角。常规的系统误差改正方法将系统误差模型简化为用沿轨道方向和垂直轨道方向的两个角度表示，并假设相机安装角偏差是主要误差来源。基于广义偏移矩阵的系统误差改正方法将相机安装误差、星历姿态测量误差等所导致的综合成像偏差，归结为一个旋转矩阵来校正。如图1所示，轨道坐标系O-XYZ的原点O位于卫星质心，Y坐标轴指向卫星飞行方向，Z坐标轴指向地球质心，X轴根据右手法则确定，令影像上某像点的投影中心为S，这里S与O重合，将其摄影光线SP的对地定位偏差分解为沿轨道和垂直于轨道(沿CCD线阵)两个方向，

和

代表摄影光线定位误差分解到这两个方向后对应的偏置角度，图中Po为真实摄影光线与地面的交点。不难发现，对地定位的偏差值越大，该方向上的偏置角越大。

实施例：取连续9个月时间内获取的共33景样本影像数据(33景样本影像成像时间各不相同，且尽量均匀分布在该连续9个月时间内)，影像覆盖了包括北京、河北、浙江、湖北、山东在内的多个省市地区。在每景影像上量取了若干检查点，检查点在每景影像上大致均匀分布。检查各景影像上检查点定位误差的大小和方向：不考虑相机的安置误差，直接基于原始轨道、姿态和相机的初始安装角等辅助数据构建共线条件方程式，计算各景内检查点沿CCD方向和沿轨道方向的定位误差，绘制各景影像在两个方向上的定位中误差图。如图2和图3所示，坐标横轴表示景编号(景编号按照时间先后顺序排列)，纵轴分别表示各景内所有检查点沿两个方向的几何定位中误差(单位：像素)。从图中可以发现，按照时间顺序排列后，各景影像的定位误差还表现出趋势性和转折点，即在一段连续时间内，定位精度有随时间线性变化的趋势，但经过某个转折点后，仍然有随时间线性变化的趋势，但变化的方向反转了。大概在第24景(位于第5个月)处存在一个明显的转折点。

第一阶段：样本图像选取。

选取前面连续5个月在轨获取的卫星影像作为样本图像，令样本图像个数为S，S≥1；计算每一景样本图像的相对成像时刻t_i，其中i＝1，2，3，...，S，t_i为从该段时间起始日期起算的天数(如果t₁为该段时间起始日期的第1天，则t₁＝1)；

第二阶段：控制点提取。

以控制点影像或地形图数据为参考，在每景样本图像上量取若干控制点；每景样本图像上的控制点数为N_i(i＝2，3，...，S)，其中N_i＞0，控制点总数

第三阶段：基于样本图像和控制点计算偏置角随时间变化的一次多项式系数。

1.将每一景样本图像摄影光线的对地定位偏差分解为沿轨道和垂直于轨道两个方向，令

和代表摄影光线定位误差分解到这两个方向后对应的偏置角度，根据星载线阵推扫影像的严格成像几何模型，对第j(j＝1，2，...，N)个控制点列共线条件方程式，如公式(1)：

$\left[\begin{array}{c}{X}_j\\ Y_j\\ Z_j\end{array}\right]=\left[\begin{array}{c}{\mathrm{Xs}}_j\\ {\mathrm{Ys}}_j\\ {\mathrm{Zs}}_j\end{array}\right]+m_j{R}_{C2\mathrm{Tj}}{R}_{\mathrm{GFj}}{R}_{\mathrm{FBj}}{R}_{\mathrm{vi}}{R}_{\mathrm{BS}}\left[\begin{array}{c}{x}_j\\ y_j\\ -f\end{array}\right]---\left(1\right)$

式中，

(X_j，Y_i，Z_j)是控制点的地心直角坐标；

(X_Sj，Y_Sj，Z_Sj)是控制点成像时刻相机投影中心的地心直角坐标；

(x_j，y_j)是控制点对应像点的相机焦平面坐标；

f是相机的焦距；

m_j是比例尺系数；

R_BS是相机坐标系到卫星本体坐标系的旋转矩阵；

R_Vi是偏置矩阵；

R_FBj是本体坐标系到轨道坐标系的旋转矩阵，即姿态矩阵；

R_GFj是轨道坐标系到地心惯性系的旋转矩阵；

R_C2Tj是从地心惯性系到地固系的旋转矩阵；

i是控制点所在样本图像的编号。

用带时间因子的一次多项式表示构造R_vi的两个偏置角

和

如公式(2)所示；

式中，t_i表示第i景样本图像的相对成像时刻；

为一次多项式系数。

2.对系数赋初值0，即R_vi初始为单位矩阵；

3.将的值代入公式(3)；

$\left[\begin{array}{c}{\stackrel{\‾}{X}}_j\\ {\stackrel{\‾}{Y}}_j\\ {\stackrel{\‾}{Z}}_j\end{array}\right]=\left[\begin{array}{c}{\mathrm{Xs}}_j\\ {\mathrm{Ys}}_j\\ {\mathrm{Zs}}_j\end{array}\right]+m_j{R}_{C2\mathrm{Tj}}{R}_{\mathrm{GFj}}{R}_{\mathrm{FBj}}{R}_v{R}_{\mathrm{BS}}\left[\begin{array}{c}{x}_j\\ y_j\\ -f\end{array}\right]---\left(3\right)$

式中，

是由

的迭代值计算得到的控制点地心直角坐标。

4.基于公式(3)，对第j(j＝1，2，...，N)个控制点列误差方程式(4)；

A_jQ＝L_j    (4)

$L_j=\left[\begin{array}{c}{X}_j-{\stackrel{\‾}{X}}_j\\ Y_j-{\stackrel{\‾}{Y}}_j\\ Z_j-{\stackrel{\‾}{Z}}_j\end{array}\right]$

式中，A_j是误差方程式的系数矩阵；

Q是由待求系数迭代累加值构成的向量；

L_j是误差向量，反映了控制点坐标计算值与真实值之间的差异。

5.基于公式(5)计算法方程系数矩阵；

$A^{T}A=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_j^T{A}_j---\left(5\right)$

$A^{T}L=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A}_j^T{L}_j$

6.基于公式(6)，平差计算待求系数

的迭代累加值向量Q；

Q＝(A^TA)^-1(A^TL)    (6)

7.计算

的累加值，返回执行步骤3进行迭代计算，当

均小于阈值10^-6时，迭代停止；

第四阶段：系统误差改正。

基于求解的一次多项式系数

对于前5个月在轨获取的任意一景卫星影像，根据其相对成像时刻构造偏置矩阵，代入共线条件方程式之后，能有效地改进影像的对地定位精度。

同样的，选取后面连续4个月在轨获取的卫星影像作为样本图像，求解其对应的一次多项式系数，从而实现对该段时间的卫星对地定位系统误差的改进。

采用该方法进行系统误差补偿后，样本图像上检查点的直接对地定位精度不再随时间变化，如图4和图5所示。试验表明，系统误差补偿之前，直接对地定位精度在百米到1公里之间，随时间呈分段线性变化，引入带时间的偏置矩阵之后，系统几何纠正产品的定位精度能够保持在150米(1σ)左右，已没有明显的时间相关性。

Claims (2)

1.一种基于带时间因子偏置矩阵的卫星影像系统误差改正方法，包括如下步骤：

一、样本图像选取：

从连续M个月在轨获取的卫星影像中，抽取若干景作为样本图像；M取不大于6的正整数；令样本图像个数为S，S≥1；计算每一景样本图像的相对成像时刻t_i，其中i＝1，2，3，...，S，t_i为从连续M个月的起始日期起算的天数；

二、控制点提取：

以控制点影像或地形图数据为参考，提取控制点并得到其在样本图像上的像点坐标；每景样本图像上的控制点数为N_i，其中，i＝1，2，3，...，S，N_i＞0，控制点总数

三、基于样本图像和控制点计算偏置角随时间变化的一次多项式系数：

1)将每一景样本图像摄影光线的对地定位偏差分解为沿轨道和垂直于轨道两个方向，令

和代表摄影光线定位误差分解到这两个方向后对应的偏置角度，根据星载线阵推扫影像的严格成像几何模型，对第j个控制点列共线条件方程式(1)，j＝1，2，...，N；

$\left[\begin{array}{c}{X}_j\\ Y_j\\ Z_j\end{array}\right]=\left[\begin{array}{c}X{s}_j\\ Y{s}_j\\ Z{s}_j\end{array}\right]+m_j{R}_{C2\mathrm{Tj}}{R}_{\mathrm{GFj}}{R}_{\mathrm{FBj}}{R}_{\mathrm{vi}}{R}_{\mathrm{BS}}\left[\begin{array}{c}{x}_j\\ y_j\\ -f\end{array}\right]---\left(1\right)$

式中，

(X_j，Y_j，Z_j)是控制点的地心直角坐标；

(X_Sj，Y_Sj，Z_Sj)是控制点成像时刻相机投影中心的地心直角坐标；

(x_j，y_j)是控制点对应像点的相机焦平面坐标；

f是相机的焦距；

m_j是比例尺系数；

R_BS是相机坐标系到卫星本体坐标系的旋转矩阵；

R_Vi是偏置矩阵；

R_FBj是本体坐标系到轨道坐标系的旋转矩阵，即姿态矩阵；

R_GFj是轨道坐标系到地心惯性系的旋转矩阵；

R_C2Tj是从地心惯性系到地固系的旋转矩阵；

i是控制点所在样本图像的编号；

用带时间因子的一次多项式(2)表示构造R_vi的两个偏置角

和

式中，t_i表示第i景样本图像的相对成像时刻；

为一次多项式系数；

2)对系数

赋初值0，即R_vi初始为单位矩阵；

3)将的值代入公式(3)；

$\left[\begin{array}{c}{\stackrel{\‾}{X}}_j\\ {\stackrel{\‾}{Y}}_j\\ {\stackrel{\‾}{Z}}_j\end{array}\right]=\left[\begin{array}{c}X{s}_j\\ Y{s}_j\\ Z{s}_j\end{array}\right]+m_j{R}_{C2\mathrm{Tj}}{R}_{\mathrm{GFj}}{R}_{\mathrm{FBj}}{R}_v{R}_{\mathrm{BS}}\left[\begin{array}{c}{x}_j\\ y_j\\ -f\end{array}\right]---\left(3\right)$

式中

是由

的迭代值计算得到的控制点地心直角坐标；

4)基于公式(3)，对第j个控制点列误差方程式(4)；

A_jQ＝L_j         (4)

$L_j=\left[\begin{array}{c}{X}_j-{\stackrel{\‾}{X}}_j\\ Y_j-{\stackrel{\‾}{Y}}_j\\ Z_j-{\stackrel{\‾}{Z}}_j\end{array}\right]$

式中，A_j是误差方程式的系数矩阵；

Q是由待求系数迭代累加值构成的向量；

L_j是误差向量，反映了控制点坐标计算值与真实值之间的差异；

5)基于公式(5)计算法方程系数矩阵；

$A^{T}A=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A_j}^T{A}_j$

                  (5)

$A^{T}L=\underset{j=1}{\overset{N}{\mathrm{\Σ}}}{A_j}^T{L}_j$

6)基于公式(6)，平差计算待求系数

的迭代累加值向量Q；

Q＝(A^TA)^-1(A^TL)         (6)

7)计算

的累加值，返回执行步骤3)进行迭代计算，当

均小于阈值10^-6时，迭代停止；

四、系统误差改正：

当系统误差在某段时间内呈近似线性变化时，基于求解的一次多项式系数

对于该段时间内的某一景卫星影像，根据其相对成像时刻构造偏置矩阵，代入共线条件方程式进行计算。

2.根据权利要求1所述的方法，其特征在于：步骤一中样本图像获取的日期各不相同，且均匀分布在M个月内。

Images