CN101907705A - 通用的多源遥感影像几何校正模型联合平差方法 - Google Patents

Abstract

本发明从遥感影像的正算几何模型出发，建立了一种通用的多源遥感影像成像模型的区域网平差，并采用改进的Levenberg-Marquardt算法解算所建立的平差方程，克服了模型病态性对解的影响。该方法与遥感影像的具体成像模型无关，可以对不同的成像模型进行联合平差。在多源遥感数据的区域网平差过程中，可以针对各数据源的特点选择最合适的成像模型，充分发挥各种模型在严密性、定位准确性、保密性等方面的特点。

Description

通用的多源遥感影像几何校正模型联合平差方法

技术领域

本发明建立了一种通用的多源遥感影像成像模型的区域网平差，该方法可以对不同的成像模型进行联合平差。在多源遥感数据的区域网平差过程中，可以针对各数据源的特点选择最合适的成像模型，充分发挥各种模型在严密性、定位准确性、保密性等方面的特点。

技术背景

卫星遥感影像的成像几何模型是描述地面点三维空间坐标与相应像点在像平面坐标系的二维坐标之间的映射关系的数学模型，它是进行影像几何校正和目标定位的基础。随着航天技术、传感器技术的发展，对于卫星遥感影像，已经出现了各种各样的成像模型，它们在严密性、复杂性及定位准确性方面都有着各自的特点。目前最常用的几何模型有严格成像模型、多项式模型、仿射变换模型和有理函数模型等，对于不同的卫星遥感数据源，严格成像模型又具有各自不同的表达形式。

对于各种成像模型，现在已有相应的区域网平差方法，但这些方法或者是针对具体的成像几何模型进行的，或者是基于有理函数模型进行的。针对具体成像模型的区域网平差方法缺乏通用性，不能用于多源遥感影像的联合平差；基于有理函数模型进行的区域网平差方法利用了有理函数模型的通用性，可以用于多源遥感影像的联合平差，但是对于一些严格成像模型已知的数据源(如Landsat系列、Spot系列等)，有理函数模型并不是最佳的成像模型，这种联合平差的方法就不能发挥不同成像几何模型的优势。本文将针对遥感影像成像几何模型建立一种统一的区域网平差方法，实现不同成像模型的联合平差。

发明内容

本发明将不同的遥感影像成像几何模型描述为通用的正算模型，建立了一种统一的区域网平差方法，实现了不同成像模型的联合平差。

按照地面点三维空间坐标与相应像点二维平面坐标之间的映射顺序，卫星遥感影像的成像几何模型可以分为两类：反算模型和正算模型。正算模型表示为：

$\left\{\begin{array}{c}x=f_x(t,X,Y,Z)\\ y=f_y(t,X,Y,Z)\end{array}\right.---(1-1)$

其中，(X，Y，Z)表示控制点的地面坐标，(x，y)表示控制点在影像上的量测坐标，t＝(t₁，t₂，…t_n)^T表示传感器几何校正模型的参数。

不同卫星的严格成像模型、仿射变换模型、多项式模型和有理函数模型等均可表示成通用的正算模型，因此可以将多种成像模型联合起来进行区域网平差。模型线性化后，每个控制点的误差方程为：

$\left\{\begin{array}{c}{-v}_x=\frac{{\∂f}_x}{\∂t}\mathrm{dt}+\frac{{\∂f}_x}{\∂X}\mathrm{dX}+\frac{{\∂f}_x}{\∂Y}\mathrm{dY}+\frac{\∂f_x}{\∂Z}\mathrm{dZ}-l_x\\ {-v}_y=\frac{{\∂f}_y}{\∂t}\mathrm{dt}+\frac{{\∂f}_y}{\∂X}\mathrm{dX}+\frac{{\∂f}_y}{\∂Y}\mathrm{dY}+\frac{\∂f_y}{\∂Z}\mathrm{dZ}-l_y\end{array}\right.---(1-2)$

其中， $\left\{\begin{array}{c}{l}_x=x-f_x(t_0,X_0,Y_0,Z_0)\\ l_y=y-f_y(t_0,X_0,Y_0,Z_0)\end{array}\right.$

将(1.2)表示成矩阵形式：

V＝Adt+Bds-L                            (1-3)

其中， $A=\left[\begin{array}{cccc}\frac{{\∂f}_x}{{\mathrm{dt}}_1}& \frac{{\∂f}_x}{{\mathrm{dt}}_2}& \·\·\·& \frac{{\∂f}_x}{{\mathrm{dt}}_n}\\ \frac{{\∂f}_y}{{\mathrm{dt}}_1}& \frac{{\∂f}_y}{{\mathrm{dt}}_2}& \·\·\·& \frac{{\mathrm{df}}_y}{{\mathrm{dt}}_n}\end{array}\right],$ $B=\left[\begin{array}{ccc}\frac{{\∂f}_x}{\mathrm{dX}}& \frac{{\∂f}_x}{\mathrm{dY}}& \frac{{\∂f}_x}{\mathrm{dZ}}\\ \frac{{\∂f}_y}{\mathrm{dX}}& \frac{{\∂f}_y}{\mathrm{dY}}& \frac{{\∂f}_y}{\mathrm{dZ}}\end{array}\right],$

dt＝(dt₁，dt₂，…，dt_n)^T，ds＝(dX，dY，dZ)^T，V＝(v_x，v_y)^T，L＝(l_x，l_y)^T。

这里没有求解函数偏导的具体解析形式，而是采用数值计算的方法来逼近各函数偏导，如：

$\frac{{\∂f}_x}{\∂X}=\underset{\mathrm{\ΔX}\→0}{\lim }\frac{f_x(t,X+\mathrm{\ΔX},Y,Z)}{\mathrm{\ΔX}}---(1-4)$

现有m景影像(可以为多源数据)，n个控制点，其中n_c个已知控制点，n_u个未知控制点(即连接点)，所有控制点所对应的像平面点个数为n′，设第l景影像的模型参数个数为p_l。

单个控制点的线性化误差方程可以表示为：

V_i＝A_idt_i+B_ids_i-L_i                      (1-5)

其中，

${\mathrm{dt}}_i={({\mathrm{dt}}_1^i,\·\·\·,{\mathrm{dt}}_{p_i}^i)}^T$

ds_i＝(dX_i，dY_i，dZ_i)^T

$V_i={(v_x^i,v_y^i)}^T$

$L_i={(l_x^i,l_y^i)}^T$

1)若第i个像平面点在第j景影像上，对应于第k个未知控制点，则：

2)若第i个像平面点在第j景影像上，对应于一个已知控制点，则：

联立所有控制点的误差方程，多模型的区域网平差方程可以记为：

V＝Adt+Bds-L            (1-6)

其中，

$A=\left[\begin{array}{c}{A}_1\\ A_2\\ \·\\ \·\\ \·\\ A_{n^{\′}}\end{array}\right],$ $B=\left[\begin{array}{c}{B}_1\\ B_2\\ \·\\ \·\\ \·\\ B_{n^{\′}}\end{array}\right],$ $\mathrm{dt}=\left[\begin{array}{c}{\mathrm{dt}}_1\\ {\mathrm{dt}}_2\\ \·\\ \·\\ \·\\ {\mathrm{dt}}_m\end{array}\right],$ $\mathrm{ds}=\left[\begin{array}{c}{\mathrm{ds}}_1\\ {\mathrm{ds}}_2\\ \·\\ \·\\ \·\\ {\mathrm{ds}}_{n_u}\end{array}\right],$ $V=\left[\begin{array}{c}{V}_1\\ V_2\\ \·\\ \·\\ \·\\ V_{n^{\′}}\end{array}\right],$ $L=\left[\begin{array}{c}{L}_1\\ L_2\\ \·\\ \·\\ \·\\ L_{n^{\′}}\end{array}\right]$

A为

的矩阵，B为2n′×3n_u的矩阵。

多模型区域网平差方程(1-6)的解算：首先采用改化法方程的方法，消去多模型的区域网平差方程中所有连接点的地面坐标未知数，仅保留各传感器几何模型的未知参数，然后利用改进的LM(Levenberg-Marquardt)算法求解出各未知参数。

在传统的LM算法中，第k次迭代计算的未知数改变量ΔX为：

${\mathrm{\ΔX}}_{\left(k\right)}=-{(J_{\left(k\right)}^T{J}_{\left(k\right)}+{\mathrm{\μ}}_{\left(k\right)}I)}^{-1}{J}_{\left(k\right)}^{T}V---(1-7)$

式中I是单位矩阵，J_(k)是系数矩阵的Jacobian矩阵，μ_(k)＞0为阻尼系数。

对LM算法的改进如下：

设的对角线元素组成向量

阻尼矩阵为

其中 ${\mathrm{\μ}}_k=\max (q_1^k,q_2^k,\·\·\·,q_{p_1}^k)\×0.001$

则(1-7)式改为：

${\mathrm{\ΔX}}_{\left(k\right)}=-{(J_{\left(k\right)}^T{J}_{\left(k\right)}+D_{\left(k\right)})}^{-1}{J}_{\left(k\right)}^{T}V---(1-8)$

附图说明

图1是使用本发明方法的具体实例处理过程示意图。

具体实施过程

下文结合说明书附图1，以Landsat-5和CBERS卫星遥感影像数据为例，对本发明的具体实施方式作详细说明。其中Landsat-5数据采用严格成像模型，CBERS数据采用有理函数模型。本发明所述方法包含但不限于所举实例。

步骤1：读取已经标准化后的卫星遥感影像数据；

步骤2：输入已知控制点和未知控制点(及连接点)；

步骤3：利用卫星影像参数获取模型参数的初始值，并计算未知控制点地面坐标的初始值，当前迭代次数置为0；

步骤4：若当前迭代次数小于最大迭代次数，转至步骤5；若当前迭代次数达到最大迭代次数，转至步骤9；

步骤5：依次对输入的卫星遥感影像数据建立误差方程，模型的线性化过程中采用数值计算的方法逼近模型对各参数的偏导；

步骤6：将所有卫星遥感影像数据的误差方程联立为整体的区域网平差方程，采用改化法方程的方法，消去其中所有连接点的地面坐标未知数，仅保留各传感器几何模型的未知参数；

步骤7：使用改进的LM算法求出未知参数的改正值向量，并更新参数的值；

步骤8：用求得的改正值向量的模与规定的限差比较，若小于限差，转至步骤9；否则，当前迭代次数增加1，转至步骤4；

步骤9：解算完成，输出各影像模型参数。

Claims (4)

1.通用的多源遥感影像几何校正模型联合平差方法的主要步骤是：

步骤(1)：读取已经标准化后的卫星遥感影像数据；

步骤(2)：输入已知控制点和未知控制点(及连接点)；

步骤(3)：利用卫星影像参数获取模型参数的初始值，并计算未知控制点地面坐标的初始值，当前迭代次数置为0；

步骤(4)：若当前迭代次数小于最大迭代次数，转至步骤5；若当前迭代次数达到最大迭代次数，转至步骤9；

步骤(5)：依次对输入的卫星遥感影像数据建立误差方程，模型的线性化过程中采用数值计算的方法逼近模型对各参数的偏导；

步骤(6)：将所有卫星遥感影像数据的误差方程联立为整体的区域网平差方程，采用改化法方程的方法，消去其中所有连接点的地面坐标未知数，仅保留各传感器几何模型的未知参数；

步骤(7)：使用改进的LM算法出未知参数的改正值向量，并更新参数的值；

步骤(8)：用求得的改正值向量的模与规定的限差比较，若小于限差，转至步骤9；否则，当前迭代次数增加1，转至步骤4；

步骤(9)：解算完成，输出各影像模型参数。

2.根据权利要求1所示方法步骤(5)，其特征是：将不同的遥感影像成像几何模型描述为通用的正算模型，采用数值计算的方法来逼近各函数偏导，是各种模型的误差方程形式一致。

正算模型表示为：

$\left\{\begin{array}{c}x=f_x(t,X,Y,Z)\\ y=f_y(t,X,Y,Z)\end{array}\right.---(2-1)$

模型线性化后，每个控制点的误差方程为：

$\left\{\begin{array}{c}{-v}_x=\frac{{\∂f}_x}{\∂t}\mathrm{dt}+\frac{{\∂f}_x}{\∂X}\mathrm{dX}+\frac{{\∂f}_x}{\∂Y}\mathrm{dY}+\frac{\∂f_x}{\∂Z}\mathrm{dZ}-l_x\\ {-v}_y=\frac{{\∂f}_y}{\∂t}\mathrm{dt}+\frac{{\∂f}_y}{\∂X}\mathrm{dX}+\frac{{\∂f}_y}{\∂Y}\mathrm{dY}+\frac{\∂f_y}{\∂Z}\mathrm{dZ}-l_y\end{array}\right.---(2-2)$

以X为例计算函数的偏导，如(2-3)式所示：

$\frac{{\∂f}_x}{\∂X}=\underset{\mathrm{\ΔX}\→0}{\lim }\frac{f_x(t,X+\mathrm{\ΔX},Y,Z)}{\mathrm{\ΔX}}---(2-3)$

3.根据权利要求1所示方法步骤(6)，将各误差方程联立成整体的区域网平差方程：

V＝Adt+Bds-L                (3-1)

式中，

$A=\left[\begin{array}{c}{A}_1\\ A_2\\ \·\\ \·\\ \·\\ A_{n^{\′}}\end{array}\right],$ $B=\left[\begin{array}{c}{B}_1\\ B_2\\ \·\\ \·\\ \·\\ B_{n^{\′}}\end{array}\right],$ $\mathrm{dt}=\left[\begin{array}{c}{\mathrm{dt}}_1\\ {\mathrm{dt}}_2\\ \·\\ \·\\ \·\\ {\mathrm{dt}}_m\end{array}\right],$ $\mathrm{ds}=\left[\begin{array}{c}{\mathrm{ds}}_1\\ {\mathrm{ds}}_2\\ \·\\ \·\\ \·\\ {\mathrm{ds}}_{n_u}\end{array}\right],$ $V=\left[\begin{array}{c}{V}_1\\ V_2\\ \·\\ \·\\ \·\\ V_{n^{\′}}\end{array}\right],$ $L=\left[\begin{array}{c}{L}_1\\ L_2\\ \·\\ \·\\ \·\\ L_{n^{\′}}\end{array}\right]$

${\mathrm{dt}}_i={({\mathrm{dt}}_1^i,\·\·\·,{\mathrm{dt}}_{p_i}^i)}^T$

ds_i＝(dX_i，dY_i，dZ_i)^T

$V_i={(v_x^i,v_y^i)}^T$

$L_i={(l_x^i,l_y^i)}^T$

1)若第i个像平面点在第j景影像上，对应于第k个未知控制点，则：

2)若第i个像平面点在第j景影像上，对应于一个已知控制点，则：

4.根据权利要求1所示，方法步骤(7)，对LM算法的改进如下：

在传统的LM算法中，第k次迭代计算的未知数改变量ΔX为：

${\mathrm{\ΔX}}_{\left(k\right)}=-{(J_{\left(k\right)}^T{J}_{\left(k\right)}+{\mathrm{\μ}}_{\left(k\right)}I)}^{-1}{J}_{\left(k\right)}^{T}V---(4-1)$

式中I是单位矩阵，J_(k)是系数矩阵的Jacobian矩阵，μ_(k)＞0为阻尼系数。

设

的对角线元素组成向量

阻尼矩阵为

其中 ${\mathrm{\μ}}_k=\max (q_1^k,q_2^k,\·\·\·,q_{p_1}^k)\×0.001$

将(4-1)式改为：

${\mathrm{\ΔX}}_{\left(k\right)}=-{(J_{\left(k\right)}^T{J}_{\left(k\right)}+D_{\left(k\right)})}^{-1}{J}_{\left(k\right)}^{T}V---(4-2)$

Images