基于1SUB1/SUB范数约束的有理函数模型参数求解和全参数优化方法.pdf

基于l1范数约束的有理函数模型参数求解和全参数优化方法：在地形无关方案中，当卫星影像内部畸变较大时，仍可以得到可靠的有理函数模型参数；在地形相关方案中，可利用少于39个控制点求解78个有理函数模型系数；在已知有理函数模型系数的条件下，利用少量地面控制点对有理函数模型进行全参数优化。本发明提高了有理函数模型参数求解的稳定性，减少了所需控制点的个数，可以在少于39个控制点的条件下求出有理函数模型参数的可靠解，并提高了影像精校正的精度。

1.  一种基于l₁范数约束的有理函数模型求解方法，所述方法利用n个控制点进行求解，控制点可以是真实地面控制点，也可以是通过严格成像模型生成的虚拟控制点；其特征是：包括以下步骤，
步骤1.1，对n个控制点的像平面坐标和地面坐标进行标准化；
步骤1.2，通过变形将有理函数模型(1)转化为关于各RPCs的线性模型(2)，
l=Nl(X,Y,Z)Ds(X,Y,Z)s=Nl(X,Y,Z)Ds(X,Y,Z)---(1)]]>
Nl(X,Y,Z)-lDl(X,Y,Z)=0Ns(X,Y,Z)-sDs(X,Y,Z)=0---(2)]]>
然后n个控制点所列的误差方程可以矩阵形式来表示，如式(3)所示，
L＝Ax+V            (3)
其中
L＝[l₁，l₂，…，l_n，s₁，s₂，…，s_n]^T，
x＝[a₀，…，a₁₉，b₁，…b₁₉，c₀，…，c₁₉，d₁，…，d₁₉]^T，
V为随机误差，
A＝[A₁，A₂，…，A_n]^T，
Ai=[1,Xi,Yi,...,Zi3,-liXi,-liYi,...,-liZi3,1,Xi,Yi,...,Zi3,-siXi,-siYi,...,-siZi3]T,]]>
(l_i，s_i)为第i个控制点标准化后的像平面行坐标和列坐标，
(X_i，Y_i，Z_i)为第i个控制点标准化后的地面经度、纬度和高程，
i＝1，2，…，n；
步骤1.3，对模型(3)加入l₁范数约束，将问题转化为带l₁范数约束的最优化问题：
min||Ax-L||22s.t.||x||1≤α---(4)]]>
其中α为l₁范数约束因子；
步骤1.4，利用LARS算法求解得到未知向量x的估计值
步骤1.5，根据估计值得到有理函数模型参数RPCs，生成有理函数模型。

2.  一种对有理函数模型全参数进行优化的精纠正方法，所述方法中已有的有理函数模型参数从卫星影像自带的RPC文件中获取或通过地形无关方案利用严格成像模型生成；需要利用n个真实地面控制点对有理函数模型参数进行优化；其特征为：包括以下步骤，
步骤2.1，对n个真实控制点的像平面坐标和地面坐标进行标准化；
步骤2.2，通过变形将有理函数模型(1)转化为关于各RPCs的线性模型(2)，然后将n个控制点所列的误差方程以矩阵形式来表示，如式(3)所示，该步骤与1.2相同；
步骤2.3，将向量x分解为两个分量的和
x＝x⁰+Δx            (5)
其中
x0=[a00,...,a190,b10,...,b190,c00,...,c190,d10,...,d190]T]]>为有理函数模型参数的初始值，从卫星影像自带的RPC文件中获取或通过地形无关方案利用严格成像模型生成，
Δx=[a0δ,...,a19δ,b1δ,...,b19δ,c0δ,...,c19δ,d1δ,...,d19δ]T]]>为有理函数模型参数的改正量。
然后模型(3)可转换为
L′＝AΔx+V              (6)
其中L′＝L-Ax⁰，V为随机误差。
步骤2.4，对模型(6)加入l₁范数约束，将问题转化为带l₁范数约束的最优化问题：
min||AΔx-L′||22s.t.||Δx||1≤α---(7)]]>
其中α为l₁范数约束因子；
步骤2.5，利用LARS算法求解得到改正向量Δx的估计值
步骤2.6，根据估计值修正有理函数模型参数，得到有理函数模型参数的估计值
x^=x0+Δx^---(8)]]>
然后生成改正后的有理函数模型。

基于1₁范数约束的有理函数模型参数求解和全参数优化方法
技术领域
本发明涉及有理函数模型参数的求解和优化方法，能够自适应地精简有理函数模型，减少所需控制点的数量，提高图像校正精度。可应用于遥感、摄影测量、测绘、图像处理等领域。 
背景技术
有理函数模型(Rational Function Model，RFM)是一种通用的几何成像模型，它具有良好的内插特性和连续性并独立于传感器和平台，此外还具有拟合精度高、成像参数及模型保密、计算速度快等优势，因而被广泛的应用于卫星遥感影像的几何校正中。 
有理函数模型将像点坐标表示为以相应地面点空间坐标为自变量的有理多项式的比值，其标准方程为如式(1)所示： 
l=Nl(X,Y,Z)Ds(X,Y,Z)s=Nl(X,Y,Z)Ds(X,Y,Z)---(1)]]>
其中l和s是像点在像平面上的行列值标准化后的结果[OGC(OpenGIS Consortium)，1999.The OpenGIS^TM Abstract Specification-Topic7：The Earth Imagery Case.]，(X，Y，Z)是物方点空间坐标标准化后的结果， 
N_l(X，Y，Z)＝a₀+a₁Z+a₂Y+a₃X+a₄ZY+a₅ZX+a₆YX+a₇Z²
+a₈Y²+a₉X²+a₁₀ZYX+a₁₁Z²Y+a₁₂Z²X+a₁₃ZY²
+a₁₄Y²X+a₁₅ZX²+a₁₆YX²+a₁₇Z³+a₁₈Y³+a₁₉X³
D_l(X，Y，Z)＝b₀+b₁Z+b₂Y+b₃X+b₄ZY+b₅ZY+b₆YX+b₇Z²
+b₈Y²+b₉X²+b₁₀ZYX+b₁₁Z²Y+b₁₂Z²X+b₁₃ZY²
+b₁₄Y²X+b₁₅ZX²+b₁₆YX²+b₁₇Z³+b₁₈Y³+b₁₉X³
N_s(X，Y，Z)＝c₀+c₁Z+c₂Y+c₃X+c₄ZY+c₅ZX+c₆YX+c₇Z²
+c₈Y²+c₉X²+c₁₀ZYX+c₁₁Z²Y+c₁₂Z²X+c₁₃ZY²
+c₁₄Y²X+c₁₅ZX²+c₁₆YX²+c₁₇Z³+c₁₈Y³+c₁₉X³
D_s(X，Y，Z)＝d₀+d₁Z+d₂Y+d₃X+d₄ZY+d₅ZX+d₆YX+d₇Z²
+d₈Y²+d₉X²+d₁₀ZYX+d₁₁Z²Y+d₁₂Z²X+d₁₃ZY²
+d₁₄Y²X+d₁₅ZX²+d₁₆YX²+d₁₇Z³+d₁₈Y³+d₁₉X³
式中，a_i，b_i，c_i，d_i(i＝0，1，…，19)即为有理多项式系数(Rational Polynomial Coefficients，RPCs)，其中b₀和d₀的值为1。 
对影像坐标行列值的标准化公式为： 
l=Line-LINE_OFFLINE_SCALEs=Sample-SAMP_OFFSAMP_SCALE---(2)]]>
其中LINE_OFF、LINE_SCALE、SAMP_OFF和SAMP_SCALE为影像坐标的标准化因子。Line表示影像行坐标，Sample表示影像的列坐标。 
对地面点的标准化公式为： 
X=Longitude-LONG_OFFLONG_SCALE]]>
Y=Latitude-LAT_OFFLAT_SCALE---(3)]]>
Z=Height-HEIGHT_OFFHEIGHT_SCALE]]>
其中LONG_OFF、LONG_SCALE、LAT_OFF、LAT_SCALE、HEIGHT_OFF、HEIGHT_SCALE为地面坐标标准化因子。Longitude表示地面点的经度、Latitude表示地面点的纬度、Height表示地面点的高程。 
有理函数模型的使用主要包括两大类：地形无关的方案和地形相关的方案： 
1)地形无关的方案是通过卫星遥感影像的严格成像模型生成多个高程面的虚拟控制点格网，再利用这些虚拟控制点计算地形无关的有理函数模型参数(RPCs)。最后利用真实地面控制点对地形无关的有理函数模型进行像方或物方补偿，实现影像的精校正。例如像方多项式模型进行误差补偿的具体形式如(4)所示： 
Δl＝a₀+a_l·l′+a_s·s′+a_ls·l′·s′ 
+a_l2·l^′2+a_s2·s^′2+… 
                        (4) 
Δs＝b₀+b_l·l′+b_c·s′+b_ls·l′·s′ 
+b_l2·l^′2+b_s2·s^′2+… 
其中，Δl和Δs是像平面坐标计算值与实际值之间的偏差，l′和s′是平面坐标的计算值，a₀，a_l，a_s，…和b₀，b_l，b_s，…是多项式模型的系数。较为常用的是像方仿射变换模型，最少需要3个控制点就能取得比较理想的改正效果。 
2)地形相关方案则是直接利用大量的真实控制点解算有理函数模型参数，得到地形相关的有理函数模型，然后利用地形相关的有理函数模型进行影像的精校正。 
其中地形无关的方案应用较为广泛，而地形相关方案的应用相对较少，主要原因是需要大量的控制点且模型的解算很不稳定。 
近年来发射的和将来要发射的高分辨率卫星(如ALOS、GeoEye、Pleiades、SPOT-6/7、ZY-3、GF-1/2等)均提供地形无关的RPCs来描述其影像的成像模型，用户可以十分方便地利用这些参数对影像进行几何处理，在精确控制点和DEM的支持下，还能进行高精度的正射校正。 
尽管有理函数模型具有很多的优势，但同时还存在一些缺点：1)模型包含78个RPCs，至少需要39个控制点才能求解；2)78个RPCs之间存在很强的相关性，这使得RPCs的求解成为一个病态性问题；3)有理函数模型是纯数学模型，缺乏几何意义；4)有理函数模型在窄视场角的高分辨率影像中应用很成功，但是对于宽视场角和定向参数不稳定的情况下适用性较差。 
法方程的病态性是求解RPCs的主要困难，常用通过正则化(例如Tikhonov正则化、岭估计、Levenberg-Marquardt方法等)或奇异值截断的方法将法方程修正为良态，然后进行求解。这些方法可以在一定程度上改善方程的病态性、防止“过拟合”现象，但是仍然需要大量的控制点来解算RPCs。当控制点的数量不足够多时，得到的解往往很不稳定。此外，对于HJ-1A/B之类的国产卫星，由于传感器视场角大、平台姿态不稳定，图像的内部畸变比较大，在解算RPCs时更容易出现模型不稳定的现象，且影像的精确纠正很难通过简单的像方改正来实现，这是的有理函数模型无法应用于这类卫星遥感影像。 
鉴于有理函数模型的通用性、简单性、保密性、计算量小等优势，克服现有方法中存在的问题，从而提高有理函数模型计算的稳定性及其适用范围具有重要的实用价值。 
发明内容
本发明的目的在于解决现有技术的不足，提出了一种基于l₁范数约束的有理函数模型求解方法和一种对有理函数模型全参数进行优化的精纠正方法。在地形无关方案中，当卫星影像内部畸变较大时，仍可以得到可靠的有理函数模型参数；在地形相关方案中，可利用少于39个控制点求解78个有理函数模型系数；在已知有理函数模型系数的条件下，利用少量地面控制点对有理函数模型进行全参数优化。本发明提高了有理函数模型参数求解的稳定性，减少了所需控制点的个数，可以在少于39个控制点的条件下求出有理函数模型参数的可靠解，并提高了影像精校正的精度。 
本发明的技术方案如下： 
一种基于l₁范数约束的有理函数模型求解方法，所述方法利用控制点进行求解，控制点可以是真实地面控制点，也可以是通过严格成像模型生成的虚拟控制点；包括以下步骤， 
步骤1.1，对所有控制点的像平面坐标和地面坐标进行标准化； 
步骤1.2，通过简单的变形将有理函数模型转化为关于各RPCs的线性模型，并将控制点所列的误差方程用矩阵形式来表示； 
步骤1.3，对线性模型模型加入l₁范数约束，将问题转化为带l₁范数约束的最优化问题； 
步骤1.4，利用LARS(Least Angle Regression)算法[Efron，Bradley，Trevor Hastie，lain Johnstone，and Robert Tibshirani.″Least angle regression.″The Annals of statistics32，no.2(2004)：407499.]求解带l₁范数约束的最优化问题； 
步骤1.5，得到有理函数模型参数RPCs，生成有理函数模型。 
一种对有理函数模型全参数进行优化的精纠正方法，所述方法中已有的有理函数模型参数从卫星影像自带的RPC文件中获取或通过地形无关方案利用严格成像模型生成；需要利用真实地面控制点对有理函数模型参数进行优化；包括以下步骤， 
步骤2.1，对所有真实地面控制点的像平面坐标和地面坐标进行标准化； 
步骤2.2，通过简单的变形将有理函数模型转化为关于各RPCs的线性模型，然后控制点所列的误差方程可以用矩阵形式来表示； 
步骤2.3，将有理函数模型参数的初始值带入误差方程，得到关于有理函数模型参数改正量的误差方程； 
步骤2.4，对模型加入l₁范数约束，将问题转化为带l₁范数约束的最优化问题； 
步骤2.5，利用LARS算法求解有理函数模型参数改正量的估计值； 
步骤2.5，修正有理函数模型参数，生成改正后有理函数模型。 
附图说明
图1是本发明实施例一流程图； 
图2是本发明实施例二流程图； 
图3是本发明实施例三流程图。 
具体实施方式
本发明提出两种方法：一种是基于l₁范数约束的有理函数模型求解方法，另一种是对有理函数模型全参数进行优化的精纠正方法。第一种方法可以分别应用于地形无关方案的有理函数参数求解和地形相关方案的有理函数模型参数求解，第二种方法可以用于已知初始有理函数模型参数前提下的有理函数模型精化。以下结合三个实施例和附图1～3详细说明本发明 技术方案。 
(1)实施例一提供基于地形无关方案的有理函数模型参数求解方法，参见附图1： 
步骤1.1，根据卫星遥感影像的严格成像模型生成虚拟控制点格网； 
卫星遥感影像的严格成像模型是根据严格的物理成像过程建立的，用来描述像素坐标(行列号)与大地经纬度坐标(经度、纬度和高程)之间的映射关系，根据像素坐标和大地经纬度坐标的映射方向可分为正算模型(式(5))和反算模型(式(6))。 
(X，Y)＝F(l，s，Z)              (5) 
(l，s)＝F^-1(X，Y，Z)             (6) 
式(5)和式(6)中l表示像平面点的行坐标，s表示像平面点的列坐标，X表示经度，Y表示纬度，Z表示高程，F(·)表示正算模型，F^-1(·)表示反算模型。 
首先利用严格成像模型计算影像的覆盖范围，再根据全球DEM数据和影像覆盖范围统计影像覆盖范围内的最大和最小高程值，并在最大和最小高程之间建立7个高程面，在每个高程面按一定的间隔利用严格成像模型的正算模型(X，Y)＝F(l，s，Z)生成21*21的均匀分布的控制点格网，这样就建立了7个不同高程面的三维控制点格网(n个控制点)。 
步骤1.2，按式(2)和式(3)对n个控制点的像平面坐标和地面坐标进行标准化； 
步骤1.3，通过变形将有理函数模型(1)转化为关于各RPCs的线性模型，如式(7)所示 
Nl(X,Y,Z)-lDl(X,Y,Z)=0Ns(X,Y,Z)-sDs(X,Y,Z)=0---(7)]]>
然后n个控制点所列的误差方程可以矩阵形式来表示，如式(8)所示， 
L＝Ax+V               (8) 
其中 
L＝[l₁，l₂，…，l_n，s₁，s₂，…，s_n]^T， 
x＝[a₀，…，a₁₉，b₁，…b₁₉，c₀，…，c₁₉，d₁，…，d₁₉]^T， 
V为随机误差， 
A＝[A₁，A₂，…，A_n]^T， 
Ai=[1,Xi,Yi,...,Zi3,-liXi,-liYi,...,-liZi3,1,Xi,Yi,...,Zi3,-siXi,-siYi,...,-siZi3]T,]]>
(l_l，s_l)为第i个控制点标准化后的像平面行列值， 
(X_i，Y_i，Z_i)为第i个控制点标准化后的地面坐标(经度、纬度和高程)， 
i＝1，2，…，n； 
步骤1.4，对模型(8)加入l₁范数约束，将问题转化为带l₁范数约束的最优化问题： 
min||Ax-L||22s.t.||x||1≤α---(9)]]>
其中 
α为l₁范数约束因子，通常可取值为10^-5； 
||·||₁和||·||₂分别表示l₁范数和l₂范数，例如对于任意向量x： 
||x||₁＝∑_i|x_l|表示向量x的l₁范数， 
表示向量x的l₂范数； 
步骤1.5，利用LARS(Least Angle Regression)算法[Efron，Bradley，Trevor Hastie，lain Johnstone，and Robert Tibshirani.″Least angle regression.″The Annals of statistics32，no.2(2004)：407-499.]求解得到未知向量x的估计值
步骤1.6，根据估计值得到有理函数模型参数RPCs，生成有理函数模型。 
(2)实施例二提供基于地形相关方案的有理函数模型参数求解方法，参见附图2： 
步骤2.1，按式(2)和式(3)对n个真实控制点的像平面坐标和地面坐标进行标准化； 
步骤2.2，通过变形将有理函数模型(1)转化为关于各RPCs的线性模型，如式(7)所示，然后将n个控制点所列的误差方程以矩阵形式来表示，如式(8)所示，该步骤与13相同； 
步骤2.3，对模型(8)加入l₁范数约束，将问题转化为带l₁范数约束的最优化问题； 
步骤2.4，利用LARS算法求解得到未知向量x的估计值
步骤2.5，根据估计值得到有理函数模型参数RPCs，生成有理函数模型。 
(3)实施例三提供已知初始有理函数模型参数条件下的有理函数模型精化，参见附图3： 
步骤3.1，按式(2)和式(3)对n个真实控制点的像平面坐标和地面坐标进行标准化； 
步骤3.2，通过简单的变形将有理函数模型(1)转化为关于各RPCs的线性模型，如式(7)所示，然后n个控制点所列的误差方程可以矩阵形式来表示，如式(8)所示，该步骤与1.3相同； 
步骤3.3，将向量x分解为两个分量的和 
x＝x⁰+Δx                (10) 
其中 
x0=[a00,...,a190,b10,...,b190,c00,...,c190,d10,...,d190]T]]>为有理函数模型参数的初始值，从卫星影像自带的RPC文件中获取或通过地形无关方案利用严格成像模型生成， 
Δx=[a0δ,...,a19δ,b1δ,...,b19δ,c0δ,...,c19δ,d1δ,...,d19δ]T]]>为有理函数模型参数的改正量。 
然后模型(8)可转换为 
L′＝AΔx+V       (11) 
其中L′＝L-Ax⁰，V为随机误差。 
步骤3.4，对模型(11)加入l₁范数约束，将问题转化为带l₁范数约束的最优化问题： 
min||AΔx-L′||22s.t.||Δx||1≤α---(12)]]>
其中α为l₁范数约束因子，通常可取值为10^-5； 
步骤3.5，利用LARS(Least Angle Regression)算法求解得到改正向量Δx的估计值
步骤3.6，根据估计值修正有理函数模型参数，得到有理函数模型参数的估计值
x^=x0+Δx^---(13)]]>
然后生成改正后的有理函数模型。