[发明专利]一种稳健的卫星遥感影像有理函数模型参数估计方法

摘要

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

主权项

一种稳健的卫星遥感影像有理函数模型参数估计方法，其步骤如下：第一步，建立控制点格网和检查点格网控制点格网和检查点格网同时建立，相比于控制点格网，检查点格网中的经纬度间隔取其一半，高程分层数取其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模型定义关系式如下：X=NumL(P,L,H)DenL(P,L,H)Y=Nums(P,L,H)Dens(P,L,H)---(3)]]>(P,L,H)为正则化的地面坐标，(X,Y)为正则化的影像坐标，计算公式如下：P=Lat-LAT_OFFLAT_SCALEL=Lon-LONG_OFFLONG_SCALEH=Height-HEIGHT_OFFHEIGHT_SCALEX=Sample-SAMP_OFFSAMP_SCALEY=Line-LINE_OFFLINE_SCALE---(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_OFF、LAT_SCALE计算为例，首先遍历控制点格网中每个格网点的纬度值，统计控制点格网纬度最大值Latmax、最小值Latmin和平均值Latavg，则正则化参数LAT_OFF、LAT_SCALE可按公式(5)得到，其他正则化参数可按与之相同方法计算得到；LAT_OFF=LatavgLAT_SCALE=max{Latmax-Latavg,Latavg-Latmin}---(5)]]>NumL(P,L,H)、DenL(P,L,H)、Nums(P,L,H)和Dens(P,L,H)为三次多项式，形式如下：NumL(P,L,H)=a1+a2L+a3P+a4H+a5LP+a6LH+a7PH+a8L2+a9P2+a10H2+a11PLH+a12L3+a13LP2+a14LH2+a15L2P+a16P3+a17PH2+a18L2H+a19P2H+a20H3DenL(P,L,H)=b1+b2L+b3P+b4H+b5LP+b6LH+b7PH+b8L2+b9P2+b10H2+b11PLH+b12L3+b13LP2+b14LH2+b15L2P+b16P3+b17PH2+b18L2H+b19P2H+b20H3Nums(P,L,H)=c1+c2L+c3P+c4H+c5LP+c6LH+c7PH+c8L2+c9P2+c10H2+c11PLH+c12L3+c13LP2+c14LH2+c15L2P+c16P3+c17PH2+c18L2H+c19P2H+c20H3Dens(P,L,H)=d1+d2L+d3P+d4H+d5LP+d6LH+d7PH+d8L2+d9P2+d10H2+d11PLH+d12L3+d13LP2+d14LH2+d15L2P+d16P3+d17PH2+d18L2H+d19P2H+d20H3---(6)]]>式(6)中，b1和d1取为1，三次多项式的系数即称为待求的RFM模型参数，共78个；将公式(3)变形为：FX=Nums(P,L,H)-X*Dens(P,L,H)=0FY=NumL(P,L,H)-Y*DenL(P,L,H)=0---(7)]]>以RFM模型参数为未知数，对公式(7)进行线性化，则误差方程为：V＝Ax‑l   (8)式中，A=∂FX∂ai∂FX∂bj∂FX∂ci∂FX∂dj∂FY∂ai∂FY∂bj∂FY∂ci∂FY∂dj,i=1,20,j=2,20,]]>l=-FX0-FY0,]]>x＝[ai bj ci dj]T；矩阵A和向量l即所需的误差方程系数矩阵和常数向量；(2)、建立新的参数求解模型：采用Householder变换方法对误差方程系数矩阵A∈Rm×n进行QR分解，有：A=QR0---(9)]]>式中，Q∈Rm×m为正交矩阵，R∈Rn×n为上三角矩阵；对正交矩阵Q分块为Q＝[Q1 Q2]，则则由此可知Rx＝c1的最小二乘解即为公式(8)的最小二乘解，新参数求解模型为：V′＝Rx‑c1    (10)式中，R为误差方程系数矩阵A进行QR分解后得到的上三角矩阵，c1为正交矩阵Q与原始常数向量l相乘后取得前n个数值所构成的新常数向量；第三步，采用Levenberg‑Marquardt方法进行RFM模型参数求解，步骤如下：(1)、参数初始值的确定：对于RFM模型参数pk，b1和d1始终取值1，在首次迭代中，仅取NumL(P,L,H)、DenL(P,L,H)、Nums(P,L,H)和Dens(P,L,H)的一次项系数为未知数，即a1～a4、b2～b4、c1～c4和d2～d4，共14个参数，对公式(7)进行线性化，并按最小二乘方法求解，获取a1～a4、b2～b4、c1～c4和d2～d4的初始值，对于RFM模型中的高阶项系数初始值均直接取0；在随后的迭代中，对RFM模型的78个参数均进行求解，参数初始值取为前一次迭代平差所计算的新参数值；对于Levenberg‑Marquardt算法中所使用的变量，k为迭代次数，λk为阻尼系数，v为阻尼系数调整常数；(2)、利用新的参数求解模型计算矩阵H：H＝R+λkI   (11)并构造增量正规方程H·δk＝c1   (12)式中，δk表示待解算的RFM模型参数pk的增量向量，c1表示利用未知数当前值所计算的新参数求解模型的误差向量，利用公式(12)进行求解，可以得到δk；(3)、阻尼系数调整：利用δk更新参数向量并计算新的误差向量，当误差向量模值大于给定收敛阈值时，将误差向量模值同上一次迭代结果进行比较，如果误差向量模值变小，令pk+1＝pk+δk，λk+1＝λk/v，并转到第三步的步骤(2)中进行迭代；如果误差向量模值变大，令λk+1＝λk·v，并返回到第三步的步骤(1)重新解算；当误差向量模值小于给定收敛阈值时，满足收敛条件，则停止迭代，输出结果；第四步，利用检查点格网对RFM模型进行精度评定利用第三步中解算出的RFM模型参数，对检查点格网中的每个格网点，通过公式(1)和公式(2)，计算格网点物方坐标(Lat,Lon)和Height所对应的影像坐标(Samplecal,Linecal)，而利用严格成像几何模型计算所获取的影像坐标(Sample,Line)可视为真值，计算(Samplecal,Linecal)与(Sample,Line)的差值，并统计出影像行列方向上最大误差、最小误差和中误差，从而对RFM模型拟合严格成像几何模型的精度进行评价。