一种稳健的卫星遥感影像有理函数模型参数估计方法
【专利摘要】本发明属于测绘科学【技术领域】，具体公开了一种稳健的卫星遥感影像有理函数模型参数估计方法，其步骤：用严格成像几何模型和全球DEM在影像覆盖范围内建立控制点和检查点格网；以RFM模型参数为未知数对RFM模型线性化得到误差方程系数矩阵，并用Householder变换的方法对误差方程系数矩阵进行QR分解，建立新的参数解算方程；用Levenberg-Marquardt方法对该方程进行参数求解得到所需RFM模型参数；用解算出的RFM模型计算检查点格网物方坐标所对应的影像坐标，并与检查点格网中原有影像坐标比较，统计影像行列方向上最大、最小和中误差，对模型拟合精度评定。该方法可获得比岭估计算法更优的拟合精度。
【专利说明】一种稳健的卫星遥感影像有理函数模型参数估计方法
【技术领域】
[0001] 本发明涉及测绘科学【技术领域】，具体涉及一种稳健的卫星遥感影像有理函数模型 参数估计方法，该方法可以有效克服模型参数求解的病态问题，实现光学、SAR等卫星遥感 影像有理函数模型参数的稳健、无偏和精确估计。
【背景技术】
[0002] 随着对地观测技术的飞速发展，光学、SAR等不同类型传感器不断涌现。为了实现 卫星遥感影像的高精度信息提取，建立卫星遥感影像的精确几何定位模型至关重要。有理 函数模型（Rational Function Model，RFM)作为一种通用的传感器成像模型，形式简单，计 算速度快，更具有一般性和保密性。目前，许多遥感卫星运营商在提供数据时均提供了相应 的RFM参数。研究遥感影像RFM模型参数稳健无偏估计方法，将有利于建立不同类型卫星 遥感影像统一的几何定位模型，降低后继数据处理难度，具有重要的意义。
[0003] 然而，RFM模型参数估计在理论研究和实际应用中主要凸显出如下问题：
[0004] ① RFM模型参数求解存在病态问题；
[0005] ②高阶RFM模型因参数过多导致解不稳定；
[0006] ③采用RFM模型拟合严格成像几何模型时会存在较为明显的残余误差。
[0007] 目前，国内外学者主要从算法改进和模型改造两个方面进行研究。如采用岭估计、 L曲线法等方法改善法方程的状态求解有理多项式系数，或略去RFM三次项通过减少参数 个数来改善法方程的病态。但是，岭估计、L曲线法不能彻底解决过度参数化问题，且属于 有偏估计。仅仅采用RFM二次项拟合严格成像模型，拟合精度往往难以保证。为此，人们进 一步提出了谱修正迭代算法、基于离差阵逐步回归的RFM参数筛选求解方法。
[0008] 总体而言，目前估计RFM模型参数的方法均基于系数矩阵法方程进行求解，不可 避免的存在如下缺陷：法方程矩阵制约数为误差方程系数矩阵制约数的平方，且阶数越高， 制约数越大。三阶RFM模型，待求解的参数将达到78个，这也是高阶RFM模型参数求解存 在病态、导致解不稳定的主要原因。
[0009] 近年来，随着计算机技术的迅猛发展，微机内存不断扩展，速度亦在不断提升，运 用正交分解技术直接基于误差方程式解算RFM模型参数，将有利于提高参数求解的稳定 性。为此，发明人融入正交分解技术、稳健最优化估计技术，提出了一种直接基于误差方程 系数矩阵的RFM模型参数无偏稳健估计方法，以解决RFM模型参数的稳健无偏估计问题。
【发明内容】
[0010] 本发明针对目前卫星遥感影像RFM模型参数解算过程中，存在法方程病态、影响 解算精度等缺点，提出了一种稳健的卫星遥感影像有理函数模型参数估计方法。该方法可 直接基于误差方程系数矩阵，大大提高了 RFM模型参数估计的稳健性和精度。
[0011] 为了达到上述目的，本发明采用如下技术措施：
[0012] (1)建立控制点格网和检查点格网。首先确定影像的近似覆盖范围，并依据全球 DEM确定高程变化范围，然后按照控制点和检查点给定的经纬度间隔、高程分层数建立控制 点立体和检查点立体格网，计算格网点物方坐标对应的影像坐标，从而建立控制点格网和 检查点格网。
[0013] (2)建立直接基于误差方程系数矩阵的参数求解模型。依据RFM模型，首先计算 RFM模型误差方程的系数矩阵和常数向量，然后采用Householder变换方法对误差方程系 数矩阵进行QR分解，利用QR分解的结果建立新的参数求解模型。
[0014] (3)采用Levenberg-Marquardt方法进行参数求解。首先确定RFM模型参数的初 始值，然后采用Levenberg-Marquardt方法进行参数求解，并不断对阻尼系数进行调整，直 至满足收敛条件。
[0015] (4)精度评定。通过建立的检查点格网，利用求出的RFM模型计算格网点对应的影 像坐标，将检查点格网中该格网点的真实影像坐标进行比较，统计出影像行列方向上最大 误差、最小误差和中误差，进行精度评定。
[0016] 与现有的技术相比，本发明的优点和显著效果主要表现在：
[0017] 1、直接基于误差方程系数矩阵求解RFM模型参数，可以有效减少方程解算的病态 性；
[0018] 2、可以稳健无偏解算出RFM模型参数，获得更好的拟合精度；
[0019] 3、无需任何人工干预，可自动计算出所需的78个RFM模型参数；
[0020] 4、本发明为卫星遥感影像的RFM模型参数解算提供了一种稳健可靠的计算方法， 创造性地直接基于误差方程系数矩阵求解，可以有效克服基于法方程矩阵求解存在病态、 解不稳定的缺陷，实现了 RFM模型参数的稳健、无偏、高精度估计。
【具体实施方式】
[0021] 下面对本发明做进一步详细描述。
[0022] 实施例1 :
[0023] 以一景SAR卫星遥感影像的RFM模型参数求解为例，采用与地形无关的计算方法， 各步骤详细阐述如下 ：
[0024] 第一步，建立控制点格网和检查点格网
[0025] 控制点格网和检查点格网一般同时建立，一般情况下，相比于控制点格网，检查点 格网中的经纬度间隔取其一半，高程分层数取其2倍。
[0026] 具体步骤如下：
[0027] (1)假设SAR卫星遥感影像严格成像几何模型的正算公式如下：
[0028] (Lat, Lon) = T(Sample, Line, Height) (1)
[0029] 式中，T表示由影像坐标（Sample, Line)和大地高Height计算大地经讳度坐标 (Lat，Lon)的转换关系。利用影像的左上角和右下角的像点坐标，通过正算公式（1)计算 对应的大地经纬度坐标，取最小包络矩形为影像的近似覆盖范围。
[0030] (2)根据影像的近似覆盖范围，以给定的控制点和检查点经纬度间隔，分别建立控 制点平面格网和检查点平面格网。
[0031] (3)利用全球DEM，如美国航空航天局提供的SRTM C波段数据、德国宇航中心提供 的SRTM X波段数据等，根据控制点平面格网中格网点的大地经纬度坐标，内插出各格网点 的高程，统计这些高程值中的最大值和最小值，得到影像覆盖范围内的高程变化范围，然后 分别根据给定的高程分层数建立控制点立体格网和检查点立体格网。
[0032] (4)假设SAR卫星遥感影像严格成像几何模型的反算公式如下：
[0033] (Sample, Line) = T_1(Lat, Lon, Height) (2)
[0034] 式中，Γ1表示由大地经纬度坐标（Lat，Lon)和大地高Height计算影像坐标 (Sampl e，Line)的转换关系。通过反算公式（2)，可以计算出上一步骤所建立的控制点立体 格网和检查点立体格网中每个格网点所对应的影像坐标，从而建立物方均匀分布的控制点 格网和检查点格网。
[0035] 第二步，建立直接基于误差方程系数矩阵的参数求解模型
[0036] 具体步骤如下：
[0037] (1)、计算RFM模型误差方程的系数矩阵和常数向量。RFM模型定义关系式如下：
[0038]
【权利要求】
1. 一种稳健的卫星遥感影像有理函数模型参数估计方法，其步骤如下： (1) 建立控制点格网和检查点格网：首先确定影像的近似覆盖范围，并依据全球DEM确 定高程变化范围，然后按照控制点和检查点给定的经纬度间隔、高程分层数建立控制点立 体和检查点立体格网，计算格网点物方坐标对应的影像坐标，从而建立控制点格网和检查 点格网； (2) 建立直接基于误差方程系数矩阵的参数求解模型：依据RFM模型，首先计算RFM模 型误差方程的系数矩阵和常数向量，然后采用Householder变换方法对误差方程系数矩阵 进行QR分解，利用QR分解的结果建立新的参数求解模型； (3) 采用Levenberg-Marquardt方法进行参数求解：首先确定RFM模型参数的初始值， 然后采用Levenberg-Marquardt方法进行参数求解，并不断对阻尼系数进行调整，直至满 足收敛条件； (4) 精度评定：通过建立的检查点格网，利用求出的RFM模型计算格网点对应的影像坐 标，与检查点格网中该格网点的真实影像坐标进行比较，统计出影像行列方向上最大误差、 最小误差和中误差，进行精度评定。
2. 根据权利要求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) 式中，Γ1表示由大地经纟韦度坐标（Lat, Lon)和大地高Height计算影像坐标（Sample, Line)的转换关系，通过反算公式（2)，可以计算出上一步骤所建立的控制点立体格网和检 查点立体格网中每个格网点所对应的影像坐标，从而建立物方均匀分布的控制点格网和检 查点格网； 第二步，建立直接基于误差方程系数矩阵的参数求解模型 (1)、计算RFM模型误差方程的系数矩阵和常数向量，RFM模型定义关系式如下： (3) (P，L，Η)为正则化的地面坐标，（X，Y)为正则化的影像坐标，计算公式如下：
(4) 式中，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_0FF、LAT_SCALE计算 为例，首先遍历控制点格网中每个格网点的纬度值，统计控制点格网纬度最大值Latmax、最 小值Lat min和平均值Latavg，则正则化参数LAT_0FF、LAT_SCALE可按公式（5)得到，其他正 则化参数可按与之相同方法计算得到； LAT_0FF = Latavg (5) LAT-SCALE - ni￡ix {L&tmax_L￡itaVg，L&taVg _L￡itmin} NumL(P，L，H)、DenL(P，L，H)、Nums(P，L，H)和 Dens(P，L，H)为三次多项式，形式如下： NumL(P, L, Η) = a1+a2L+a3P+a4H+a5LP+a 6LH+a7PH+a8L2+agP2 +a10H2+anPLH+a12L3+a 13LP2+a14LH2+a15L 2P+a16P3+a17PH2 +a18L2H+a19P2H+a2〇H 3 DenL(P, L, H) = b1+b2L+b3P+b4H+b5LP+b 6LH+b7PH+b8L2+b9P 2 +b10H2+bnPLH+b12L3+b 13LP2+b14LH2+b15L 2P+b16P3+b17PH2 (6) +b18L2H+b19P2H+b20H 3 Nums(P, L, H) = c1+c2L+c3P+C4H+c5LP+c 6LH+c7PH+c8L2+CgP2 +c10H2+cnPLH+c12L3+c 13LP2+c14LH2+c15L 2P+c16P3+c17PH2 +c18L2H+c19P2H+c2〇H 3 Dens(P, L, H) = d1+d2L+d3P+d4H+d5LP+d 6LH+d7PH+d8L2+d9P 2 +d10H2+d11PLH+d12L3+d 13LP2+d14LH2+d15L 2P+d16P3+d17PH2 +d18L2H+d19P2H+d20H 3 式（6)中，^和屯取为1，三次多项式的系数即称为待求的RFM模型参数，共78个； 将公式（3)变形为： Fx = Nums (P，L，H) -X*Dens (P，L，Η) = Ο (7) FY = NumL(P，L，H)-Y*DenL(P，L，Η) = Ο 以RFM模型参数为未知数，对公式（7)进行线性化，则误差方程为： V = Ax-1 (8) 式中，
(i = 1,20, j = 2,20),
x = [a-j bj Cj dj]； 矩阵A和向量1即所需的误差方程系数矩阵和常数向量； (2)、建立新的参数求解模型：采用Householder变换方法对误差方程系数矩阵 A e RmXn进行QR分解，有：
(9) 式中，Q e ΓΧπ为正交矩阵，R e RnXn为上三角矩阵； 对正交矩阵Q分块为Q = [Qi Q2]，则
则
由此可知Rx = Cl的最小二乘解即为公 式（8)的最小二乘解，新参数求解模型为： Y' = Rx-Cl (10) 式中，R为误差方程系数矩阵A进行QR分解后得到的上三角矩阵，Cl为正交矩阵Q与 原始常数向量1相乘后取得前η个数值所构成的新常数向量； 第三步，采用Levenberg-Marquardt方法进行RFM模型参数求解，步骤如下：（1)、参数 初始值的确定：对于RFM模型参数pk，h和屯始终取值1，在首次迭代中，仅取NumJP，L， Η)、DenJP，L，Η)、Nums(P，L，Η)和 Dens(P，L，Η)的一次项系数为未知数，SP ?a4、b2 ? b4、Cl?c4和d2?d4，共14个参数，对公式（7)进行线性化，并按最小二乘方法求解，获取 at?a4、b2?bpq?c4和d2?d 4的初始值，对于RFM模型中的高阶项系数初始值均直接 取〇 ;在随后的迭代中，对RFM模型的78个参数均进行求解，参数初始值取为前一次迭代平 差所计算的新参数值；对于Levenberg-Marquardt算法中所使用的变量，k为迭代次数，λ k 为阻尼系数，v为阻尼系数调整常数； (2)、利用新的参数求解模型计算矩阵Η : H = R+AkI (11) 并构造增量正规方程 Η · δ k = Cl (12) 式中，S k表示待解算的RFM模型参数pk的增量向量，Cl表示利用未知数当前值所计算 的新参数求解模型的误差向量，利用公式（12)进行求解，可以得到Sk; (3)、阻尼系数调整：利用δ k更新参数向量并计算新的误差向量，当误差向量模值大于 给定收敛阈值时，将误差向量模值同上一次迭代结果进行比较，如果误差向量模值变小，令 Pk+1 = Pk+ δ k，λ k+1 = λ k/v，并转到第三步的步骤（2)中进行迭代；如果误差向量模值变大， 令A k+1 = λ, · v，并返回到第三步的步骤（1)重新解算；当误差向量模值小于给定收敛阈 值时，满足收敛条件，则停止迭代，输出结果； 第四步，利用检查点格网对RFM模型进行精度评定 利用第三步中解算出的RFM模型参数，对检查点格网中的每个格网点，通过公式（1) 和公式（2)，计算格网点物方坐标（Lat，Lon)和Height所对应的影像坐标（Sample。# Line。」，而利用严格成像几何模型计算所获取的影像坐标（Sample，Line)可视为真值，计 算（SampleupLine^)与（Sample，Line)的差值，并统计出影像行列方向上最大误差、最小 误差和中误差，从而对RFM模型拟合严格成像几何模型的精度进行评价。
【文档编号】G06F19/00GK104123457SQ201410348719
【公开日】2014年10月29日 申请日期:2014年7月21日 优先权日:2014年7月21日
【发明者】吴颖丹, 郑列, 朱永松, 朱莹 申请人:湖北工业大学