CN104123457B - 一种稳健的卫星遥感影像有理函数模型参数估计方法 - Google Patents

Abstract

本发明属于测绘科学技术领域，具体公开了一种稳健的卫星遥感影像有理函数模型参数估计方法，其步骤：用严格成像几何模型和全球DEM在影像覆盖范围内建立控制点和检查点格网；以RFM模型参数为未知数对RFM模型线性化得到误差方程系数矩阵，并用Householder变换的方法对误差方程系数矩阵进行QR分解，建立新的参数解算方程；用Levenberg‑Marquardt方法对该方程进行参数求解得到所需RFM模型参数；用解算出的RFM模型计算检查点格网物方坐标所对应的影像坐标，并与检查点格网中原有影像坐标比较，统计影像行列方向上最大、最小和中误差，对模型拟合精度评定。该方法可获得比岭估计算法更优的拟合精度。

Description

一种稳健的卫星遥感影像有理函数模型参数估计方法

技术领域

本发明涉及测绘科学技术领域，具体涉及一种稳健的卫星遥感影像有理函数模型参数估计方法，该方法可以有效克服模型参数求解的病态问题，实现光学、SAR等卫星遥感影像有理函数模型参数的稳健、无偏和精确估计。

背景技术

随着对地观测技术的飞速发展，光学、SAR等不同类型传感器不断涌现。为了实现卫星遥感影像的高精度信息提取，建立卫星遥感影像的精确几何定位模型至关重要。有理函数模型(Rational Function Model，RFM)作为一种通用的传感器成像模型，形式简单，计算速度快，更具有一般性和保密性。目前，许多遥感卫星运营商在提供数据时均提供了相应的RFM参数。研究遥感影像RFM模型参数稳健无偏估计方法，将有利于建立不同类型卫星遥感影像统一的几何定位模型，降低后继数据处理难度，具有重要的意义。

然而，RFM模型参数估计在理论研究和实际应用中主要凸显出如下问题：

①RFM模型参数求解存在病态问题；

②高阶RFM模型因参数过多导致解不稳定；

③采用RFM模型拟合严格成像几何模型时会存在较为明显的残余误差。

目前，国内外学者主要从算法改进和模型改造两个方面进行研究。如采用岭估计、L曲线法等方法改善法方程的状态求解有理多项式系数，或略去RFM三次项通过减少参数个数来改善法方程的病态。但是，岭估计、L曲线法不能彻底解决过度参数化问题，且属于有偏估计。仅仅采用RFM二次项拟合严格成像模型，拟合精度往往难以保证。为此，人们进一步提出了谱修正迭代算法、基于离差阵逐步回归的RFM参数筛选求解方法。

总体而言，目前估计RFM模型参数的方法均基于系数矩阵法方程进行求解，不可避免的存在如下缺陷：法方程矩阵制约数为误差方程系数矩阵制约数的平方，且阶数越高，制约数越大。三阶RFM模型，待求解的参数将达到78个，这也是高阶RFM模型参数求解存在病态、导致解不稳定的主要原因。

近年来，随着计算机技术的迅猛发展，微机内存不断扩展，速度亦在不断提升，运用正交分解技术直接基于误差方程式解算RFM模型参数，将有利于提高参数求解的稳定性。为此，发明人融入正交分解技术、稳健最优化估计技术，提出了一种直接基于误差方程系数矩阵的RFM模型参数无偏稳健估计方法，以解决RFM模型参数的稳健无偏估计问题。

发明内容

本发明针对目前卫星遥感影像RFM模型参数解算过程中，存在法方程病态、影响解算精度等缺点，提出了一种稳健的卫星遥感影像有理函数模型参数估计方法。该方法可直接基于误差方程系数矩阵，大大提高了RFM模型参数估计的稳健性和精度。

为了达到上述目的，本发明采用如下技术措施：

(1)建立控制点格网和检查点格网。首先确定影像的近似覆盖范围，并依据全球DEM确定高程变化范围，然后按照控制点和检查点给定的经纬度间隔、高程分层数建立控制点立体和检查点立体格网，计算格网点物方坐标对应的影像坐标，从而建立控制点格网和检查点格网。

(2)建立直接基于误差方程系数矩阵的参数求解模型。依据RFM模型，首先计算RFM模型误差方程的系数矩阵和常数向量，然后采用Householder变换方法对误差方程系数矩阵进行QR分解，利用QR分解的结果建立新的参数求解模型。

(3)采用Levenberg-Marquardt方法进行参数求解。首先确定RFM模型参数的初始值，然后采用Levenberg-Marquardt方法进行参数求解，并不断对阻尼系数进行调整，直至满足收敛条件。

(4)精度评定。通过建立的检查点格网，利用求出的RFM模型计算格网点对应的影像坐标，将检查点格网中该格网点的真实影像坐标进行比较，统计出影像行列方向上最大误差、最小误差和中误差，进行精度评定。

与现有的技术相比，本发明的优点和显著效果主要表现在：

1、直接基于误差方程系数矩阵求解RFM模型参数，可以有效减少方程解算的病态性；

2、可以稳健无偏解算出RFM模型参数，获得更好的拟合精度；

3、无需任何人工干预，可自动计算出所需的78个RFM模型参数；

4、本发明为卫星遥感影像的RFM模型参数解算提供了一种稳健可靠的计算方法，创造性地直接基于误差方程系数矩阵求解，可以有效克服基于法方程矩阵求解存在病态、解不稳定的缺陷，实现了RFM模型参数的稳健、无偏、高精度估计。

具体实施方式

下面对本发明做进一步详细描述。

实施例1：

以一景SAR卫星遥感影像的RFM模型参数求解为例，采用与地形无关的计算方法，各步骤详细阐述如下：

第一步，建立控制点格网和检查点格网

控制点格网和检查点格网一般同时建立，一般情况下，相比于控制点格网，检查点格网中的经纬度间隔取其一半，高程分层数取其2倍。

具体步骤如下：

(1)假设SAR卫星遥感影像严格成像几何模型的正算公式如下：

(Lat，Lon)＝T(Sample，Line，Height) (1)

式中，T表示由影像坐标(Sample，Line)和大地高Height计算大地经纬度坐标(Lat，Lon)的转换关系。利用影像的左上角和右下角的像点坐标，通过正算公式(1)计算对应的大地经纬度坐标，取最小包络矩形为影像的近似覆盖范围。

(2)根据影像的近似覆盖范围，以给定的控制点和检查点经纬度间隔，分别建立控制点平面格网和检查点平面格网。

(3)利用全球DEM，如美国航空航天局提供的SRTM C波段数据、德国宇航中心提供的SRTM X波段数据等，根据控制点平面格网中格网点的大地经纬度坐标，内插出各格网点的高程，统计这些高程值中的最大值和最小值，得到影像覆盖范围内的高程变化范围，然后分别根据给定的高程分层数建立控制点立体格网和检查点立体格网。

(4)假设SAR卫星遥感影像严格成像几何模型的反算公式如下：

(Sample，Line)＝T^-1(Lat，Lon，Height) (2)

式中，T^-1表示由大地经纬度坐标(Lat，Lon)和大地高Height计算影像坐标(Sample，Line)的转换关系。通过反算公式(2)，可以计算出上一步骤所建立的控制点立体格网和检查点立体格网中每个格网点所对应的影像坐标，从而建立物方均匀分布的控制点格网和检查点格网。

第二步，建立直接基于误差方程系数矩阵的参数求解模型

具体步骤如下：

(1)、计算RFM模型误差方程的系数矩阵和常数向量。RFM模型定义关系式如下：

(P，L，H)为正则化的地面坐标，(X，Y)为正则化的影像坐标，计算公式如下：

式中，LAT_OFF、LAT_SCALE、LONG_OFF、LONG_SCALE、HEIGHT_OFF和HEIGHT_SCALE为地面点坐标的正则化参数。SAMP_OFF、SAMP_SCALE、LINE_OFF和LINE_SCALE为影像坐标的正则化参数。(Lat，Lon)和Height分别表示控制点格网中某格网点所对应的大地经纬度坐标和大地高，统称为格网点物方坐标；(Sample，Line)为该格网点所对应的影像坐标。

正则化参数的计算方法如下，以大地坐标纬度的正则化参数LAT_OFF、LAT_SCALE计算为例，首先遍历控制点格网中每个格网点的纬度值，统计控制点格网纬度最大值Lat_max、最小值Lat_min和平均值Lat_avg，则正则化参数LAT_OFF、LAT_SCALE可按公式(5)得到，其他正则化参数可按与之相同方法计算得到。

Num_L(P，L，H)、Den_L(P，L，H)、Num_s(P，L，H)和Den_s(P，L，H)为三次多项式，形式如下：

式(6)中，b₁和d₁取为1，三次多项式的系数即称为待求的RFM模型参数，共78个。

将公式(3)变形为：

F_X＝Num_s(P，L，H)-X*Den_s(P，L，H)＝0

(7)

F_Y＝Num_L(P，L，H)-Y*Den_L(P，L，H)＝0

以RFM模型参数为未知数，对公式(7)进行线性化，则误差方程为：

V＝Ax-1 (8)

式中，

x＝[a_i b_j c_i d_j]^T。

矩阵A和向量l即我们所需的误差方程系数矩阵和常数向量。

(2)、建立新的参数求解模型。采用Householder变换方法对误差方程系数矩阵A∈R^m×n进行QR分解，有：

式中，Q∈R^m×m为正交矩阵，R∈R^n×n为上三角矩阵。

对正交矩阵Q分块为Q＝[Q₁ Q₂]，则

则由此可知Rx＝c₁的最

小二乘解即为公式(8)的最小二乘解，新参数求解模型为：

V′＝Rx-c₁ (10)

式中，R为误差方程系数矩阵A进行QR分解后得到的上三角矩阵，c₁为正交矩阵Q与原始常数向量l相乘后取得前n个数值所构成的新常数向量。

第三步，采用Levenberg-Marquardt方法进行RFM模型参数求解

Levenberg-Marquardt算法，是数学中最常见的优化算法，是一种介于牛顿法与梯度下降法之间的一种非线性优化算法，其优点在于对过参数化问题不敏感，能够有效处理冗余参数带来的强相关问题，使迭代优化陷于局部极小值的机会大大减少。

利用Levenberg-Marquardt算法RFM模型参数求解的步骤如下：

(1)、参数初始值的确定：对于RFM模型参数p_k，b₁和d₁始终取值1，在首次迭代中，仅取Num_L(P，L，H)、Den_L(P，L，H)、Num_s(P，L，H)和Den_s(P，L，H)的一次项系数为未知数，即a₁～a₄、b₂～b₄、c₁～c₄和d₂～d₄，共14个参数，对公式(7)进行线性化，并按最小二乘方法求解，获取a₁～a₄、b₂～b₄、c₁～c₄和d₂～d₄的初始值，对于RFM模型中的高阶项系数初始值均直接取0。在随后的迭代中，对RFM模型的78个参数均进行求解，参数初始值取为前一次迭代平差所计算的新参数值。对于Levenberg-Marquardt算法中所使用的变量，k为迭代次数，λ_k为阻尼系数，v为阻尼系数调整常数，本实施例中k初始化为0，λ_k初始值λ₀设为0.01，v取值为10。

(2)、利用新的参数求解模型计算矩阵H：

H＝R+λ_kI (11)

并构造增量正规方程

H·δ_k＝c₁ (12)

式中，δ_k表示待解算的RFM模型参数p_k的增量向量，c₁表示利用未知数当前值所计算的新参数求解模型的误差向量。利用公式(12)进行求解，可以得到δ_k。

(3)、阻尼系数调整：利用δ_k更新参数向量并计算新的误差向量，当误差向量模值大于给定收敛阈值(本实施例中取1e-8)，进一步将误差向量模值同上一次迭代结果进行比较，如果误差向量模值变小，令p_k+1＝p_k+δ_k，λ_k+1＝λ_k/v，并转到第三步的步骤(2)中进行迭代；如果误差向量模值变大，令λ_k+1＝λ_k·v，并返回到第三步的步骤(1)重新解算。当误差向量模值小于给定收敛阈值，满足收敛条件，则停止迭代，输出结果。

第四步，利用检查点格网对RFM模型进行精度评定

利用第三步中解算出的RFM模型参数，对检查点格网中的每个格网点，通过公式(1)和公式(2)，计算格网点物方坐标(Lat，Lon)和Height所对应的影像坐标(Sample_cal，Line_cal)，而利用严格成像几何模型计算所获取的影像坐标(Sample，Line)可视为真值，计算(Sample_cal，Line_cal)与(Sample，Line)的差值，并统计出影像行列方向上最大误差、最小误差和中误差，从而可以对RFM模型拟合严格成像几何模型的精度进行评价。

下面以一景EnvisatSAR卫星遥感影像的RFM模型参数计算为例，将本实施例1提出的方法与经典的岭估计方法进行对比试验。在RFM模型的实际运用过程中，为了降低RFM模型参数求解的病态性，常用的一种简化方法是令Den_L(P，L，H)＝Den_s(P，L，H)≠1，此时RFM模型参数个数减少为59个。下面将分别针对78参数RFM模型和59参数RFM模型进行对比试验，结果如表1所示。

表1 本实施例1的方法与岭估计算法的对比

从表中数据容易发现，当RFM模型参数减少为59个参数时，岭估计算法与实施例1方法的拟合精度相当，平面拟合误差都达到0.008像素左右。但是，岭估计算法受岭参数的影响较大，当使用相同的岭参数对78参数RFM模型进行计算时，岭估计算法检查点的平面拟合误差达到12.7992像素，计算失败。而与此同时，本专实施例方法仍然可以成功解算出78参数RFM模型，并且获得比59参数RFM模型更优的拟合结果，平面拟合误差降低为0.006像素左右。

算法效率方面，在Intel(R)Xeon(R)CPU X5550@2.67GHz 2.66GHz(2处理器)上运行，岭估计算法耗时18.02s，本实施例提出的算法为19.59s，耗时上只是略微提高，但是免除了选取合适岭参数的环节，完全可满足实际应用的需求。

Claims (1)

1.一种稳健的卫星遥感影像有理函数模型参数估计方法，其步骤如下：

第一步，建立控制点格网和检查点格网

控制点格网和检查点格网同时建立，相比于控制点格网，检查点格网中的经纬度间隔取其一半，高程分层数取其2倍；

具体步骤如下：

(1)假设SAR卫星遥感影像严格成像几何模型的正算公式如下：

(Lat,Lon)＝T(Sample,Line,Height) (1)

式中，T表示由影像坐标(Sample,Line)和大地高Height计算大地经纬度坐标(Lat,Lon)的转换关系，利用影像的左上角和右下角的像点坐标，通过正算公式(1)计算对应的大地经纬度坐标，取最小包络矩形为影像的近似覆盖范围；

(2)根据影像的近似覆盖范围，以给定的控制点和检查点经纬度间隔，分别建立控制点平面格网和检查点平面格网；

(3)利用全球DEM，根据控制点平面格网中格网点的大地经纬度坐标，内插出各格网点的高程，统计这些高程值中的最大值和最小值，得到影像覆盖范围内的高程变化范围，然后分别根据给定的高程分层数建立控制点立体格网和检查点立体格网；

(4)假设SAR卫星遥感影像严格成像几何模型的反算公式如下：

(Sample,Line)＝T^-1(Lat,Lon,Height) (2)

式中，T^-1表示由大地经纬度坐标(Lat,Lon)和大地高Height计算影像坐标(Sample,Line)的转换关系，通过反算公式(2)，可以计算出上一步骤所建立的控制点立体格网和检查点立体格网中每个格网点所对应的影像坐标，从而建立物方均匀分布的控制点格网和检查点格网；

第二步，建立直接基于误差方程系数矩阵的参数求解模型

(1)、计算RFM模型误差方程的系数矩阵和常数向量，RFM模型定义关系式如下：

$\left\{\begin{array}{c}X=\frac{{\mathrm{Num}}_L(P,L,H)}{{\mathrm{Den}}_L(P,L,H)}\\ Y=\frac{{\mathrm{Num}}_s(P,L,H)}{{\mathrm{Den}}_s(P,L,H)}\end{array}\right.---\left(3\right)$

(P,L,H)为正则化的地面坐标，(X,Y)为正则化的影像坐标，计算公式如下：

$\begin{array}{c}P=\frac{Lat-LAT\_OFF}{LAT\_SCALE}\\ L=\frac{Lon-LONG\_OFF}{LONG\_SCALE}\\ H=\frac{Height-HEIGHT\_OFF}{HEIGHT\_SCALE}\\ X=\frac{Sample-SAMP\_OFF}{SAMP\_SCALE}\\ Y=\frac{Line-LINE\_OFF}{LINE\_SCALE}\end{array}---\left(4\right)$

式中，LAT_OFF、LAT_SCALE、LONG_OFF、LONG_SCALE、HEIGHT_OFF和HEIGHT_SCALE为地面点坐标的正则化参数；SAMP_OFF、SAMP_SCALE、LINE_OFF和LINE_SCALE为影像坐标的正则化参数；(Lat,Lon)和Height分别表示控制点格网中某格网点所对应的大地经纬度坐标和大地高，统称为格网点物方坐标；(Sample,Line)为该格网点所对应的影像坐标；

正则化参数的计算方法如下，以大地坐标纬度的正则化参数LAT_OFF、LAT_SCALE计算为例，首先遍历控制点格网中每个格网点的纬度值，统计控制点格网纬度最大值Lat_max、最小值Lat_min和平均值Lat_avg，则正则化参数LAT_OFF、LAT_SCALE可按公式(5)得到，其他正则化参数可按与之相同方法计算得到；

$\begin{array}{c}LAT\_OFF={\mathrm{Lat}}_{avg}\\ LAT\_SCALE=max\{{\mathrm{Lat}}_{max}-{\mathrm{Lat}}_{avg},{\mathrm{Lat}}_{avg}-{\mathrm{Lat}}_{min}\}\end{array}---\left(5\right)$

Num_L(P,L,H)、Den_L(P,L,H)、Num_s(P,L,H)和Den_s(P,L,H)为三次多项式，形式如下：

$\begin{array}{c}{\mathrm{Num}}_L(P,L,H)=a_1+a_{2}L+a_{3}P+a_{4}H+a_{5}LP+a_{6}LH+a_{7}PH+a_8{L}^2+a_9{P}^2\\ +a_{10}{H}^2+a_{11}PLH+a_{12}{L}^3+a_{13}{\mathrm{LP}}^2+a_{14}{\mathrm{LH}}^2+a_{15}{L}^{2}P+a_{16}{P}^3+a_{17}{\mathrm{PH}}^2\\ +a_{18}{L}^{2}H+a_{19}{P}^{2}H+a_{20}{H}^3\\ {\mathrm{Den}}_L(P,L,H)=b_1+b_{2}L+b_{3}P+b_{4}H+b_{5}LP+b_{6}LH+b_{7}PH+b_8{L}^2+b_9{P}^2\\ +b_{10}{H}^2+b_{11}PLH+b_{12}{L}^3+b_{13}{\mathrm{LP}}^2+b_{14}{\mathrm{LH}}^2+b_{15}{L}^{2}P+b_{16}{P}^3+b_{17}{\mathrm{PH}}^2\\ +b_{18}{L}^{2}H+b_{19}{P}^{2}H+b_{20}{H}^3\\ {\mathrm{Num}}_s(P,L,H)=c_1+c_{2}L+c_{3}P+c_{4}H+c_{5}LP+c_{6}LH+c_{7}PH+c_8{L}^2+c_9{P}^2\\ +c_{10}{H}^2+c_{11}PLH+c_{12}{L}^3+c_{13}{\mathrm{LP}}^2+c_{14}{\mathrm{LH}}^2+c_{15}{L}^{2}P+c_{16}{P}^3+c_{17}{\mathrm{PH}}^2\\ +c_{18}{L}^{2}H+c_{19}{P}^{2}H+c_{20}{H}^3\\ {\mathrm{Den}}_s(P,L,H)=d_1+d_{2}L+d_{3}P+d_{4}H+d_{5}LP+d_{6}LH+d_{7}PH+d_8{L}^2+d_9{P}^2\\ +d_{10}{H}^2+d_{11}PLH+d_{12}{L}^3+d_{13}{\mathrm{LP}}^2+d_{14}{\mathrm{LH}}^2+d_{15}{L}^{2}P+d_{16}{P}^3+d_{17}{\mathrm{PH}}^2\\ +d_{18}{L}^{2}H+d_{19}{P}^{2}H+d_{20}{H}^3\end{array}---\left(6\right)$

式(6)中，b₁和d₁取为1，三次多项式的系数即称为待求的RFM模型参数，

共78个；

将公式(3)变形为：

$\begin{array}{c}{F}_X={\mathrm{Num}}_s(P,L,H)-X*{\mathrm{Den}}_s(P,L,H)=0\\ F_Y={\mathrm{Num}}_L(P,L,H)-Y*{\mathrm{Den}}_L(P,L,H)=0\end{array}---\left(7\right)$

以RFM模型参数为未知数，对公式(7)进行线性化，则误差方程为：

V＝Ax-l (8)

式中，

$A=\left[\begin{array}{cccc}\frac{\∂F_X}{\∂a_i}& \frac{\∂F_X}{\∂b_j}& \frac{\∂F_X}{\∂c_i}& \frac{\∂F_X}{\∂d_j}\\ \frac{\∂F_Y}{\∂a_i}& \frac{\∂F_Y}{\∂b_j}& \frac{\∂F_Y}{\∂c_i}& \frac{\∂F_Y}{\∂d_j}\end{array}\right],i=1,20,j=2,20,$

$l=\left[\begin{array}{c}-F_X^0\\ -F_Y^0\end{array}\right],$

x＝[a_i b_j c_i d_j]^T；

矩阵A和向量l即所需的误差方程系数矩阵和常数向量；

(2)、建立新的参数求解模型：采用Householder变换方法对误差方程系数矩阵A∈R^m×n进行QR分解，有：

$A=Q\left[\begin{array}{c}R\\ 0\end{array}\right]---\left(9\right)$

式中，Q∈R^m×m为正交矩阵，R∈R^n×n为上三角矩阵；

对正交矩阵Q分块为Q＝[Q₁ Q₂]，则

则由此可知Rx＝c₁的最小二乘解即为公式(8)的最小二乘解，新参数求解模型为：

V′＝Rx-c₁ (10)

式中，R为误差方程系数矩阵A进行QR分解后得到的上三角矩阵，c₁为正交矩阵Q与原始常数向量l相乘后取得前n个数值所构成的新常数向量；

第三步，采用Levenberg-Marquardt方法进行RFM模型参数求解，步骤如下：(1)、参数初始值的确定：对于RFM模型参数p_k，b₁和d₁始终取值1，在首次迭代中，仅取Num_L(P,L,H)、Den_L(P,L,H)、Num_s(P,L,H)和Den_s(P,L,H)的一次项系数为未知数，即a₁～a₄、b₂～b₄、c₁～c₄和d₂～d₄，共14个参数，对公式(7)进行线性化，并按最小二乘方法求解，获取a₁～a₄、b₂～b₄、c₁～c₄和d₂～d₄的初始值，对于RFM模型中的高阶项系数初始值均直接取0；在随后的迭代中，对RFM模型的78个参数均进行求解，参数初始值取为前一次迭代平差所计算的新参数值；对于Levenberg-Marquardt算法中所使用的变量，k为迭代次数，λ_k为阻尼系数，v为阻尼系数调整常数；

(2)、利用新的参数求解模型计算矩阵H：

H＝R+λ_kI (11)

并构造增量正规方程

H·δ_k＝c₁ (12)

式中，δ_k表示待解算的RFM模型参数p_k的增量向量，c₁表示利用未知数当前值所计算的新参数求解模型的误差向量，利用公式(12)进行求解，可以得到δ_k；

(3)、阻尼系数调整：利用δ_k更新参数向量并计算新的误差向量，当误差向量模值大于给定收敛阈值时，将误差向量模值同上一次迭代结果进行比较，如果误差向量模值变小，令p_k+1＝p_k+δ_k，λ_k+1＝λ_k/v，并转到第三步的步骤(2)中进行迭代；如果误差向量模值变大，令λ_k+1＝λ_k·v，并返回到第三步的步骤(1)重新解算；当误差向量模值小于给定收敛阈值时，满足收敛条件，则停止迭代，输出结果；

第四步，利用检查点格网对RFM模型进行精度评定

利用第三步中解算出的RFM模型参数，对检查点格网中的每个格网点，通过公式(1)和公式(2)，计算格网点物方坐标(Lat,Lon)和Height所对应的影像坐标(Sample_cal,Line_cal)，而利用严格成像几何模型计算所获取的影像坐标(Sample,Line)可视为真值，计算(Sample_cal,Line_cal)与(Sample,Line)的差值，并统计出影像行列方向上最大误差、最小误差和中误差，从而对RFM模型拟合严格成像几何模型的精度进行评价。