CN103148870B - 基于高精度配准控制点的卫星ccd阵列影像几何检校方法 - Google Patents

Abstract

本发明提出了一种基于高精度配准控制点的卫星CCD阵列影像几何检校方法，所述方法利用满足亚像素级配准精度的配准方法按影像列等间隔地提取大量控制点，基于构建的卫星影像定标严密几何模型，求解偏置矩阵；进一步基于定标几何模型求解CCD阵列中的部分探元指向角，并采用多项式模型对解求的探元指向角进行拟合，进而获取CCD阵列中的所有探元指向角。本发明在严密几何模型中引入偏置矩阵，可以很好的补偿由于载荷安装偏差、姿轨偏差引起的几何定位误差；建立探元指向角模型，直接标定探元指向角，避免了像点偏移的模型选择难题，具有较强的适用性，实用价值高。

Description

基于高精度配准控制点的卫星CCD阵列影像几何检校方法

技术领域

本发明属于摄影成像技术领域，特别涉及一种基于高精度配准控制点的卫星CCD阵列影像几何检校方法。

背景技术

目前，国内研制的国产卫星在图像视觉效果和辐射质量上都得到提高，但是几何质量仍和国外同类卫星存在较大的差距。例如，国产卫星的无控制点定位精度都在百米以上，而国外卫星多数在几十米量级；在有足够控制点的条件下，国外卫星影像正射纠正精度约为1个像元，而国产卫星影像正射纠正的精度为3-4个像元。

卫星CCD阵列影像的严密几何成像模型可以用下式简略表示(每颗卫星均会略有差异，因此该式仅为简略表达式)，该式建立了影像上像点与地面点的一一对应关系。

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{\mathrm{WGS}84}={\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]}_{\mathrm{WGS}84}+m{R}_{\mathrm{body}2\mathrm{WGS}84}{R}_{\mathrm{camera}2\mathrm{body}}{\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]}_{\mathrm{camera}}$

式中：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{\mathrm{WGS}84}$ 为像点对应的地物点在WGS84坐标系下的地面坐标；

${\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]}_{\mathrm{WGS}84}$ 为成像时刻卫星在WGS84坐标系下的坐标，其是星上下传的WGS84坐标系下的GPS数据；

R_body2WGS84由星上下传的定姿设备测量数据获得，描述卫星成像姿态；

R_camera2body为载荷在卫星本体坐标系下的安装矩阵，描述载荷与卫星本体的相对关系，其值会在卫星发射前进行精密测量；

${\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]}_{\mathrm{camera}}$ 为影像点在相机坐标系下的坐标，其中x₀、y₀描述相机主点，f描述相机主距。

m为成像比例；

相机内参数误差(镜头畸变、CCD排列变形等)是影响卫星影像正射纠正精度(内部精度)的主要因素。目前，针对内畸变造成的像点偏移研究非常多，并且建立了各种误差模型，如像主点偏移、主距误差、CCD阵列缩放旋转等线性模型，偏心畸变、径向畸变等非线性模型；对于某颗特定卫星传感器，人们无法预先掌握其相机内畸变情况，难以正确选择合适的模型表述其像点偏移；另外，对于相机内部可能存在的非线性模型，也很难正确选择模型阶数，容易因选择阶数过高而造成参数过度化，引起解的不稳定性。

因此，亟需一种新的卫星CCD阵列影像几何检校方法，以避开选择畸变模型这一难题。

发明内容

本发明目的在于解决现有技术不足，提出了一种基于高精度配准控制点的卫星CCD阵列影像几何检校方法，采用偏置矩阵补偿载荷安装误差、姿态及轨道系统误差，采用解求CCD阵列探元的本体坐标系指向角的方法补偿相机畸变误差，避免了像点偏移模型的选择难题。

本发明另外的优点、目的和特性，一部分将在下面的说明书中得到阐明，而另一部分对于本领域的普通技术人员通过对下面的说明的考察将是明显的或从本发明的实施中学到。通过在文字的说明书和权利要求书及附图中特别地指出的结构可实现和获得本发明目的和优点。

本发明提供了一种基于高精度配准控制点的卫星CCD阵列影像几何检校方法，其特征在于，所述方法包括以下步骤：

步骤1.1，针对待检校的卫星CCD阵列影像与高精度数字正射影像(DOM)分辨率的差异，利用所述待检校的卫星CCD阵列影像的内方位元素、辅助姿态、轨道、影像行扫描时间文件和高精度数字高程模型(DEM)影像，构建严密几何成像模型，对所述高精度DOM影像进行重采样，生成一幅与所述待检校的卫星CCD阵列影像相同大小的模拟影像。所述步骤1.1具体包括以下步骤：

步骤1.1.1，利用卫星发射前测量的相机内方位元素，待检校卫星CCD阵列影像的辅助姿态、轨道、影像行扫描时间文件，结合高精度DEM影像，构建严密几何成像模型如下：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{\mathrm{WGS}84}{\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]}_{\mathrm{WGS}84}+m{R}_{\mathrm{body}2\mathrm{WGS}84}{R}_U{R}_{\mathrm{camera}2\mathrm{body}}{\left[\begin{array}{c}x-x_0-\mathrm{\Δx}\\ y-y_0-\mathrm{\Δy}\\ -f\end{array}\right]}_{\mathrm{camera}}$

式中：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{\mathrm{WGS}84}$ 为像点对应的物点在WGS84坐标系下的地面坐标；

${\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]}_{\mathrm{WGS}84}$ 为成像时刻卫星在WGS84坐标系下的坐标，其是星上下传的WGS84坐标系下的GPS数据，从星上下传的轨道文件中获得；

R_body2WGS84由星上下传的姿态文件获得，描述卫星成像姿态；

R_camera2body为载荷在卫星本体坐标系下的安装矩阵，描述载荷与卫星本体的相对关系，其值会在卫星发射前进行精密测量；

m为成像比例；

R_U为偏置矩阵；

${\left[\begin{array}{c}x-x_0-\mathrm{\Δx}\\ y-y_0-\mathrm{\Δy}\\ -f\end{array}\right]}_{\mathrm{camera}}$ 为像点在相机坐标系下的坐标，其由主点(x₀，y₀)、主距f及卫星发射前测量的相机畸变值(Δx，Δy)共同决定；

步骤1.1.2，对模拟影像上的任一个像素(x，y)，其中x为影像上像素的列坐标，y为影像上像素的行坐标，通过所述步骤1.1.1构建的严密几何模型建立对应关系(X，Y，Z)_WGS84＝Tran(x，y)，其中Tran(x，y)表示在步骤1.1.1中构建的几何模型，(X，Y，Z)_WGS84为WGS84椭球下的空间直角坐标，按地图投影转换公式进行坐标转换 ${(X,Y,Z)}_{\mathrm{WGS}84}\⇒(\mathrm{lat},\mathrm{lon},h)\⇒{(\mathrm{east},\mathrm{north})}_{\mathrm{refDom}},$ 将(X，Y，Z)_WGS84转到高精度DOM的投影坐标系下，得到(east，north)_refDom；

步骤1.1.3，将(east，north)_refDom在高精度DOM上进行定位。具体的，获取高精度DOM左上角点坐标(hStart，vStart)和DOM分辨率(hRes，vRes)，则(east，north)_refDom对应的DOM像素点为((east-hStart)/hRes，(north-vStart)/vRes)，对所述DOM像素点进行灰度重采样，得到模拟影像像素(x，y)的灰度；

步骤1.1.4，重复步骤1.1.1-1.1.3，直至生成一幅与检校影像相同大小的模拟影像。

步骤1.2，采用相位相关的高精度配准算法，对所述待检校的卫星CCD阵列影像与模拟影像进行列等间隔配准，获取控制点；所述步骤1.2具体包括以下步骤：

步骤1.2.1，根据所述待检校的卫星CCD阵列影像的CCD阵列大小，确定配准列间隔；具体的，假定所述待检校的卫星CCD阵列影像的大小为w·h，w为影像宽，h为影像高，取配准列间隔为n，n为一正整数，则所述配准将在所述待检校的卫星CCD阵列影像的n、2n，3n......in列上进行，它们对应的探元号记为x_n、x_2n，x_3n......x_i·n(x_i·n≤w)；

步骤1.2.2，针对所有所述待检校的卫星CCD阵列影像的配准列x_i·n(x_i·n≤w)，对每个配准列下的h个影像点逐一与模拟影像配准，获取配准点对(x_i·n，y，x′，y′，其中(x_i·n，y)为模拟影像上的第x_i·n列，第y行上的像素，(x′，y′)为对应的所述待检校的卫星CCD阵列影像；

步骤1.2.3，所述待检校的卫星CCD阵列影像与模拟影像之间优选的用仿射模型 $\begin{array}{c}{x}_{i\·n}=a_0+a_1{x}^{\′}+a_2{y}^{\′}\\ y=b_0+b_1{x}^{\′}+b_2{y}^{\′}\end{array}$ 建立对应关系；依据该对应关系，利用步骤1.2.2获取的所有配准点通过最小二乘解求得系数(a₀ a₁ a₂ b₀ b₁ b₂)、模型单位权中误差、以及各配准点平差后的残差，将残差大于所述模型单位权中误差的各配准点剔除；重复进行该操作，直至单位权中误差小于一个预定的阈值；

步骤1.2.4，利用步骤1.1.1构建的严密几何成像模型，将各配准点对(x_i·n，y，x′，y′)对应的模拟影像像素(x_i·n，y)计算到地面坐标：(X，Y，Z)_WGS84＝Tran(x_i·n，y)，从而得到所述待检校的卫星CCD阵列影像的控制点(x′，y′，X，Y，Z)；

步骤1.3，利用通过步骤1.2.4得到的所有控制点来解求偏置矩阵，以利用求得的偏置矩阵来补偿载荷安装误差、姿态和轨道系统误差；求解所述偏置矩阵R_U的步骤具体如下：

所述偏置矩阵R_U记为：

采用指向角模型表述内方位元素：

令：

${\left[\begin{array}{c}\stackrel{\‾}{x}\\ \stackrel{\‾}{y}\\ \stackrel{\‾}{z}\end{array}\right]}_{\mathrm{body}}=R_{\mathrm{camera}2\mathrm{body}}{\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]}_{\mathrm{camera}}$

上式中：

R_camera2body为载荷在卫星本体坐标系下的安装矩阵，描述载荷与卫星本体的相对关系，其值会在卫星发射前进行精密测量；

${\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]}_{\mathrm{camera}}$ 为影像点在相机坐标系下的坐标，其中x₀、y₀描述相机主点，f描述相机主距。

由上式得出坐标后，

易知：

${\mathrm{\ψ}}_x=\arctan \left(\frac{\stackrel{\‾}{x}}{\stackrel{\‾}{z}}\right)$

${\mathrm{\ψ}}_y=\arctan \left(\frac{\stackrel{\‾}{y}}{\stackrel{\‾}{z}}\right)$

则步骤1中严密几何成像模型可变换为：

$\left[\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\right)\\ \tan \left({\mathrm{\ψ}}_y\right)\\ 1\end{array}\right]=\frac{1}{m}{R}_U^{-1}{R}_{\mathrm{body}2\mathrm{WGS}84}^{-1}{\left[\begin{array}{c}X-X_S\\ Y-Y_S\\ Z-Z_S\end{array}\right]}_{\mathrm{WGS}84}$

令 ${\left[\begin{array}{c}\stackrel{\‾}{X}\\ \stackrel{\‾}{Y}\\ \stackrel{\‾}{Z}\end{array}\right]}_{\mathrm{body}}=R_{\mathrm{body}2\mathrm{WGS}84}^{-1}{\left[\begin{array}{c}X-X_S\\ Y-Y_S\\ Z-Z_S\end{array}\right]}_{\mathrm{WGS}84}$ 则上式可化为：

$\frac{a_1\stackrel{\‾}{X}+b_1\stackrel{\‾}{Y}+c_1\stackrel{\‾}{Z}}{a_3\stackrel{\‾}{X}+b_3\stackrel{\‾}{Y}+c_3\stackrel{\‾}{Z}}-\tan \left({\mathrm{\ψ}}_x\right)=0$

$\frac{a_2\stackrel{\‾}{X}+b_2\stackrel{\‾}{Y}+c_2\stackrel{\‾}{Z}}{a_3\stackrel{\‾}{X}+b_3\stackrel{\‾}{Y}+c_3\stackrel{\‾}{Z}}-\tan \left({\mathrm{\ψ}}_y\right)=0$

令

$f\left(x\right)=\frac{a_1\stackrel{\‾}{X}+b_1\stackrel{\‾}{Y}+c_1\stackrel{\‾}{Z}}{a_3\stackrel{\‾}{X}+b_3\stackrel{\‾}{Y}+c_3\stackrel{\‾}{Z}}-\tan \left({\mathrm{\ψ}}_x\right)$

$f\left(y\right)=\frac{a_2\stackrel{\‾}{X}+b_2\stackrel{\‾}{Y}+c_2\stackrel{\‾}{Z}}{a_3\stackrel{\‾}{X}+b_3\stackrel{\‾}{Y}+c_3\stackrel{\‾}{Z}}-\tan \left({\mathrm{\ψ}}_y\right)$

列出误差方程：V＝Bx-l，W

式中， $V=\left[\begin{array}{c}{V}_{f\left(x\right)}\\ V_{f\left(y\right)}\end{array}\right]$ $l=\left(\begin{array}{c}-f_x^0\\ -f_y^0\end{array}\right),$

W为权矩阵；

利用最小二乘原理：x＝(B^TB)^-1B^Tl，由此可计算出偏置矩阵三个角元素从而求解出偏置矩阵。

步骤1.4，将解求得到的偏置矩阵引入严密成像几何模型，利用控制点解求待检校的卫星CCD阵列探元指向角，最终获取CCD阵列中每个探元在卫星本体坐标系下的探元指向角；所述步骤1.4具体包括以下步骤：

步骤1.4.1，针对配准列x_i·n(x_i·n≤w)，利用该列下经过步骤1.2.4配准获得的多个控制点来解求探元x_i·n在本体坐标系下的指向角；具体的，

将步骤1中的严密成像几何模型转换为：

${\left[\begin{array}{c}X-X_S\\ Y-Y_S\\ Z-Z_S\end{array}\right]}_{\mathrm{WGS}84}=m{R}_{\mathrm{body}2\mathrm{WGS}84}{R}_U{\left[\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\right)\\ \tan \left({\mathrm{\ψ}}_y\right)\\ 1\end{array}\right]}_{\mathrm{body}}$

令 $\left[\begin{array}{c}\mathrm{XX}\\ \mathrm{YY}\\ \mathrm{ZZ}\end{array}\right]={\left[\begin{array}{c}X-X_S\\ Y-Y_S\\ Z-Z_S\end{array}\right]}_{\mathrm{WGS}84},$ $R=R_{\mathrm{body}2\mathrm{WGS}84}{R}_U\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]$

$\frac{a_1\tan \left({\mathrm{\ψ}}_x\right)+a_2\tan \left({\mathrm{\ψ}}_y\right)+a_3}{c_1\tan \left({\mathrm{\ψ}}_x\right)+c_2\tan \left({\mathrm{\ψ}}_y\right)+c_3}-\frac{\mathrm{XX}}{\mathrm{ZZ}}=0$

上式可化为： $\frac{b_1\tan \left({\mathrm{\ψ}}_x\right)+b_2\tan \left({\mathrm{\ψ}}_y\right)+b_3}{c_1\tan \left({\mathrm{\ψ}}_x\right)+c_2\tan \left({\mathrm{\ψ}}_y\right)+c_3}-\frac{\mathrm{YY}}{\mathrm{ZZ}}=0$

$f_x=\frac{a_1\tan \left({\mathrm{\ψ}}_x\right)+a_2\tan \left({\mathrm{\ψ}}_y\right)+a_3}{c_1\tan \left({\mathrm{\ψ}}_x\right)+c_2\tan \left({\mathrm{\ψ}}_y\right)+c_3}-\frac{\mathrm{XX}}{\mathrm{ZZ}}$

令 $f_y=\frac{b_1\tan \left({\mathrm{\ψ}}_x\right)+b_2\tan \left({\mathrm{\ψ}}_y\right)+b_3}{c_1\tan \left({\mathrm{\ψ}}_x\right)+c_2\tan \left({\mathrm{\ψ}}_y\right)+c_3}-\frac{\mathrm{YY}}{\mathrm{ZZ}}$

列出误差方程：V＝Bx-l，W

式中， $V=\left[\begin{array}{c}{V}_{f_x}\\ V_{f_y}\end{array}\right]$ $B=\left(\begin{array}{cc}\frac{\∂f_x}{\∂\left(\tan \left({\mathrm{\ψ}}_x\right)\right)}& \frac{\∂f_x}{\∂\left(\tan \left({\mathrm{\ψ}}_y\right)\right)}\\ \frac{\∂f_y}{\∂\left(\tan \left({\mathrm{\ψ}}_x\right)\right)}& \frac{\∂f_y}{\∂\left(\tan \left({\mathrm{\ψ}}_y\right)\right)}\end{array}\right)$ $l=\left(\begin{array}{c}-f_x^0\\ -f_y^0\end{array}\right),$ x＝(tan(ψ_x)tan(ψ_y))，

W为权矩阵。

利用最小二乘原理：x＝(B^TB)^-1B^Tl，由此可计算出x_i·n探元的指向角正切值(tan(ψ_x)tan(ψ_y))，进而得到x_i·n探元的探元指向角(ψ_x ψ_y)。

步骤1.4.2，对各配准列x_i·n(x_i·n≤w)执行步骤1.4.1，根据步骤1.2.3删除误匹配点后保留下来的高精度配准点分布情况，解求得到多个探元的探元指向角：(tan(ψ_x)tan(ψ_y))_n、(tan(ψ_x)tan(ψ_y))_2n……，n为配准列间隔，对解求获取的所有探元指向角进行多项式拟合，一般取不超过5次多项式。

tan(ψ_x)＝a₀+a₁x+a₂x²+……+a_ixⁱ

tan(ψ_y)＝b₀+b₁x+b₂x²+……+b_ixⁱ，i≤5

其中x表示探元编号，对应于待检校图像列号；基于最小二乘解求得到参数a₀、a₁……a_i；b₀、b₁……b_i，i≤5

步骤1.4.3，取0≤x＜w-1，逐一将x代入步骤1.4.2中的多项式，得到所有探元的探元指向角。

优选的，在步骤1.1.1中，所述相机内方位元素例如可为主点、主距和畸变。

本发明在严密几何模型中引入偏置矩阵，可以很好的补偿由于载荷安装偏差、姿轨偏差引起的几何定位误差；建立探元指向角模型，直接标定探元指向角，避免了像点偏移的模型选择难题，具有较强的适用性，实用价值高。

附图说明

图1是根据本发明实施例的、卫星本体坐标系下的探元指向角的示意图。

图2是根据本发明实施例的、一种基于高精度配准控制点的卫星CCD阵列影像几何检校方法的流程图。

具体实施方式

本发明针对卫星CCD阵列影像，采用亚像素级的高精度配准算法提取大量控制点，通过在严密几何成像模型中引入偏置矩阵补偿姿轨外方位元素误差、解求CCD阵列探元指向角补偿内方位元素误差的方法，消除影响几何定位精度的外方位元素、内方位元素系统误差，提高卫星几何定位精度。

本发明中载荷泛指星上设备，如星敏、相机等；“探元”与“CCD探元”含义相同，不做区分。

本发明在严密几何成像模型中引入偏置矩阵补偿姿轨外方位元素误差，下面对此作出详细解释。卫星在轨运行时，由于外界物理环境的变化，会导致载荷安装值(如R_camera2body)偏离发射前的地面标定值，造成几何定位的系统偏差；另外，针对卫星成像误差的研究表明，在卫星成像的较短时段内，姿态、轨道误差引起的几何定位偏差呈现系统性，且两者具有较强的相关性。为此，在星历姿态上引入偏置矩阵R_U，综合补偿载荷安装、姿态轨道误差，便可大幅度地提高卫星影像的无控制点定位精度。

下式描述在严密几何模型中引入偏置矩阵：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{\mathrm{WGS}84}={\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]}_{\mathrm{WGS}84}+m{R}_{\mathrm{body}2\mathrm{WGS}84}{R}_U{R}_{\mathrm{camera}2\mathrm{body}}{\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]}_{\mathrm{camera}}$

如前面所描述的，在上式中：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{\mathrm{WGS}84}$ 为像点对应的物点在WGS84坐标系下的地面坐标；

${\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]}_{\mathrm{WGS}84}$ 为成像时刻卫星在WGS84坐标系下的坐标，其是星上下传的WGS84坐标系下的GPS数据；

R_body2WGS84由星上下传的定姿设备测量数据获得，描述卫星成像姿态；

R_camera2body为载荷在卫星本体坐标系下的安装矩阵，描述载荷与卫星本体的相对关系，其值会在卫星发射前进行精密测量；

${\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]}_{\mathrm{camera}}$ 为影像点在相机坐标系下的坐标，其中x₀、y₀描述相机主点，f描述相机主距。

m为成像比例；

R_U为偏置矩阵。

下面结合图1来描述卫星本体坐标系下的探元指向角概念：

如图1所示，其中，表示卫星本体坐标系，矢量表示CCD阵列上某探元在坐标系中的指向，ψ_x、ψ_y即为该探元的沿轨指向角和垂轨指向角。

令：

${\left[\begin{array}{c}\stackrel{\‾}{x}\\ \stackrel{\‾}{y}\\ \stackrel{\‾}{z}\end{array}\right]}_{\mathrm{body}}=R_{\mathrm{camera}2\mathrm{body}}{\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]}_{\mathrm{camera}}$

上式中：

R_camera2body为载荷在卫星本体坐标系下的安装矩阵，描述载荷与卫星本体的相对关系，其值会在卫星发射前进行精密测量；

${\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]}_{\mathrm{camera}}$ 为影像点在相机坐标系下的坐标，其中x₀、y₀描述相机主点，f描述相机主距。

由上式得出坐标后，

易知：

${\mathrm{\ψ}}_x=\arctan \left(\frac{\stackrel{\‾}{x}}{\stackrel{\‾}{z}}\right)$

${\mathrm{\ψ}}_y=\arctan \left(\frac{\stackrel{\‾}{y}}{\stackrel{\‾}{z}}\right)$

为了尽可能正确地获取所有探元的本体指向，需先解求一定数量的探元指向角，进而通过多项式拟合出所有探元指向角。这要求有足够多的地面控制点。通过常规的布设靶标方式获取控制点成本较高，不宜采用。因此，本发明选择已有的高精度数字正射影像(Digital Orthophoto Map，DOM)、数字高程模型(DigitalElevation Model，DEM)，通过亚像素级的高精度配准算法将卫星影像与参考DOM影像进行配准，获得满足需求的控制点，用以解求一定数量的探元指向角。

下面结合附图2，对本发明做详细地说明。本发明所提出的基于高精度配准控制点的卫星CCD阵列影像几何检校方法具体包括以下步骤：

步骤1.1，针对待检校的卫星CCD阵列影像与高精度数字正射影像(DOM)分辨率的差异，利用所述待检校的卫星CCD阵列影像的内方位元素、辅助姿态、轨道、影像行扫描时间文件和高精度数字高程模型(DEM)影像，构建严密几何成像模型，对所述高精度DOM影像进行重采样，生成一幅与所述待检校的卫星CCD阵列影像相同大小的模拟影像。所述步骤1.1具体包括以下步骤：

步骤1.1.1，利用卫星发射前测量的相机内方位元素(例如，主点、主距和畸变)，待检校卫星CCD阵列影像的辅助姿态、轨道、影像行扫描时间文件，结合高精度DEM影像，构建严密几何成像模型如下：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{\mathrm{WGS}84}{\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]}_{\mathrm{WGS}84}+m{R}_{\mathrm{body}2\mathrm{WGS}84}{R}_U{R}_{\mathrm{camera}2\mathrm{body}}{\left[\begin{array}{c}x-x_0-\mathrm{\Δx}\\ y-y_0-\mathrm{\Δy}\\ -f\end{array}\right]}_{\mathrm{camera}}$

式中：

${\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]}_{\mathrm{WGS}84}$ 为像点对应的物点在WGS84坐标系下的地面坐标；

${\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right]}_{\mathrm{WGS}84}$ 为成像时刻卫星在WGS84坐标系下的坐标，其是星上下传的WGS84坐标系下的GPS数据，从星上下传的轨道文件中获得；

R_body2WGS84由星上下传的姿态文件获得，描述卫星成像姿态；

R_camera2body为载荷在卫星本体坐标系下的安装矩阵，描述载荷与卫星本体的相对关系，其值会在卫星发射前进行精密测量；

m为成像比例；

R_U为偏置矩阵。

${\left[\begin{array}{c}x-x_0-\mathrm{\Δx}\\ y-y_0-\mathrm{\Δy}\\ -f\end{array}\right]}_{\mathrm{camera}}$ 为像点在相机坐标系下的坐标，其由主点(x₀，y₀)、主距f及卫星发射前测量的相机畸变值(Δx，Δy)共同决定。

步骤1.1.2，对模拟影像上的任一个像素(x，y)，其中x为影像上像素的列坐标，y为影像上像素的行坐标，通过所述步骤1.1.1构建的严密几何模型建立对应关系(X，Y，Z)_WGS84＝Tran(x，y)，其中Tran(x，y)表示在步骤1.1.1中构建的几何模型，(X，Y，Z)_WGS84为WGS84椭球下的空间直角坐标，按地图投影转换公式进行坐标转换 ${(X,Y,Z)}_{\mathrm{WGS}84}\⇒(\mathrm{lat},\mathrm{lon},h)\⇒{(\mathrm{east},\mathrm{north})}_{\mathrm{refDom}},$ 将(X，Y，Z)_WGS84转到高精度DOM的投影坐标系下，得到(east，north)_refDom。

利用地图投影转换公式进行的坐标转换对于本领域技术人员来说是熟知的，在此不做详细描述。

步骤1.1.3，将(east，north)_refDom在高精度DOM上进行定位。具体的，获取高精度DOM左上角点坐标(hStart，vStart)和DOM分辨率(hRes，vRes)，则(east，north)_refDom对应的DOM像素点为((east-hStart)/hRes，(north-vStart)/vRes)，对所述DOM像素点进行灰度重采样，得到模拟影像像素(x，y)的灰度；

步骤1.1.4，重复步骤1.1.1-1.1.3，直至生成一幅与所述待检校的卫星CCD阵列影像相同大小的模拟影像。

步骤1.2，采用相位相关的高精度配准算法，对所述待检校的卫星CCD阵列影像与模拟影像进行列等间隔配准，获取控制点。所述步骤1.2具体包括以下步骤：

步骤1.2.1，根据所述待检校的卫星CCD阵列影像的CCD阵列大小，确定配准列间隔。假定所述待检校的卫星CCD阵列影像的大小为w.h，w为影像宽，h为影像高，取配准列间隔为n，n为一正整数，则所述配准将在所述待检校的卫星CCD阵列影像的n、2n，3n……in列上进行，它们对应的探元号(等同于所述待检校的卫星CCD阵列影像的列坐标)记为x_n、x_2n，x_3n……x_i·n(x_i·n≤w)；

步骤1.2.2，针对所有所述待检校的卫星CCD阵列影像的配准列x_i·n(x_i·n≤w)，对每个配准列下的h个影像点逐一与模拟影像配准，获取配准点对(x_i·n，y，x′，y′)，其中(x_i·n，y)为模拟影像上的第x_i·n列，第y行上的像素，(x′，y′)为对应的所述待检校的卫星CCD阵列影像；

步骤1.2.3，所述待检校的卫星CCD阵列影像与模拟影像之间用仿射模型 $\begin{array}{c}{x}_{i\·n}=a_0+a_1{x}^{\′}+a_2{y}^{\′}\\ y=b_0+b_1{x}^{\′}+b_2{y}^{\′}\end{array}$ 建立对应关系；依据该对应关系，利用步骤1.2.2获取的所有配准点通过最小二乘解求得系数(a₀ a₁ a₂ b₀ b₁ b₂)、模型单位权中误差、以及各配准点平差后的残差，将残差大于所述模型单位权中误差的各配准点剔除；重复进行该操作，直至单位权中误差小于一个预定的阈值；

所述阈值根据检校精度要求而定，例如检校精度要求误差小于1个像素，那么阈值即为1个像素，也就是说，单位权中误差小于1个像素。

步骤1.2.4，利用步骤1.1.1构建的严密几何成像模型，将各配准点对(x_i·n，y，x′，y′)对应的模拟影像像素(x_i·n，y)计算到地面坐标：(X，Y，Z)_WGS84＝Trar(x_i·n，y)，从而得到所述待检校的卫星CCD阵列影像的控制点(x′，y′，X，Y，Z)。

步骤1.3，利用通过步骤1.2.4得到的所有控制点来解求偏置矩阵，以利用求得的偏置矩阵来补偿载荷安装误差、姿态和轨道系统误差。

偏置矩阵R_U的解求方法按摄影测量中后方交会方法进行。求解所述偏置矩阵R_U的步骤具体如下：

所述偏置矩阵R_U记为：

采用指向角模型表述内方位元素：

令：

${\left[\begin{array}{c}\stackrel{\‾}{x}\\ \stackrel{\‾}{y}\\ \stackrel{\‾}{z}\end{array}\right]}_{\mathrm{body}}=R_{\mathrm{camera}2\mathrm{body}}{\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]}_{\mathrm{camera}}$

上式中：

R_camera2body为载荷在卫星本体坐标系下的安装矩阵，描述载荷与卫星本体的相对关系，其值会在卫星发射前进行精密测量；

${\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]}_{\mathrm{camera}}$ 为影像点在相机坐标系下的坐标，其中x₀、y₀描述相机主点，f描述相机主距。

由上式得出坐标后，

易知：

${\mathrm{\ψ}}_x=\arctan \left(\frac{\stackrel{\‾}{x}}{\stackrel{\‾}{z}}\right)$

${\mathrm{\ψ}}_y=\arctan \left(\frac{\stackrel{\‾}{y}}{\stackrel{\‾}{z}}\right)$

则步骤1中严密几何成像模型可变换为：

$\left[\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\right)\\ \tan \left({\mathrm{\ψ}}_y\right)\\ 1\end{array}\right]=\frac{1}{m}{R}_U^{-1}{R}_{\mathrm{body}2\mathrm{WGS}84}^{-1}{\left[\begin{array}{c}X-X_S\\ Y-Y_S\\ Z-Z_S\end{array}\right]}_{\mathrm{WGS}84}$

令 ${\left[\begin{array}{c}\stackrel{\‾}{X}\\ \stackrel{\‾}{Y}\\ \stackrel{\‾}{Z}\end{array}\right]}_{\mathrm{body}}=R_{\mathrm{body}2\mathrm{WGS}84}^{-1}{\left[\begin{array}{c}X-X_S\\ Y-Y_S\\ Z-Z_S\end{array}\right]}_{\mathrm{WGS}84}$ 则上式可化为：

$\frac{a_1\stackrel{\‾}{X}+b_1\stackrel{\‾}{Y}+c_1\stackrel{\‾}{Z}}{a_3\stackrel{\‾}{X}+b_3\stackrel{\‾}{Y}+c_3\stackrel{\‾}{Z}}-\tan \left({\mathrm{\ψ}}_x\right)=0$

$\frac{a_2\stackrel{\‾}{X}+b_2\stackrel{\‾}{Y}+c_2\stackrel{\‾}{Z}}{a_3\stackrel{\‾}{X}+b_3\stackrel{\‾}{Y}+c_3\stackrel{\‾}{Z}}-\tan \left({\mathrm{\ψ}}_y\right)=0$

令

$f\left(x\right)=\frac{a_1\stackrel{\‾}{X}+b_1\stackrel{\‾}{Y}+c_1\stackrel{\‾}{Z}}{a_3\stackrel{\‾}{X}+b_3\stackrel{\‾}{Y}+c_3\stackrel{\‾}{Z}}-\tan \left({\mathrm{\ψ}}_x\right)$

$f\left(y\right)=\frac{a_2\stackrel{\‾}{X}+b_2\stackrel{\‾}{Y}+c_2\stackrel{\‾}{Z}}{a_3\stackrel{\‾}{X}+b_3\stackrel{\‾}{Y}+c_3\stackrel{\‾}{Z}}-\tan \left({\mathrm{\ψ}}_y\right)$

列出误差方程：V＝Bx-l，W

式中， $V=\left[\begin{array}{c}{V}_{f\left(x\right)}\\ V_{f\left(y\right)}\end{array}\right]$ $l=\left(\begin{array}{c}-f_x^0\\ -f_y^0\end{array}\right),$

W为权矩阵。

利用最小二乘原理：x＝(B^TB)^-1B^Tl，由此可计算出偏置矩阵三个角元素从而求解出偏置矩阵。

若实施时由于参数的相关性，使得求解过程中方程出现病态，一般的解决方式是岭参数估计法。

步骤1.4，将解求得到的偏置矩阵引入严密成像几何模型，利用控制点解求待检校的卫星CCD阵列探元指向角，最终获取CCD阵列中每个探元在卫星本体坐标系下的探元指向角；所述步骤1.4具体包括以下步骤：

步骤1.4.1，针对配准列x_i·n(x_i·n≤w)，利用该列下经过步骤1.2.4配准获得的多个控制点来解求探元x_i·n在本体坐标系下的指向角；具体的，

将步骤1中的严密成像几何模型转换为：

${\left[\begin{array}{c}X-X_S\\ Y-Y_S\\ Z-Z_S\end{array}\right]}_{\mathrm{WGS}84}=m{R}_{\mathrm{body}2\mathrm{WGS}84}{R}_U{\left[\begin{array}{c}\tan \left({\mathrm{\ψ}}_x\right)\\ \tan \left({\mathrm{\ψ}}_y\right)\\ 1\end{array}\right]}_{\mathrm{body}}$

令 $\left[\begin{array}{c}\mathrm{XX}\\ \mathrm{YY}\\ \mathrm{ZZ}\end{array}\right]={\left[\begin{array}{c}X-X_S\\ Y-Y_S\\ Z-Z_S\end{array}\right]}_{\mathrm{WGS}84},$ $R=R_{\mathrm{body}2\mathrm{WGS}84}{R}_U\left[\begin{array}{ccc}{a}_1& a_2& a_3\\ b_1& b_2& b_3\\ c_1& c_2& c_3\end{array}\right]$

$\frac{a_1\tan \left({\mathrm{\ψ}}_x\right)+a_2\tan \left({\mathrm{\ψ}}_y\right)+a_3}{c_1\tan \left({\mathrm{\ψ}}_x\right)+c_2\tan \left({\mathrm{\ψ}}_y\right)+c_3}-\frac{\mathrm{XX}}{\mathrm{ZZ}}=0$

上式可化为： $\frac{b_1\tan \left({\mathrm{\ψ}}_x\right)+b_2\tan \left({\mathrm{\ψ}}_y\right)+b_3}{c_1\tan \left({\mathrm{\ψ}}_x\right)+c_2\tan \left({\mathrm{\ψ}}_y\right)+c_3}-\frac{\mathrm{YY}}{\mathrm{ZZ}}=0$

$f_x=\frac{a_1\tan \left({\mathrm{\ψ}}_x\right)+a_2\tan \left({\mathrm{\ψ}}_y\right)+a_3}{c_1\tan \left({\mathrm{\ψ}}_x\right)+c_2\tan \left({\mathrm{\ψ}}_y\right)+c_3}-\frac{\mathrm{XX}}{\mathrm{ZZ}}$

令 $f_y=\frac{b_1\tan \left({\mathrm{\ψ}}_x\right)+b_2\tan \left({\mathrm{\ψ}}_y\right)+b_3}{c_1\tan \left({\mathrm{\ψ}}_x\right)+c_2\tan \left({\mathrm{\ψ}}_y\right)+c_3}-\frac{\mathrm{YY}}{\mathrm{ZZ}}$

列出误差方程：V＝Bx-l，W

式中， $V=\left[\begin{array}{c}{V}_{f_x}\\ V_{f_y}\end{array}\right]$ $B=\left(\begin{array}{cc}\frac{\∂f_x}{\∂\left(\tan \left({\mathrm{\ψ}}_x\right)\right)}& \frac{\∂f_x}{\∂\left(\tan \left({\mathrm{\ψ}}_y\right)\right)}\\ \frac{\∂f_y}{\∂\left(\tan \left({\mathrm{\ψ}}_x\right)\right)}& \frac{\∂f_y}{\∂\left(\tan \left({\mathrm{\ψ}}_y\right)\right)}\end{array}\right)$ $l=\left(\begin{array}{c}-f_x^0\\ -f_y^0\end{array}\right),$ x＝(tan(ψ_x)tan(ψ_y))，

W为权矩阵。

利用最小二乘原理：x＝(B^TB)^-1B^Tl，由此可计算出x_i·n探元的指向角正切值(tan(ψ_x)tan(ψ_y))，进而得到x_i·n探元的探元指向角(ψ_x ψ_y)。

步骤1.4.2，对各配准列x_i·n(x_i·n≤w)执行步骤1.4.1，根据步骤1.2.3删除误匹配点后保留下来的高精度配准点分布情况，解求得到多个探元的探元指向角：(tan(ψ_x)tan(ψ_y))_n、(tan(ψ_x)tan(ψ_y))_2n……，n为配准列间隔，对解求获取的所有探元指向角进行多项式拟合，一般取不超过5次多项式。

tan(ψ_x)＝a₀+a₁x+a₂x²+……+a_ixⁱ

tan(ψ_y)＝b₀+b₁x+b₂x²+……+b_ixⁱ，i≤5

其中x表示探元编号，对应于待检校图像列号；基于最小二乘解求得到参数a₀、a₁……a_i；b₀、b₁……b_i，i≤5

步骤1.4.3，取0≤x＜w-1，逐一将x代入步骤1.4.2中的多项式，得到所有探元的探元指向角。

本发明在严密几何模型中引入偏置矩阵，可以很好的补偿由于载荷安装偏差、姿轨偏差引起的几何定位误差；建立探元指向角模型，直接标定探元指向角，避免了像点偏移的模型选择难题，具有较强的适用性，实用价值高。

以上内容仅为本发明的较佳实施例，对于本领域的普通技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，本说明书内容不应理解为对本发明的限制。

Claims (2)

1.一种基于高精度配准控制点的卫星CCD阵列影像几何检校方法，其特征在于，所述方法包括以下步骤：

步骤1.1，针对待检校的卫星CCD阵列影像与高精度数字正射影像DOM分辨率的差异，利用所述待检校的卫星CCD阵列影像的内方位元素、辅助姿态、轨道、影像行扫描时间文件和高精度数字高程模型DEM影像，构建严密几何成像模型，对所述高精度DOM影像进行重采样，生成一幅与所述待检校的卫星CCD阵列影像相同大小的模拟影像；所述步骤1.1具体包括以下步骤：

步骤1.1.1，利用卫星发射前测量的相机内方位元素，待检校卫星CCD阵列影像的辅助姿态、轨道、影像行扫描时间文件，结合高精度DEM影像，构建严密几何成像模型如下：

式中：

为像点对应的物点在WGS84坐标系下的地面坐标；

为成像时刻卫星在WGS84坐标系下的坐标，其是星上下传的WGS84坐标系下的GPS数据，从星上下传的轨道文件中获得；

R_body2WGS84由星上下传的姿态文件获得，描述卫星成像姿态；

R_camera2body为载荷在卫星本体坐标系下的安装矩阵，描述载荷与卫星本体的相对关系，其值会在卫星发射前进行精密测量；

m为成像比例；

R_U为偏置矩阵；

为像点在相机坐标系下的坐标，其由主点(x₀，y₀)、主距f及卫星发射前测量的相机畸变值(Δx，Δy)共同决定；

步骤1.1.2，对模拟影像上的任一个像素(x，y)，其中x为影像上像素的列坐标，y为影像上像素的行坐标，通过所述步骤1.1.1构建的严密几何模型建立对应关系(X，Y，Z)_WGS84＝Tran(x，y)，其中Tran(x，y)表示在步骤1.1.1中构建的几何模型，(X，Y，Z)_WGS84为WGS84椭球下的空间直角坐标，按地图投影转换公式进行坐标转换将(X，Y，Z)_WGS84转到高精度DOM的投影坐标系下，得到(east，north)_refDom；

步骤1.1.3，将(east，north)_refDom在高精度DOM上进行定位；具体的，获取高精度DOM左上角点坐标(hStart，vStart)和DOM分辨率(hRes，vRes)，则(east，north)_refDom对应的DOM像素点为((east-hStart)/hRes，(north-vStart)/vRes)，对所述DOM像素点进行灰度重采样，得到模拟影像像素(x，y)的灰度；

步骤1.1.4，重复步骤1.1.1-1.1.3，直至生成一幅与待检校的卫星CCD阵列影像相同大小的模拟影像；

步骤1.2，采用相位相关的高精度配准算法，对所述待检校的卫星CCD阵列影像与模拟影像进行列等间隔配准，获取控制点；所述步骤1.2具体包括以下步骤：

步骤1.2.1，根据所述待检校的卫星CCD阵列影像的CCD阵列大小，确定配准列间隔；具体的，假定所述待检校的卫星CCD阵列影像的大小为w·h，w为影像宽，h为影像高，取配准列间隔为n，n为一正整数，则所述配准将在所述待检校的卫星CCD阵列影像的n、2n，3n......in列上进行，它们对应的探元号记为x_n、x_2n，x_3n......x_i·n(x_i·n≤w)；

步骤1.2.2，针对所有所述待检校的卫星CCD阵列影像的配准列x_i·n(x_i·n≤w)，对每个配准列下的h个影像点逐一与模拟影像配准，获取配准点对(x_i·n，y，x′，y′，其中(x_i·n，y)为模拟影像上的第x_i·n列，第y行上的像素，(x′，y′)为对应的所述待检校的卫星CCD阵列影像；

步骤1.2.3，所述待检校的卫星CCD阵列影像与模拟影像之间用仿射模型x_i·n＝a₀+a₁x′+a₂y′

y＝b₀+b₁x′+b₂y′建立对应关系；依据该对应关系，利用步骤1.2.2获取的所有 配准点通过最小二乘解求得系数(a₀ a₁ a₂ b₀ b₁ b₂)、模型单位权中误差、以及各配准点平差后的残差，将残差大于所述模型单位权中误差的各配准点剔除；重复进行该操作，直至单位权中误差小于一个预定的阈值；

步骤1.2.4，利用步骤1.1.1构建的严密几何成像模型，将各配准点对(x_i·n，y，x′，y′)对应的模拟影像像素(x_i·n，y)计算到地面坐标：(X，Y，Z)_WGS84＝Tran(xi_·n，y)，从而得到所述待检校的卫星CCD阵列影像的控制点(x′，y′，X，Y，Z)；

步骤1.3，利用通过步骤1.2.4得到的所有控制点来解求偏置矩阵，以利用求得的偏置矩阵来补偿载荷安装误差、姿态和轨道系统误差；求解所述偏置矩阵R_U的步骤具体如下：

所述偏置矩阵R_U记为：

采用指向角模型表述内方位元素：

令：

上式中：

R_camera2body为载荷在卫星本体坐标系下的安装矩阵，描述载荷与卫星本体的相对关系，其值会在卫星发射前进行精密测量；

为影像点在相机坐标系下的坐标，其中x₀、y₀描述相机主点，f描述相机主距；

由上式得出坐标后，

易知：

则步骤1中严密几何成像模型可变换为：

令则上式可化为：

令

列出误差方程：V＝Bx-l，W

式中，

W为权矩阵；

利用最小二乘原理：x＝(B^TB)^-1B^Tl，由此可计算出偏置矩阵三个角元素 从而求解出偏置矩阵；

步骤1.4，将解求得到的偏置矩阵引入严密成像几何模型，利用控制点解求待检校的卫星CCD阵列探元指向角，最终获取CCD阵列中每个探元在卫星本体坐标系下的探元指向角；所述步骤1.4具体包括以下步骤：

步骤1.4.1，针对配准列x_i·n(x_i·n≤w)，利用该列下经过步骤1.2.4配准获得的多 个控制点来解求探元x_i·n在本体坐标系下的指向角；具体的，

将步骤1中的严密成像几何模型转换为：

令

上式可化为：

令

列出误差方程：V＝Bx-l，W

式中，x＝(tan(ψ_x)tan(ψ_y))，W为权矩阵；

利用最小二乘原理：x＝(B^TB)^-1B^Tl，由此可计算出x_i·n探元的指向角正切值(tan(ψ _x)tan(ψ_y))，进而得到x_i·n探元的探元指向角(ψ_x ψ_y)；

步骤1.4.2，对各配准列x_i·n(x_i·n≤w)执行步骤1.4.1，根据步骤1.2.3删除误匹配点后保留下来的高精度配准点分布情况，解求得到多个探元的探元指向角：(tan(ψ_x)tan(ψ_y))_n、(tan(ψ_x)tan(ψ_y))_2n……，n为配准列间隔，对解求获取的所有探元指向角进行多项式拟合，取不超过5次多项式；

tan(ψ_x)＝a₀+a₁x+a₂x²+……+a_ixⁱ

tan(ψ_y)＝b₀+b₁x+b₂x²+……+b_ixⁱ，i≤5

其中x表示探元编号，对应于待检校图像列号；基于最小二乘解求得到参数a₀、a₁……a_i；b₀、b₁……b_i，i≤5

步骤1.4.3，取0≤x＜w-1，逐一将x代入步骤1.4.2中的多项式，得到所有探元的探元指向角。

2.根据权利要求1所述的基于高精度配准控制点的卫星CCD阵列影像几何检校方法，其特征在于，在步骤1.1.1中，所述相机内方位元素例如可为主点、主距和畸变。