CN102855221A - 一种基于嵌套回归的有理函数模型参数自动优选及求解方法 - Google Patents

Abstract

本发明将有理多项式系数(Rational Polynomial Coefficients，RPCs)求解问题转化为具有复共线性的多元线性回归问题，提出一种以可决系数为评价准则的嵌套回归方法，自动优选有理函数参数。通过逐步引入有理函数模型中影响显著的系数，舍去影响不显著的冗余系数，实现RPCs的自动优选及解算，得到精简的有理函数模型。利用本发明方法能够自动简化有理函数模型，降低控制点使用数量，提高图像校正和配准精度。

Description

一种基于嵌套回归的有理函数模型参数自动优选及求解方法

技术领域

本发明涉及一种有理函数模型参数的自动优选及求解方法，能够简化有理函数模型，降低控制点使用数量，提高图像校正精度。可应用于遥感、摄影测量、测绘、图像处理等领域。

背景技术

有理函数模型(Rational Function Model，RFM)是一种完全的数学模型，其实质就是利用大量的控制点，来拟合严格成像关系，不要求了解传感器的实际特性和成像过程，它适用于不同类型的传感器，并且可以隐藏传感器在成像过程中的各种物理参数，因此在摄影测量和遥感领域得到了广泛的应用。利用有理函数模型处理卫星遥感影像的关键在于精确求解其参数，即有理多项式系数(Rational Polynomial Coefficients，RPCs)。由于有理函数系数过多，系数之间产生强相关，结果很不稳定，甚至不能得到可靠的解[LIN Xianyong，YUAN Xiuxiao，Improvement of the stability solving rational polynomial coefficients.In the Intemational Archivesof Photogrammetry，Remote Sensing and Spatial Information Sciences，ISPRS，Beijing，2008，vol.XXXVII，part B2.]。

十多年来，国内外学者对RFM做了广泛的研究。开放地理信息协会(OpenGIS Consortium，OGC)将RFM中的坐标值标准化为-1到+1之间的数值[OGC(OpenGIS Consortium)，1999.TheOpenGISTM Abstract Specification-Topic 7：The Earth Imagery Case.]，减小了法方程系数矩阵的条件数。Tao等人深入地研究了RFM的理论，提出用Tikhonov正则化和L曲线法来克服RFM的病态性[TAO C V，HU Y.Investigation of the rational function model.Proceedings of ASPRSAnnual Convention，Washington D.C.，May 22-26，2000.]，袁修孝等人比较了岭迹法、L曲线法、经验公式法以及广义岭估计等岭估计的方法用于RPCs求解的效果，证实了L曲线法的有效性[袁修孝，林先勇.基于岭估计的有理多项式参数求解方法.武汉大学学报(信息科学版).2008，33(11)：1130-1133.)。也可以通过Levenberg-Marquardt方法、奇异值分解的方法减小法方程系数矩阵的条件数，此外，将RFM中的分母项设置为常数也能在一定程度上降低模型的病态性[FRASERA C S，DIALB G，GRODECKI J.Sensor orientation via RPCs.ISPRS Journal ofPhotogrammetry And Remote Sensing，2006，60(2006)：182-194.]。目前这些方法中岭估计(尤其是L曲线法)的应用最为广泛，对于克服病态性也十分有效，但是RPCs的解算仍然存在一些问题，如岭估计属于有偏估计、解算模型需要大量的控制点、模型缺乏物理意义等。

鉴于到岭估计方法的不足，可以考虑通过优选RPCs的方法简化有理函数模型，从而克服模型的病态性。袁修孝等利用复共线性分析的方法对RPCs进行了优选，提高了RPCs求解的稳定性[袁修孝，曹金山.一种基于复共线性分析的RPC参数优选法.武汉大学学报(信息科学版)，2011，36(6)：665-669.]，这种方法是以复共线性为评价准则的剔除优选法。然而由于变量间多重相关性的形式非常复杂，而且还缺乏十分可靠的检验方法，删除部分相关变量的做法常导致增大模型的解释误差，将本应保留的信息舍弃[王惠文.偏最小二乘回归方法及其应用.北京：国防工业出版社，1999.]。

发明内容

本发明的目的是克服现有方法的不足，将RPCs求解问题转化为具有复共线性的多元线性回归问题，提出一种以可决系数为评价准则的嵌套回归方法，自动优选有理函数参数。通过逐步引入有理函数模型中影响显著的系数，舍去影响不显著的冗余系数，得到精简的有理函数模型，进而求得稳定可靠的RPCs。利用本发明方法能够自动简化有理函数模型，降低控制点使用数量，提高图像校正精度。

为了实现上述目的，本发明的技术路线如下：

第一步、将RPCs求解问题转化为多元线型回归问题

(1-1)消去标准RFM中的分母项并移项，得到关于RPCs的线性模型

(1-2)利用线性化的RFM对每个控制点建立观测方程

(1-3)联立所有观测方程，得到初始多元线性回归方程

第二步、用嵌套回归方法优选并求解RPCs

(2-1)设置待选RPC参数集为标准RFM中的所有RPC参数，已选RPC参数集为空，目标向量为初始多元线性回归方程的目标向量，设置迭代停止条件

(2-2)分别计算待选RPC参数集中各参数对目标向量的可决系数

(2-3)选出(2-2)中最大可决系数对应的RPC参数，将其从待选RPC参数集移至已选RPC参数集，并记录相应的回归系数

(2-4)如果满足停止条件，则转至(2-5)；否则，更新目标向量，并转至(2-2)

(2-5)将已选RPC参数集中的RPCs设置为其对应的回归系数，将待选RPC参数集中的RPCs设为0，合并已选RPC参数集和待选RPC参数集，得到RPCs的解

附图说明

图1是本发明的主流程示意图

图2是控制点充足条件下均方根误差和法方程系数矩阵条件数随RPCs个数的变化曲线，其中图2(a)是行RPCs的情况，图2(b)是列RPCs的情况

图3是10个控制点时均方根误差和法方程系数矩阵条件数随RPCs个数的变化曲线，其中图3(a)是行RPCs的情况，图3(b)是列RPCs的情况

具体实施方式

本发明方法的主要流程为：首先将RPCs求解问题转化为多元线性回归问题；然后采用嵌套回归的方法，逐步引入有理函数模型中影响显著的系数，自动优选并求解出RPCs。

1)建立RPCs的多元线性回归方程

有理函数模型将像点坐标表示为以相应地面点空间坐标为自变量的有理多项式的比值，其标准方程为如式(1)所示：

$\left\{\begin{array}{c}r=\frac{N_r(X,Y,Z)}{D_r(X,Y,Z)}\\ c=\frac{N_c(X,Y,Z)}{D_c(X,Y,Z)}\end{array}\right.---\left(1\right)$

其中r和c是像点在像平面上的行列值标准化后的结果[OGC(OpenGIS Consortium)，1999.The OpenGIS^TM Abstract Specification-Topic 7：The Earth Imagery Case.]，(X，Y，Z)是物方点空间坐标标准化后的结果，

N_r(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_r(X，Y，Z)＝b₀+b₁Z+b₂Y+b₃X+b₄ZY+b₅ZX+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_c(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_c(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)即为有理多项式系数(RPCs)，其中b₀和d₀的值为1。

通过简单的变形将有理函数模型转化为关于各RPCs的线性模型，如式(2)所示：

$\left\{\begin{array}{c}{N}_r(X,Y,Z)-r{D}_r(X,Y,Z)=0\\ N_c(X,Y,Z)-c{D}_c(X,Y,Z)=0\end{array}\right.---\left(2\right)$

由式(2)可知行、列RPCs是相互独立的，因此可分开单独求解。以行RPCs为例，(2)式中的第一式可写为：

r＝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³-b₁rZ-b₂rY-b₃rX-b₄rZY(3)-b₅rZX-b₆rYX-b₇rZ²-b₈rY²-b₉rX²-b₁₀rZYX-b₁₁rZ²Y-b₁₂rZ²X-b₁₃rZY²-b₁₄rY²X-b₁₅rZX²-b₁₆rYX²-b₁₇rZ³-b₁₈rY³-b₁₉rX³

联立所有观测方程则有：

r＝Xβ_r+ε       (4)

其中β_r＝[a₀，a₁，…，a₁₉，b₁，b₂，…，b₁₉]^T， $r=\left[\begin{array}{c}{r}_1\\ r_2\\ .\\ .\\ .\\ r_n\end{array}\right],$ n为观测方程的个数，ε为随机误差，X＝[1，X₁，X₂，…，X₃₈]，1表示各行均为1的列向量，X_i(i＝1，2，…，n)为矩阵X的各列向量，如X₁＝[Z₁，Z₂，…，Z_n]^T，X₃₈＝[-r₁X₁ ³，-r₂X₂ ³，…，-r_nX_n ³]^T等。

2)用嵌套回归方法优选并求解RPCs

本发明以可决系数为优选RPCs的评价准则。在统计学中，所谓可决系数就是回归平方和在总变差中所占的比重[CAMERON A C，WINDMEIJER F A G.An R-squared Measure ofGoodness of Fit for Some Common Nonlinear Regression Models.Journal of Econometrics，1997，77(2)：329-342.]，通过式(5)来计算：

$R^2=\frac{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{({\hat{y}}_i-\stackrel{\‾}{y})}^2}{\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{(y_i-\stackrel{\‾}{y})}^2}---\left(5\right)$

其中n为观测值的个数，y_i为各观测值，

为各观测值的均值，

为各观测量的回归结果。R²的值在0～之间，R²越接近1，说明在总变差中由模型作出了解释的部分占的比重越大，模型拟合程度越好；反之R²越接近0，说明模型对样本观测值的拟合程度越差[庞皓.计量经济学.成都：西南财经大学出版社，2001：132-141.]。

(4)式所示的多元线性回归问题的完整形式如下：

r＝Xβr₊ε          (6)

其中r、X和β_r同(4)式中各项，ε为随机误差，记X＝[1，X₁，X₂，…，X₃₈]，1表示各行均为1的列向量，X_i为矩阵X的各列向量，如：X₃₈＝[-r₁X₁ ³，-r₂X₂ ³，…，-r_nX_n ³]^T，X₁＝[Z₁，Z₂，…，Z_n]^T等。

采用常规回归分析方法得到的结果的完整形式如下：

$\hat{r}=X{\widehat{\mathrm{\β}}}_r---\left(7\right)$

其中

为β_r的估计量，

为拟合结果。

基于嵌套回归的RPCs自动优选及求解算法如下：

步骤1：设置收敛阈值t₁＝0.5/scale，t₂＝0.05/scale，其中scale为坐标标准化过程中的缩放倍数。令k＝2，r₁＝r，no为控制点个数，nv为行(列)RPCs最大个数，转入步骤2；

步骤2：依次用X_i和r_k建立线性回归模型

r_k＝β_k01+β_k1X_i+ε_k       (8)

分别计算各模型的可决系数，假设在所有X_i中X_(m)对应的模型的可决系数最大，利用常规最小二乘法对该模型进行回归分析可得到

${\hat{r}}_k={\widehat{\mathrm{\β}}}_{k0}1+{\widehat{\mathrm{\β}}}_{k1}{X}_{\left(k\right)}---\left(9\right)$

转入步骤3；

步骤3：若k≥no或k≥nv，转入步骤4；否则，计算残差向量

和均方根误差如果σ_k＜t₁且|σ_k-σ_k-1|＜t₂，则转到步骤4；否则，令k＝k+1，转入步骤2；

步骤4：得到优选变量集{X_(m)|m＝1，2，…，k}，对线性回归模型(10)进行常规最小二乘回归，结果为式(11)。

r＝β₀1+β₀X₍₁₎+β₂X₍₂₎+…+β_kX_(k)+ε     (10)

$\hat{r}={\widehat{\mathrm{\β}}}_{0}1+{\widehat{\mathrm{\β}}}_1{X}_{\left(1\right)}+{\widehat{\mathrm{\β}}}_2{X}_{\left(2\right)}+...+{\widehat{\mathrm{\β}}}_k{X}_{\left(k\right)}---\left(11\right)$

将β_r中X_(m)所对应的各系数设为

其余系数设为0即得到待求的RPCs，算法结束。

该算法的流程见附图1。

这里优选变量的停止条件为均方根小于阈值t₁且最后两次迭代的均方根差小于阈值t₂，通常可将t₁设为0.5个像素，t₂设为0.05个像素。这种停止条件的物理意义非常明确：如果控制点的拟合误差达到0.5个像素，且加入一个变量后控制点的拟合精度没有明显提高，则停止加入新的变量。

应用实例

下面以两个实例来说明本发明方法的具体应用。

应用实例一：

选择西安地区的一景QuickBird P2AS影像进行试验，该影像空间分辨率为0.61米，地面最大高程差约为650米，影像范围内均匀分布80个地面控制点。利用本文方法和80个地面控制点，逐步从39个行(列)RPCs中选出重要的系数，记录各个步骤中模型拟合的均方根误差和法方程系数矩阵的条件数，可以得到均方根误差和法方程系数矩阵条件数随RPCs个数的变化曲线，如图2所示。图2(a)为行RPCs的情况，图2(b)为列RPCs的情况，图2(a)和图2(b)均由均方根误差变化曲线和法方程系数矩阵条件数变化曲线组成，横坐标为RPCs的个数，纵坐标分别为均方根误差和条件数的对数值。这里使用对数的目的是缩小各数值间的数量级差，从而提高数据的区分度。

从图2可以看出，当RPCs个数较少时，均方根误差较大、法方程系数矩阵条件数较小，随着RPCs个数的增加，均方根误差逐步减小，法方程系数矩阵条件数逐步增大。但是当RPCs的个数达到10个左右时，均方根误差减小的幅度变得很小，继续不断加入新的系数时均方根误差甚至会增大。不难分析出原因：当RPCs个数较少时，系数间的相关性小，因此法方程系数矩阵条件数较小；但由于模型简单，不能较好地拟合观测数据，因此均方根误差较大。随着RPCs个数的增加，拟合精度逐步提高，但系数间的相关性也在变大，当系数个数较多时，法方程系数矩阵条件数急剧增大，模型出现明显的病态性，非常不稳定。

从80个地面控制点中均匀的选出10个控制点，然后利用本发明方法，逐步从39个行(列)RPCs中选出重要的系数，记录各个步骤中模型拟合的均方根误差和法方程系数矩阵的条件数，可以得到均方根误差和法方程系数矩阵条件数随RPCs个数的变化曲线，如图3所示，图3的结构与图2相同。

由图3可知，当RPCs个数小于控制点的个数时，均方根误差随着系数个数的增加逐渐变小；系数个数达到控制点的个数后，均方根误差急剧减小，法方程系数矩阵条件数急剧增大，模型变得很不稳定。

实例一的结果说明两个问题：

1)在控制点充足的条件下，当优选了一定数量的RPCs后，继续增加RPCs并不能明显改善拟合精度，甚至会降低拟合精度，而模型的病态性会显著提高。

2)当控制点较少(10个)时，多于控制点个数的RPCs无法求解。

应用实例二：

使用Quickbird P2AS、ALOS PRISM 1B2、SPOT5HRG 1A和Landsat5L2影像各一景进行RPCs求解试验，各影像均有分布良好的80个控制点，分别采用普通最小二乘法、岭估计法(L曲线)及基于嵌套回归的RPCs自动优选法求解各影像的RPCs，其中嵌套回归方法所选定的阈值分别为t₁＝0.5像素t₂＝0.05像素。试验结果如表1所示。

表13种方法解算RPCs时法方程系数矩阵的条件数

表1列出了不同卫星影像数据采用3种方法解算RPCs时法方程的条件数，以及嵌套回归方法所优选的行(列)RPCs的个数。分析表1中的数据，可以看出：

1)直接采用最小二乘方法求解RPCs，法方程系数矩阵的条件数高达10¹⁸，病态性非常严重；使用岭估计方法后，法方程系数矩阵的条件数降到10⁵～10¹⁰级别，病态性得到减弱，但是系数间的相关性仍然很强；采用本文的嵌套回归方法，从39个行(列)RPCs中优选出少量的系数，法方程的条件数下降到10¹～10³级别，病态性基本消除。

2)对于不同的影像，用本文方法所优选出的RPCs个数也不同，这说明不同影像的变形程度及成像几何复杂度是不同的。其中，Landsat5L2的模型最简单，所须的系数最少，SPOT5HRG 1A的模型最复杂，所须的系数最多。一般来说图像产品的处理级别越高，图像内部的变形程度越小，其几何校正模型就越简单。

3)由表中数据可知，同一影像所需的行RPCs一般多于所需的列RPCs，这说明影像在行方向和列方向的变形程度并不相同，在行方向的变形更复杂，且由下文的试验结果(表2)可知，通常情况下，图像在列方向的纠正精度也高于行方向。

4)嵌套回归方法优选出的RPCs越少，最小二乘解法中法方程系数矩阵的条件数就越高。选出的RPCs越少，说明该影像的几何校正模型越简单，78个RPCs之间的相关性就越大，因此模型的病态性就越明显。

从Quickbird、ALOS、SPOT5和Landsat5影像的60个控制点中分别均匀的选出20个和10个控制点，分别使用最小二乘法、L曲线法和嵌套回归方法求解RPCs，用控制点评价模型的拟合精度，用60个检查点来评价几何处理精度，结果如表2所示。在不同数量控制点的条件下，各影像所自动优选出的行(列)RPCs个数与表1中的个数基本一致，差别不超过2个系数。

表2RPCs解算精度

分析表2中的数据可以看出：

1)优选部分RPCs可以减小参数冗余并增加模型的稳定性。在控制点充足(60个)的条件下，最小二乘法和L曲线法均能得到较高的拟合精度，但是由图1的分析可知，使用39个行(列)RPCs拟合控制点数据很不稳定，因此，部分情况(ALOS和SPOT5)下，检查点的精度会很差。使用L曲线法可以增加模型的稳定性，检查点的精度得到了保证。使用嵌套回归的方法可以得到稳定的解，由于只选取了部分RPCs，控制点的拟合精度并没有最小二乘法和L曲线法高，但是检查点的精度不差于L曲线法的结果。当然也可以通过调整阈值来选取更多的RPCs，从而得到更高的拟合精度，但这样并不一定能提高检查点的精度。

2)嵌套回归的方法可以解决在控制点数不足39个时RPCs的求解问题，一般情况下，不到20个控制点就可以得到稳定可靠的解，建立精确的有理函数模型。。当控制点数量不足39个时，最小二乘法的拟合误差为0，但是检查点的误差非常大，通过L曲线法引入岭参数可以在一定程度上减小检查点的误差，但误差仍然在10³个像素以上。这就是多元回归分析中所谓的过度拟合(over-fitting)现象，此时，最小二乘法和岭估计方法都无法得到正确的RPCs。使用嵌套回归的方法优选RPCs，当所选系数个数少于控制点的个数时，可以得到稳定可靠的RPCs解，检查点的精度也得到了保证。

3)当控制点个数少于模型所需的最少RPCs个数时，本文方法不能得到稳定地解。利用嵌套回归方法和10个控制点解算SPOT5影像的RPCs时，检查点在行方向的误差达到20个像素(见表2)，这是因为该影像所需要的行RPCs为14个(见表1)，仅用10个控制点不能得到稳定的解。

4)使用嵌套回归的方法求解RPCs时，列方向的几何精度一般高于行方向。

5)不同的遥感影像所需的控制点个数也不一样。对于Quickbird、ALOS、SPOT5数据，当控制点数减少到10个时，检查点的精度会明显下降，但是对于Landsat5数据，10个控制点和60个控制点的校正精度基本一致，这说明Landsat5影像的成像模型比较简单，10个控制点就足够了，加入更多的控制点并不会显著地提高校正的精度。

Claims (3)

1.一种基于嵌套回归的有理函数模型参数自动优选及求解方法，其步骤为：

第一步、将RPCs求解问题转化为多元线型回归问题

第二步、用嵌套回归方法优选并求解RPCs。

2.如权利要求1所述的方法，其中，第一步包括以下三个步骤：

(2-1)消去标准RFM中的分母项并移项，得到关于RPCs的线性模型

(2-2)利用线性化的RFM对每个控制点建立观测方程

(2-3)联立所有观测方程，得到初始多元线性回归方程。

3.如权利要求2所述的方法，其中，第二步包括以下五个步骤：

(3-1)设置待选RPC参数集为标准RFM中的所有RPC参数，已选RPC参数集为空，目标向量为初始多元线性回归方程的目标向量，设置迭代停止条件

(3-2)分别计算待选RPC参数集中各参数对目标向量的可决系数

(3-3)选出(3-2)中最大可决系数对应的RPC参数，将其从待选RPC参数集移至已选RPC参数集，并记录相应的回归系数

(3-4)如果满足停止条件，则转至(3-5)；否则，更新目标向量，并转至(3-2)

(3-5)将已选RPC参数集中的RPCs设置为对应的回归系数，将待选RPC参数集中的RPCs设为0，合并已选RPC参数集和待选RPC参数集，得到RPCs的解 。

Images