CN103778610A - 一种星载线阵传感器垂轨摆扫影像的几何预处理方法 - Google Patents

Abstract

一种星载线阵传感器垂轨摆扫影像的几何预处理方法，包括根据原始单帧图像的成像几何，建立原始单帧图像的共线方程模型；对每个原始单帧图像进行几何纠正处理，包括构建物方局部坐标与地心直角坐标的相互换算关系，构建切平面图像坐标与物方局部坐标之间的相互换算关系，对每个原始单帧图像划分虚拟立体格网，解算原始单帧图像对应的有理多项式模型系数，建立每个原始单帧图像坐标与切平面图像坐标之间的正反计算关系，基于有理多项式模型对各原始单帧图像进行几何纠正，得到物方切平面坐标系下的帧图像；对所得物方切平面坐标系下的所有帧图像进行基于坐标的拼接，得到拼接后的图像；解算拼接后图像对应的有理多项式模型系数。

Description

一种星载线阵传感器垂轨摆扫影像的几何预处理方法

技术领域

本发明属于航天与航空摄影测量领域，涉及星载线阵传感器在垂轨摆扫成像方式下的高精度几何预处理方法。

背景技术

线阵CCD（电荷耦合元件图像传感器）是目前获取高分辨率光学卫星影像的主要传感器。在高分辨率光学成像过程中，由于镜头焦距长导致CCD观测视场窄。采用垂轨摆扫成像方式代替沿轨推扫成像方式可以有效增大观测视场角。然而，在垂轨摆扫成像方式下，由于反光镜摆扫运动导致成像关系比较复杂、所得图像的几何畸变很大，这在高精度几何应用中亟需解决。

常规的线阵摆扫影像通过格网点坐标内插来实现系统几何纠正，如果用户希望利用地面控制点对系统几何纠正产品进行更高精度的几何处理就只能采用常规的多项式模型，这样，产品的几何精度（特别是内部几何精度）对地面控制点的分布、数量有较大的依赖性，产品的几何精度难以得到保障。

发明内容

本发明所要解决的问题是：针对星载线阵传感器摆扫成像方式，提供一种高精度的几何预处理方法，以解决图像内部几何畸变大的问题；而且，通过给出图像的有理多项式模型，能够使用户可以利用地面控制点对图像进行后续更高精度的几何处理。

本发明的技术方案为一种星载线阵传感器垂轨摆扫影像的几何预处理方法，进行以下步骤，

步骤1，根据原始单帧图像的成像几何，建立原始单帧图像的共线方程模型；

步骤2，对每个原始单帧图像进行几何纠正处理，包括以下子步骤，

步骤2.1，建立物方局部坐标系和物方切平面坐标系，构建物方局部坐标与地心直角坐标的相互换算关系；

所述物方局部坐标系O-XYZ定义为，经过原始图像中心点的光线与物方平均高程面H₀相交的交点为物方局部坐标系的原点O；XY为物方切平面，与参考椭球面相切，Z轴垂直于XY切平面向上为正；X轴沿飞行方向，Y轴垂直于飞行方向，XYZ组成右手系；物方局部坐标系下的坐标记为物方局部坐标；

物方切平面坐标系O-XY定义为，物方切平面坐标系的原点O、X轴、Y轴都与物方局部坐标系的原点O、X轴、Y轴重合，X轴、Y轴的刻度单位为米；物方切平面坐标系下的坐标记为物方切平面坐标；

步骤2.2，建立切平面图像坐标系，并构建切平面图像坐标与物方局部坐标之间的相互换算关系；

所述切平面图像坐标系O’-X’Y’定义为，位于物方切平面上，切平面图像坐标系的X’轴、Y’轴与物方切平面坐标系的X轴、Y轴平行，X’轴、Y’轴的刻度单位为像素，原点O’为物方切平面坐标系下的图像范围的最小外接矩形的左上角点；切平面图像坐标系下的坐标记为切平面图像坐标；

步骤2.3，对每个原始单帧图像划分虚拟立体格网，计算虚拟立体格网点的地心直角坐标，并将虚拟立体格网点作为控制点，解算原始单帧图像对应的有理多项式模型系数；

步骤2.4，根据步骤2.3所得每个原始单帧图像的有理多项式模型和步骤2.1、步骤2.2，建立每个原始单帧图像坐标与切平面图像坐标之间的正反计算方式；

步骤2.5，基于步骤2.4根据有理多项式模型建立的正反计算方式，对各原始单帧图像进行几何纠正，得到物方切平面坐标系下的帧图像；

步骤3，对步骤2.5所得物方切平面坐标系下的所有帧图像进行基于坐标的拼接，得到拼接后的图像；

步骤4，解算步骤3所得拼接后图像对应的有理多项式模型系数。

而且，步骤1中，建立原始单帧图像的共线方程模型如下，

设R_t、R_GF、R_FB、R_BS分别为t时刻下从地心惯性坐标系到地固地心直角坐标系的旋转矩阵、从轨道坐标系到地心惯性坐标系的旋转矩阵、从本体坐标系到轨道坐标系的旋转矩阵、从相机坐标系到本体坐标系的旋转矩阵,[X_t Y_t Z_t]^T为t时刻卫星质心在地心直角坐标系下的坐标矢量，[x(l) y(l) 1]^T为通过像点p(s,l)的光线在相机坐标系下的矢量，则其对应物点P在地心直角坐标系下的坐标[X_P Y_P Z_P]^T为：

$\left[\begin{array}{c}{X}_P\\ Y_P\\ Z_P\end{array}\right]=\left[\begin{array}{c}{X}_t\\ Y_t\\ Z_t\end{array}\right]+{\mathrm{mR}}_t{R}_{\mathrm{GF}}{R}_{\mathrm{FB}}{R}_{\mathrm{BS}}\left[\begin{array}{c}x\left(l\right)\\ y\left(l\right)\\ 1\end{array}\right]---\left(1\right)$

其中，(s,l)为像点的原始单帧图像坐标，l为原始单帧图像的行号，s为原始单帧图像的列号。

而且，步骤2.1中，构建物方局部坐标与地心直角坐标的相互换算关系如下，

设物方局部坐标系原点的地心直角坐标为(X₀,Y₀,Z₀)，地理坐标为(L₀,B₀)，B₀为纬度，L₀为经度，任一物点P的地心直角坐标为(X_P,Y_P,Z_P)，物点P的物方局部坐标为(x_P,y_P,z_P)，

从地心直角坐标(X_P,Y_P,Z_P)到物方局部坐标(x_P,y_P,z_P)的换算公式为，

$\left[\begin{array}{c}{x}_P\\ y_P\\ z_P\end{array}\right]=R_{\mathrm{WGS}84}^{\mathrm{LOC}}\left[\begin{array}{c}{X}_P-X_0\\ Y_P-Y_0\\ Z_P-Z_0\end{array}\right]---\left(2\right)$

其中， $R_{\mathrm{WGS}84}^{\mathrm{LOC}}=\left[\begin{array}{ccc}-\sin \left(L_0\right)& \cos \left(L_0\right)& 0\\ -\sin \left(B_0\right)\cos \left(L_0\right)& -\sin \left(B_0\right)\sin \left(L_0\right)& \cos \left(B_0\right)\\ \cos \left(B_0\right)\cos \left(L_0\right)& \cos \left(B_0\right)\sin \left(L_0\right)& \sin \left(B_0\right)\end{array}\right]---\left(3\right)$

从物方局部坐标(x_P,y_P,z_P)到地心直角坐标(X_P,Y_P,Z_P)的换算公式为，

$\left[\begin{array}{c}{X}_P\\ Y_P\\ Z_P\end{array}\right]=\left[\begin{array}{c}{X}_0\\ Y_0\\ Z_0\end{array}\right]+{\left(R_{\mathrm{WGS}84}^{\mathrm{LOC}}\right)}^T\left[\begin{array}{c}{x}_P\\ y_P\\ z_P\end{array}\right]---\left(4\right)$

其中，

为从地心直角坐标系到物方局部坐标系的旋转矩阵。

而且，步骤2.2中，构建切平面图像坐标与物方局部坐标之间的相互换算关系如下，

设切平面图像坐标系原点的物方局部坐标为(x₀,y₀,z₀)，任一物点P的物方局部坐标为(x_P,y_P,z_P)，相应物方切平面图像坐标为(s_P,l_P)，

从物方局部坐标(x_P,y_P,z_P)到切平面图像坐标(s_P,l_P)的换算公式为：

s_P＝(x_p-x₀)/d，l_P＝(y_p-y₀)/d       （5）

从切平面图像坐标(s_P,l_P)到物方局部坐标(x_P,y_P,z_P)的换算公式为：

x_P＝x₀+s_p×d，y_P＝y₀+l_p×d      （6）

其中，d为图像分辨率。

而且，步骤2.3中，对某原始单帧图像划分虚拟立体格网，计算虚拟立体格网点的地心直角坐标实现方式如下，

设星载线阵传感器的探元数为N，每次扫描的子帧数为W，对应地面的高程范围为[H_min,H_max]，则原始单帧图像大小为N行×W列；假设经过某像点p(s,l)的光线与K个高程面H₁,H₂,...,H_K相交，(s,l)为像点的原始单帧图像坐标，0≤l≤N-1，0≤s≤W-1，H_min≤H_k≤H_max，1≤k≤K，得到K个交点P₁,P₂,...,P_K，根据公式（1）求出这些交点对应的地心直角坐标[X_P1 Y_P1 Z_P1]^T,[X_P2 Y_P2 Z_P2]^T,...,[X_PK Y_PK Z_PK]^T，进而求得这些交点对应的地理坐标(L_P1,B_P1),(L_P2,B_P2),...,(L_Pk,B_PK)。

而且，步骤2.4中，根据有理多项式模型，建立每个原始单帧图像坐标与切平面图像坐标之间的正反计算方式如下，

设原始单帧图像对应的有理多项式模型为：

s＝f_x(B,L,H),l＝f_y(B,L,H)       （7）

其中f_x,f_y为有理多项式函数，(B,L,H)为物点的地理坐标，(s,l)为对应像点的原始单帧图像坐标；

基于有理多项式模型的坐标反算，包括首先由切平面图像坐标按照公式（6）求得物方切平面坐标(x_P,y_P)，物方切平面坐标(x_P,y_P)加上z_P=0得到物方局部坐标(x_P,y_P,z_P)；然后，先根据公式（4）将物方局部坐标(x_P,y_P,z_P)转换为地心直角坐标(X_P,Y_P,Z_P)，再计算出地理坐标作为(B,L,H)代入公式（7）求出原始单帧图像坐标；

基于有理多项式模型的坐标正算，包括以下子步骤，

步骤a，设物点的初始高程值H⁽⁰⁾为平均高程，将原始单帧图像坐标和物点的初始高程值代入公式（1）求得物点的地心直角坐标[X_P Y_P Z_P]^T，再根据公式（2）求得物方局部坐标的初始结果(x⁽⁰⁾,y⁽⁰⁾,z⁽⁰⁾)；

步骤b，令迭代次数k=1；

步骤c，令物点的高程值H^(k)＝H^(k-1)-z^(k-1)，再代入公式（1）求得物点的地心直角坐标作为新的[X_P Y_P Z_P]^T，再根据公式（2）求得物方局部坐标的第k次迭代结果(x^(k),y^(k),z^(k))；

步骤d，判断是否|z^(k)|＜T_Z，T_Z为预设的阈值，是则迭代终止，输出当前的物方局部坐标(x^(k),y^(k),z^(k))，否则令k=k+1，返回步骤c继续迭代；

步骤e，根据公式（5）将物方局部坐标换算成切平面图像坐标。

而且，步骤4中，解算步骤3所得拼接后图像对应的有理多项式模型系数实现如下，

首先，对步骤2.3所得各原始单帧图像的虚拟立体格网点取并集后，重新计算这些虚拟立体格网点的图像坐标，包括先根据公式（2）、（3）将各虚拟立体格网点的地理坐标换算成物方局部坐标，再根据所得公式（5）将物方局部坐标换算成切平面图像坐标；

然后，将所有坐标变换后的虚拟立体格网点作为控制点，解算出拼接后图像对应的有理多项式模型系数。

而且，所述地心直角坐标和地理坐标都是定义在以地球质心为原点的参考椭球上的，以地球质心为原点的参考椭球为WGS84椭球。

通过以上步骤就可以实现星载线阵传感器垂轨摆扫影像的几何预处理。本发明针对星载线阵传感器的摆扫成像方式，基于物方切平面坐标系对各摆扫帧图像（即原始单帧图像）进行几何预纠正，不仅消除各摆扫帧图像的内部几何畸变，实现同一个平面坐标系下的帧图像拼接，而且还能够解算出拼接后图像对应的高精度有理多项式模型系数，得到高精度的有理多项式模型，为用户进一步利用地面控制点进行高精度几何处理创造了必要条件。

附图说明

图1为本发明实施例的星载线阵传感器摆扫成像示意图；

图2为本发明实施例的光线与高程面和切平面相交示意图；

图3为本发明实施例的原始单帧图像的虚拟立体格网点示意图；

图4为本发明实施例的从原始单帧图像到物方切平面的坐标正算迭代示意图；

图5为本发明实施例的虚拟立体格网点的切平面图像坐标示意图。

图6为本发明实施例的原始图像中心点示意图。

具体实施方式

本发明提出的一种星载线阵传感器垂轨摆扫影像的几何预处理，所要解决的问题是通过对原始单帧图像进行基于物方切平面坐标系的几何纠正与拼接处理，实现对垂轨摆扫影像的高精度几何预处理。以下结合附图和实施例详细说明本发明技术方案。

实施例针对星载线阵传感器垂轨摆扫影像，以参考椭球采用WGS84椭球（涉及地心直角坐标系和地理坐标系），地心惯性坐标系采用J2000惯性坐标系为例，执行步骤如下，可采用计算机软件技术实现自动运行流程：

步骤1，根据原始单帧图像的成像几何，建立原始单帧图像的共线方程模型，实现从原始单帧图像坐标到地心直角坐标的正算。

如图1为星载线阵传感器摆扫成像示意图。图上方粗黑线为线阵CCD，放置方向与飞行方向一致；线阵每次成像得到一个子帧，在一个摆扫、回扫周期内，线阵摆扫（成像）得到一个单帧图像，线阵回扫不成像，仅使得线阵回到初始位置，为下一个周期的摆扫做好准备；线阵摆扫的方向与飞行方向垂直，线阵回扫的方向与摆扫方向相反。图1中线阵摆扫方向从右往左，线阵回扫方向从左往右，具体实施时，线阵摆扫方向也可能从左往右，线阵回扫方向与线阵摆扫方向相反即可。图下方为每个单帧图像的地面覆盖范围，相邻单帧图像的地面覆盖范围有少许重叠，每个单帧图像存在全景畸变，即：由中心向两边，子帧对应的地面范围逐渐变大。

实施例建立原始单帧图像的共线方程模型如下：

假设R_t、R_GF、R_FB、R_BS分别为t时刻下从J2000惯性坐标系到地固地心直角坐标系的旋转矩阵、从轨道坐标系到J2000惯性坐标系的旋转矩阵、从本体坐标系到轨道坐标系的旋转矩阵、从相机坐标系到本体坐标系的旋转矩阵,[X_t Y_t Z_t]^T为t时刻卫星质心在WGS84地心直角坐标系下的坐标矢量，[x(l) y(l) 1]^T为通过像点p(s,l)的光线在相机坐标系下的矢量，则其对应物点P在WGS84地心直角坐标系下的坐标[X_P Y_P Z_P]^T为：

$\left[\begin{array}{c}{X}_P\\ Y_P\\ Z_P\end{array}\right]=\left[\begin{array}{c}{X}_t\\ Y_t\\ Z_t\end{array}\right]+{\mathrm{mR}}_t{R}_{\mathrm{GF}}{R}_{\mathrm{FB}}{R}_{\mathrm{BS}}\left[\begin{array}{c}x\left(l\right)\\ y\left(l\right)\\ 1\end{array}\right]---\left(1\right)$

其中，(s,l)为像点的原始单帧图像坐标，每个原始单帧图像有多个子帧，每个子帧对应一列像素，每列像素个数=CCD线阵的探元个数。l为原始单帧图像的行号，s为原始单帧图像的列号，利用现有技术，可从s计算出成像时刻t；m为比例系数，可利用现有技术，通过光线与参考椭球的某高程面H相交，求解一元二次多项式，得到m的值，再代入上式，求得P点的地心直角坐标[X_P Y_P Z_P]^T。

为便于实施参考起见，提供详细说明如下：

设 ${\left[\begin{array}{ccc}\stackrel{\‾}{X}& \stackrel{\‾}{Y}& \stackrel{\‾}{Z}\end{array}\right]}^T$ 为经过像点p(s,l)的光线在WGS84地心直角坐标系中的矢量表示，其值仅代表了矢量的方向，矢量长度由m决定。

令 $\left[\begin{array}{c}\stackrel{\‾}{X}\\ \stackrel{\‾}{Y}\\ \stackrel{\‾}{Z}\end{array}\right]=R_t{R}_{\mathrm{GF}}{R}_{\mathrm{FB}}{R}_{\mathrm{BS}}\left[\begin{array}{c}x\left(l\right)\\ y\left(l\right)\\ 1\end{array}\right]$

则： $\left[\begin{array}{c}{X}_P\\ Y_P\\ Z_P\end{array}\right]=\left[\begin{array}{c}{X}_t\\ Y_t\\ Z_t\end{array}\right]+m\left[\begin{array}{c}\stackrel{\‾}{X}\\ \stackrel{\‾}{Y}\\ \stackrel{\‾}{Z}\end{array}\right],$ 即 $\begin{array}{c}{X}_P=X_t+m\stackrel{\‾}{X}\\ Y_P=Y_t+m\stackrel{\‾}{Y}\\ Z_p=Z_t+m\stackrel{\‾}{Z}\end{array}$

设物点P的高程为H，则位于长半轴为a+H，短半轴为b+H的参考椭球面上，由椭球面方程得：

$\frac{{X_P}^2+{Y_P}^2}{{(a+H)}^2}+\frac{{Z_P}^2}{{(b+H)}^2}=1$

$\frac{{(X_t+m\stackrel{\‾}{X})}^2+{(Y_t+m\stackrel{\‾}{Y})}^2}{{(a+H)}^2}+\frac{{(Z_t+m\stackrel{\‾}{Z})}^2}{{(b+H)}^2}=1$

$(\frac{{\stackrel{\‾}{X}}^2+{\stackrel{\‾}{Y}}^2}{{(a+H)}^2}+\frac{{\stackrel{\‾}{Z}}^2}{{(b+H)}^2})m^2+(\frac{\stackrel{\‾}{X}{X}_t+\stackrel{\‾}{Y}{Y}_t}{{(a+H)}^2}+\frac{\stackrel{\‾}{Z}{Z}_t}{{(b+H)}^2})2m+(\frac{{X_t}^2+{Y_t}^2}{{(a+H)}^2}+\frac{{Z_t}^2}{{(b+H)}^2}-1)=0$

通过解算一元二次方程可求得两个根，为光线与参考椭球面相交的两个交点，一个位于参考椭球正面，另一个位于参考椭球反面（舍去不用）。

最后将m代入（1）式可求得P点的地心直角坐标。

步骤2，对每个原始单帧图像进行几何纠正处理，包括以下子步骤，

步骤2.1，建立物方局部坐标系和物方切平面坐标系，构建物方局部坐标与地心直角坐标的相互换算关系。

物方局部坐标系O-XYZ定义如下：经过原始图像中心点的光线与物方平均高程面H₀相交的交点为物方局部坐标系的原点O；XY为物方切平面，与参考椭球面相切，Z轴垂直于XY切平面向上为正；X轴沿飞行方向，Y轴垂直于飞行方向，XYZ组成右手系。物方局部坐标系下的坐标记为物方局部坐标。

其中原始图像中心点如图6所示，每个单帧图像各有一个中心点，多个单帧图像构成的整体原始图像，原始图像的中心点如图所示。假设线阵的探元数为N，原始图像包含M个单帧图像（如图中第1、2…M个原始单帧图像），每个单帧图像包含W个子帧（如图中第1、2、3…W个子帧），则单帧图像大小为宽=W、高=N，原始图像大小为宽=W、高=N×M，原始图像中心点就是（W/2，N×M/2）。

所述物方切平面坐标系O-XY定义为，原点O、X轴、Y轴均与物方局部坐标系重合，X轴、Y轴的刻度单位为米；物方切平面坐标系下的坐标记为物方切平面坐标。

假设物方局部坐标系原点的WGS84地心直角坐标为(X₀,Y₀,Z₀)，WGS84地理坐标为(L₀,B₀)，B₀为纬度，L₀为经度，任一物点P的WGS84地心直角坐标为(X_P,Y_P,Z_P)，则P点的物方局部坐标为(x_P,y_P,z_P)：

$\left[\begin{array}{c}{x}_P\\ y_P\\ z_P\end{array}\right]=R_{\mathrm{WGS}84}^{\mathrm{LOC}}\left[\begin{array}{c}{X}_P-X_0\\ Y_P-Y_0\\ Z_P-Z_0\end{array}\right]---\left(2\right)$

上式为从WGS84地心直角坐标(X_P,Y_P,Z_P)到物方局部坐标(x_P,y_P,z_P)的换算公式。

其中，

为从WGS84地心直角坐标系到物方局部坐标系的旋转矩阵，

$R_{\mathrm{WGS}84}^{\mathrm{LOC}}=\left[\begin{array}{ccc}-\sin \left(L_0\right)& \cos \left(L_0\right)& 0\\ -\sin \left(B_0\right)\cos \left(L_0\right)& -\sin \left(B_0\right)\sin \left(L_0\right)& \cos \left(B_0\right)\\ \cos \left(B_0\right)\cos \left(L_0\right)& \cos \left(B_0\right)\sin \left(L_0\right)& \sin \left(B_0\right)\end{array}\right]---\left(3\right)$

同理，从物方局部坐标(x_P,y_P,z_P)到WGS84地心直角坐标(X_P,Y_P,Z_P)的换算公式为：

$\left[\begin{array}{c}{X}_P\\ Y_P\\ Z_P\end{array}\right]=\left[\begin{array}{c}{X}_0\\ Y_0\\ Z_0\end{array}\right]+{\left(R_{\mathrm{WGS}84}^{\mathrm{LOC}}\right)}^T\left[\begin{array}{c}{x}_P\\ y_P\\ z_P\end{array}\right]---\left(4\right)$

步骤2.2，建立切平面图像坐标系，并构建切平面图像坐标与物方局部坐标之间的相互换算关系。

如图2所示，有地形表面、高程面、切平面，各成像光线与高程面、切平面相交，可得到相应交点。

切平面图像坐标系O’-X’Y’定义如下：位于物方切平面上，切平面图像坐标系的X’轴、Y’轴与物方切平面坐标系的X轴、Y轴平行，坐标轴刻度单位为像素，原点O'为物方切平面坐标系下的图像范围的最小外接矩形的左上角点。切平面图像坐标系下的坐标记为切平面图像坐标。

假设切平面图像坐标系原点的物方局部坐标为(x₀,y₀,z₀)，图像分辨率为d，任一物点P的物方局部坐标为(x_P,y_P,z_P)，则其切平面图像坐标为(s_P,l_P)，

从物方局部坐标(x_P,y_P,z_P)到切平面图像坐标(s_P,l_P)的换算公式为：

s_P＝(x_p-x₀)/d，l_P＝(y_p-y₀)/d       （5）

从切平面图像坐标(s_P,l_P)到物方局部坐标(x_P,y_P,z_P)的换算公式为：

x_P＝x₀+s_p×d，y_P＝y₀+l_p×d       （6）

步骤2.3，对每个原始单帧图像划分虚拟立体格网，计算虚拟立体格网点的WGS84地理坐标；并将虚拟立体格网点作为控制点，解算原始单帧图像对应的有理多项式模型系数。

如图3所示，经过原始单帧图像上的光线与各个高程面H₁,H₂,...,H_K相交，得到虚拟立体格网点。

假设星载线阵传感器的探元数为N，每次扫描的子帧数为W，对应地面的高程范围为[H_min,H_max]，则原始单帧图像大小为N行×W列（像素）。假设经过某像点p(s,l)的光线与K个高程面H₁,H₂,...,H_K相交（其中，0≤l≤N-1，0≤s≤W-1，H_min≤H_k≤H_max，1≤k≤K），将得到K个交点P₁,P₂,...,P_K，根据公式（1），可求出这些交点对应的WGS84地心直角坐标[X_P1 Y_P1 Z_P1]^T,[X_P2 Y_P2 Z_P2]^T,...,[X_PK Y_PK Z_PK]^T，进而根据现有技术求得这些交点对应的WGS84地理坐标(L_P1,B_P1),(L_P2,B_P2),...,(L_Pk,B_PK)。按照此过程，经过原始单帧图像上各像点分别得到交点作为虚拟立体格网点，计算出所有虚拟立体格网点的WGS84地理坐标。

根据现有技术，将上述虚拟立体格网点当作控制点，可解算出原始所有单帧图像对应的有理多项式模型系数。

步骤2.4，根据每个原始单帧图像的有理多项式模型，建立每个原始单帧图像坐标与切平面图像坐标之间的正反计算方式。

假设某原始单帧图像对应的有理多项式模型为：

s＝f_x(B,L,H),l＝f_y(B,L,H)       （7）

其中f_x,f_y为有理多项式函数，(B,L,H)为物点的WGS84地理坐标，(s,l)为对应像点的原始单帧图像坐标。

基于有理多项式模型的坐标反算，是已知点的切平面图像坐标求原始单帧图像坐标。首先由切平面图像坐标(s_P,l_P)按照公式（6）求得物方切平面坐标(x_P,y_P)，其中物方切平面坐标(x_P,y_P)加上z_P=0就可得到物方局部坐标(x_P,y_P,z_P)。然后，根据公式（4）将物方局部坐标(x_P,y_P,z_P)转换为地心直角坐标(X_P,Y_P,Z_P)，再计算出WGS84地理坐标，作为(B,L,H)代入公式（7）求出原始单帧图像坐标。

基于有理多项式模型的坐标正算，是已知点的原始单帧图像坐标求切平面图像坐标。具体先假设物点的初始高程值H⁽⁰⁾为平均高程，代入公式（1）求得物点的地心直角坐标[X_P Y_P Z_P]^T，再根据公式（2）求得物方局部坐标，设为(x⁽⁰⁾,y⁽⁰⁾,z⁽⁰⁾)。由于z⁽⁰⁾不一定等于0，即光线与高程面的交点不一定位于物方切平面上，所以坐标正算需要迭代。

令迭代次数k=1，物点的高程值H^(k)＝H^(k-1)-z^(k-1)，再将原始单帧图像坐标和物点的初始高程值代入公式（1）求得物点的地心直角坐标，即新的[X_P Y_P Z_P]^T，再根据公式（2）求得物方局部坐标的第k次迭代结果(x^(k),y^(k),z^(k))。判断若|z^(k)|＜T_Z，T_Z为预设的阈值，迭代终止，输出物方局部坐标，否则令k=k+1，继续迭代。具体实施时，本领域技术人员可根据精度要求自行设置一个正数作为T_Z，例如0.01米。

具体实施时，可采用包括以下子步骤的流程实现坐标正算：

步骤a，设物点的初始高程值H⁽⁰⁾为平均高程，将原始单帧图像坐标和物点的初始高程值代入公式（1）求得物点的地心直角坐标[X_P Y_P Z_P]^T，再根据公式（2）求得物方局部坐标的初始结果(x⁽⁰⁾,y⁽⁰⁾,z⁽⁰⁾)；

步骤b，令迭代次数k=1；

步骤c，令物点的高程值H^(k)＝H^(k-1)-z^(k-1)，再代入公式（1）求得物点的地心直角坐标，即新的[X_P Y_P Z_P]^T，再根据公式（2）求得物方局部坐标的第k次迭代结果(x^(k),y^(k),z^(k))；

步骤d，判断是否|z^(k)|＜T_Z，T_Z为预设的阈值，是则迭代终止，输出当前的物方局部坐标(x^(k),y^(k),z^(k))，否则令k=k+1，返回步骤c继续迭代。

步骤e，根据公式（5）将物方局部坐标换算成切平面图像坐标。

如图4所示为像点从原始单帧图像坐标到切平面图像坐标的正算迭代示意图，经过原始单帧图像上的光线首先与平均高程面相交，得到的交点1偏离物方切平面较远；通过修正高程值，依次得到交点2、交点3，逐渐接近于物方切平面。

步骤2.5，基于步骤2.4根据有理多项式模型建立的正反计算方式，对每个原始单帧图像进行几何纠正，得到物方切平面坐标系下的相应帧图像。

可按照现有技术中的间接法进行纠正，为便于实施参考起见，提供间接法几何纠正流程说明如下：

假设原始单帧图像的坐标系为坐标系①，切平面图像坐标系作为坐标系②，

1)将坐标系①的原始单帧图像四个角点坐标，通过步骤2.4所得坐标正算关系，得到输出切平面图像在坐标系②中的范围；

2)对坐标系②输出切平面图像范围内的每个像素，通过步骤2.4所得坐标反算关系，得到其在坐标系①中的坐标；

3)最后，根据坐标系①中的坐标，在原始单帧图像上进行灰度重采样（灰度重采样为现有技术），给坐标系②上的每个像素赋灰度值。

上述过程中，坐标正算和坐标反算采用步骤2.4所建立的关系。

步骤3，对步骤2.5所得物方切平面坐标系下的所有帧图像进行基于坐标的拼接，得到拼接后的图像，可输出。

此步骤为现有技术。其中，拼接后图像范围为各个帧图像范围的并集的最小外接矩阵。

步骤4，解算拼接后图像对应的有理多项式模型系数。

其中有理多项式模型系数的解算方法为现有技术。有理多项式模型系数的解算需要虚拟立体格网控制点。对于本发明步骤3得到的拼接后图像，虚拟立体格网控制点来自步骤2.3所有单帧图像的虚拟立体格网点（取并集）。其中，所有单帧图像的虚拟立体格网点取并集后，这些虚拟立体格网点的图像坐标需要重新计算，即先根据步骤2.1所得公式（2）、（3）将各虚拟立体格网点的WGS84地理坐标换算成物方局部坐标，再根据步骤2.2所得公式（5）将物方局部坐标换算成切平面图像坐标。

如图5所示，为虚拟立体格网点的切平面图像坐标示意图。

最后，将所有坐标变换后的虚拟立体格网点作为控制点，用现有技术解算出拼接后图像对应的有理多项式模型系数。

通过以上步骤就可以实现星载线阵传感器垂轨摆扫影像的高精度几何预处理。该方法通过将帧图像纠正到物方切平面坐标系，并进行基于坐标的图像拼接，不仅消除了图像内部几何畸变，还给出了拼接后图像的有理多项式模型系数。该方法可以全自动执行，适合星载线阵传感器影像的几何预处理过程。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (8)

1.一种星载线阵传感器垂轨摆扫影像的几何预处理方法，其特征在于，进行以下步骤，

步骤1，根据原始单帧图像的成像几何，建立原始单帧图像的共线方程模型；

步骤2，对每个原始单帧图像进行几何纠正处理，包括以下子步骤，

步骤2.1，建立物方局部坐标系和物方切平面坐标系，构建物方局部坐标与地心直角坐标的相互换算关系；

所述物方局部坐标系O-XYZ定义为，经过原始图像中心点的光线与物方平均高程面H₀相交的交点为物方局部坐标系的原点O；XY为物方切平面，与参考椭球面相切，Z轴垂直于XY切平面向上为正；X轴沿飞行方向，Y轴垂直于飞行方向，XYZ组成右手系；物方局部坐标系下的坐标记为物方局部坐标；

物方切平面坐标系O-XY定义为，物方切平面坐标系的原点O、X轴、Y轴都与物方局部坐标系的原点O、X轴、Y轴重合，X轴、Y轴的刻度单位为米；物方切平面坐标系下的坐标记为物方切平面坐标；

步骤2.2，建立切平面图像坐标系，并构建切平面图像坐标与物方局部坐标之间的相互换算关系；

所述切平面图像坐标系O’-X’Y’定义为，位于物方切平面上，切平面图像坐标系的X’轴、Y’轴与物方切平面坐标系的X轴、Y轴平行，X’轴、Y’轴的刻度单位为像素，原点O’为物方切平面坐标系下的图像范围的最小外接矩形的左上角点；切平面图像坐标系下的坐标记为切平面图像坐标；

步骤2.3，对每个原始单帧图像划分虚拟立体格网，计算虚拟立体格网点的地心直角坐标，并将虚拟立体格网点作为控制点，解算原始单帧图像对应的有理多项式模型系数；

步骤2.4，根据步骤2.3所得每个原始单帧图像的有理多项式模型和步骤2.1、步骤2.2，建立每个原始单帧图像坐标与切平面图像坐标之间的正反计算方式；

步骤2.5，基于步骤2.4根据有理多项式模型建立的正反计算方式，对各原始单帧图像进行几何纠正，得到物方切平面坐标系下的帧图像；

步骤3，对步骤2.5所得物方切平面坐标系下的所有帧图像进行基于坐标的拼接，得到拼接后的图像；

步骤4，解算步骤3所得拼接后图像对应的有理多项式模型系数。

2.根据权利要求1所述星载线阵传感器垂轨摆扫影像的几何预处理方法，其特征在于：步骤1中，建立原始单帧图像的共线方程模型如下，

设R_t、R_GF、R_FB、R_BS分别为t时刻下从地心惯性坐标系到地固地心直角坐标系的旋转矩阵、从轨道坐标系到地心惯性坐标系的旋转矩阵、从本体坐标系到轨道坐标系的旋转矩阵、从相机坐标系到本体坐标系的旋转矩阵,[X_t Y_t Z_t]^T为t时刻卫星质心在地心直角坐标系下的坐标矢量，[x(l) y(l) 1]^T为通过像点p(s,l)的光线在相机坐标系下的矢量，则其对应物点P在地心直角坐标系下的坐标[X_P Y_P Z_P]^T为：

$\left[\begin{array}{c}{X}_P\\ Y_P\\ Z_P\end{array}\right]=\left[\begin{array}{c}{X}_t\\ Y_t\\ Z_t\end{array}\right]+{\mathrm{mR}}_t{R}_{\mathrm{GF}}{R}_{\mathrm{FB}}{R}_{\mathrm{BS}}\left[\begin{array}{c}x\left(l\right)\\ y\left(l\right)\\ 1\end{array}\right]---\left(1\right)$

其中，(s,l)为像点的原始单帧图像坐标，l为原始单帧图像的行号，s为原始单帧图像的列号。

3.根据权利要求2所述星载线阵传感器垂轨摆扫影像的几何预处理方法，其特征在于：步骤2.1中，构建物方局部坐标与地心直角坐标的相互换算关系如下，

设物方局部坐标系原点的地心直角坐标为(X₀,Y₀,Z₀)，地理坐标为(L₀,B₀)，B₀为纬度，L₀为经度，任一物点P的地心直角坐标为(X_P,Y_P,Z_P)，物点P的物方局部坐标为(x_P,y_P,z_P)，

从地心直角坐标(X_P,Y_P,Z_P)到物方局部坐标(x_P,y_P,z_P)的换算公式为，

$\left[\begin{array}{c}{x}_P\\ y_P\\ z_P\end{array}\right]=R_{\mathrm{WGS}84}^{\mathrm{LOC}}\left[\begin{array}{c}{X}_P-X_0\\ Y_P-Y_0\\ Z_P-Z_0\end{array}\right]---\left(2\right)$

其中， $R_{\mathrm{WGS}84}^{\mathrm{LOC}}=\left[\begin{array}{ccc}-\sin \left(L_0\right)& \cos \left(L_0\right)& 0\\ -\sin \left(B_0\right)\cos \left(L_0\right)& -\sin \left(B_0\right)\sin \left(L_0\right)& \cos \left(B_0\right)\\ \cos \left(B_0\right)\cos \left(L_0\right)& \cos \left(B_0\right)\sin \left(L_0\right)& \sin \left(B_0\right)\end{array}\right]---\left(3\right)$

从物方局部坐标(x_P,y_P,z_P)到地心直角坐标(X_P,Y_P,Z_P)的换算公式为，

$\left[\begin{array}{c}{X}_P\\ Y_P\\ Z_P\end{array}\right]=\left[\begin{array}{c}{X}_0\\ Y_0\\ Z_0\end{array}\right]+{\left(R_{\mathrm{WGS}84}^{\mathrm{LOC}}\right)}^T\left[\begin{array}{c}{x}_P\\ y_P\\ z_P\end{array}\right]---\left(4\right)$

其中，

为从地心直角坐标系到物方局部坐标系的旋转矩阵。

4.根据权利要求3所述星载线阵传感器垂轨摆扫影像的几何预处理方法，其特征在于：步骤2.2中，构建切平面图像坐标与物方局部坐标之间的相互换算关系如下，

设切平面图像坐标系原点的物方局部坐标为(x₀,y₀,z₀)，任一物点P的物方局部坐标为(x_P,y_P,z_P)，相应物方切平面图像坐标为(s_P,l_P)，

从物方局部坐标(x_P,y_P,z_P)到切平面图像坐标(s_P,l_P)的换算公式为：

s_P＝(x_p-x₀)/d，l_P＝(y_p-y₀)/d       （5）

从切平面图像坐标(s_P,l_P)到物方局部坐标(x_P,y_P,z_P)的换算公式为：

x_P＝x₀+s_p×d，y_P＝y₀+l_p×d       （6）

其中，d为图像分辨率。

5.根据权利要求4所述星载线阵传感器垂轨摆扫影像的几何预处理方法，其特征在于：步骤2.3中，对某原始单帧图像划分虚拟立体格网，计算虚拟立体格网点的地心直角坐标实现方式如下，

设星载线阵传感器的探元数为N，每次扫描的子帧数为W，对应地面的高程范围为[H_min,H_max]，则原始单帧图像大小为N行×W列；假设经过某像点p(s,l)的光线与K个高程面H₁,H₂,...,H_K相交，(s,l)为像点的原始单帧图像坐标，0≤l≤N-1，0≤s≤W-1，H_min≤H_k≤H_max，1≤k≤K，得到K个交点P₁,P₂,...,P_K，根据公式（1）求出这些交点对应的地心直角坐标[X_P1 Y_P1 Z_P1]^T,[X_P2 Y_P2 Z_P2]^T,...,[X_PK Y_PK Z_PK]^T，进而求得这些交点对应的地理坐标(L_P1,B_P1),(L_P2,B_P2),...,(L_Pk,B_PK)。

6.根据权利要求5所述星载线阵传感器垂轨摆扫影像的几何预处理方法，其特征在于：步骤2.4中，根据有理多项式模型，建立每个原始单帧图像坐标与切平面图像坐标之间的正反计算方式如下，

设原始单帧图像对应的有理多项式模型为：

s＝f_x(B,L,H),l＝f_y(B,L,H)      （7）

其中f_x,f_y为有理多项式函数，(B,L,H)为物点的地理坐标，(s,l)为对应像点的原始单帧图像坐标；

基于有理多项式模型的坐标反算，包括首先由切平面图像坐标按照公式（6）求得物方切平面坐标(x_P,y_P)，物方切平面坐标(x_P,y_P)加上z_P=0得到物方局部坐标(x_P,y_P,z_P)；然后，先根据公式（4）将物方局部坐标(x_P,y_P,z_P)转换为地心直角坐标(X_P,Y_P,Z_P)，再计算出地理坐标作为(B,L,H)代入公式（7）求出原始单帧图像坐标；

基于有理多项式模型的坐标正算，包括以下子步骤，

步骤a，设物点的初始高程值H⁽⁰⁾为平均高程，将原始单帧图像坐标和物点的初始高程值代入公式（1）求得物点的地心直角坐标[X_P Y_P Z_P]^T，再根据公式（2）求得物方局部坐标的初始结果(x⁽⁰⁾,y⁽⁰⁾,z⁽⁰⁾)；

步骤b，令迭代次数k=1；

步骤c，令物点的高程值H^(k)＝H^(k-1)-z^(k-1)，再代入公式（1）求得物点的地心直角坐标作为新的[X_P Y_P Z_P]^T，再根据公式（2）求得物方局部坐标的第k次迭代结果(x^(k),y^(k),z^(k))；

步骤d，判断是否|z^(k)|＜T_Z，T_Z为预设的阈值，是则迭代终止，输出当前的物方局部坐标(x^(k),y^(k),z^(k))，否则令k=k+1，返回步骤c继续迭代；

步骤e，根据公式（5）将物方局部坐标换算成切平面图像坐标。

7.根据权利要求6所述星载线阵传感器垂轨摆扫影像的几何预处理方法，其特征在于：步骤4中，解算步骤3所得拼接后图像对应的有理多项式模型系数实现如下，

首先，对步骤2.3所得各原始单帧图像的虚拟立体格网点取并集后，重新计算这些虚拟立体格网点的图像坐标，包括先根据公式（2）、（3）将各虚拟立体格网点的地理坐标换算成物方局部坐标，再根据所得公式（5）将物方局部坐标换算成切平面图像坐标；

然后，将所有坐标变换后的虚拟立体格网点作为控制点，解算出拼接后图像对应的有理多项式模型系数。

8.根据权利要求1或2或3或4或5或6或7所述星载线阵传感器垂轨摆扫影像的几何预处理方法，其特征在于：所述地心直角坐标和地理坐标都是定义在以地球质心为原点的参考椭球上的，以地球质心为原点的参考椭球为WGS84椭球。

Images