CN100481133C - 基于有理多项式成像模型的正变换方法 - Google Patents

Abstract

本发明涉及一种RPC模型正变换方法，该方法将反变换结果引入正变换过程进行迭代变换，其过程包括以下步骤：设置高程、和RPC模型参数，设定地面矩形区域；计算地面矩形区域4角点和中心点的影像像素坐标；计算仿射变换系数；计算量测像点坐标对应的近似经纬度坐标；计算理论像点坐标；计算残差e，如果残差e满足设定的差值范围，输出近似经纬度坐标作为量测像点坐标对应的经纬度坐标，否则调整地面矩形区域，然后返回迭代，直到残差e满足差值范围。本发明的技术方案能够确保在影像供应商提供的有限条件下获得精确的地面点位置，在遥感影像应用领域具有重要意义。

Description

基于有理多项式成像模型的正变换方法

技术领域

本发明属于摄影成像技术领域，特别涉及基于有理多项式成像模型的正变换方法。

背景技术

由于严密成像几何模型具有较高的定位精度，因此一直是摄影测量学的首选，卫星遥感影像处理常常要利用严密成像几何模型。然而，一些高性能的传感器***虽然实现了商业化，但是由于卫星平台和有效载荷等关键技术参数是航天大国的技术机密，高分辨率传感器的核心信息和卫星轨道参数并未公开，如IKONOS等卫星，无法利用严格成像几何模型进行处理；同时传感器成像方式的多样化也对摄影测量软件提出了模型通用化的要求，这样在软件框架设计时就不必为现有的或将来可能会出现的各种类型的传感器一一确定其严格几何模型的形式，而可以统一采用通用模型进行处理，从而大大降低了程序设计的复杂性，更易于软件升级和维护，特别在同时处理多源卫星遥感影像数据时，通用成像几何模型更能显示出其优势。

常见的通用成像几何模型有以下几种：多项式、直接线性变换、仿射变换、有理函数模型等。RPC模型是卫星遥感影像的通用成像几何模型其中的一种，它适用于各类传感器包括最新的航空和航天传感器模型。基于RPC的成像模型并不要求了解传感器的实际特性和成像过程，是一种能获得和严格成像模型近似一致精度的形式简单的广义成像模型。

RPC模型具体定义关系式如下：

$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)}$

其中(X，Y)为正则化的影像坐标，(P，L，H)为正则化的地面点坐标。

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

Num_L(P，L，H)＝a₁+a₂L+a₃P+a₄H+a₅LP+a₆LH+a₇PH+a₈L²+a₉P²

+a₁₀H²+a₁₁PLH+a₁₂L³+a₁₃LP²+a₁₄LH²+a₁₅L²P+a₁₆P³+a₁₇PH²

+a₁₈L²H+a₁₉P²H+a₂₀H³

Den_L(P，L，H)＝b₁+b₂L+b₃P+b₄H+b₅LP+b₆LH+b₇PH+b₈L²+b₉P²

+b₁₀H²+b₁₁PLH+b₁₂L³+b₁₃LP²+b₁₄LH²+b₁₅L²P+b₁₆P³+b₁₇PH²

+b₁₈L²H+b₁₉P²H+b₂₀H³

Num_s(P，L，H)＝c₁+c₂L+c₃P+c₄H+c₅LP+c₆LH+c₇PH+c₈L²+c₉P²

+c₁₀H²+c₁₁PLH+c₁₂L³+c₁₃LP²+c₁₄LH²+c₁₅L²P+c₁₆P³+c₁₇PH²

+c₁₈L²H+c₁₉P²H+c₂₀H³

Den_s(P，L，H)＝d₁+d₂L+d₃P+d₄H+d₅LP+d₆LH+d₇PH+d₈L²+d₉P²

+d₁₀H²+d₁₁PLH+d₁₂L³+d₁₃LP²+d₁₄LH²+d₁₅L²P+d₁₆P³+d₁₇PH²

+d₁₈L²H+d₁₉P²H+d₂₀H³

其中，三次多项式的系数a₁，...，a₂₀，b₁，...，b₂₀，c₁，...，c₂₀是RPC文件中提供的模型参数，b₁和d₁通常为1，l，s分别为影像列数值和行数值。

所谓正则化，是一项处理无限大、发散以及一些不合理表示式的方法，其方法透过引入一项辅助性的概念——正则化因子。在RPC模型中运用这种方法对地面点和影像点的坐标进行处理。

正则化地面点坐标定义公式为

$P=\frac{\mathrm{Latitude}-\mathrm{LAT}\_\mathrm{OFF}}{\mathrm{LAT}\_\mathrm{SCALE}}$

$L=\frac{\mathrm{Longitude}-\mathrm{LONG}\_\mathrm{OFF}}{\mathrm{LONG}\_\mathrm{SCALE}}$

$H=\frac{\mathrm{Height}-\mathrm{HEIGHT}\_\mathrm{OFF}}{\mathrm{HEIGHT}\_\mathrm{SCALE}}$

其中，正则化因子LAT_OFF、LAT_SCALE、LONG_OFF、LONG_SCALE、HEIGHT_OFF、HEIGHT_SCALE是RPC文件中包含的地面点坐标正则化模型参数。Latitude表示经度、Longitude表示纬度、Height表示高程，此三项即可代表地面点的空间坐标。注：某点沿铅垂线方向到大地水准面的距离，通常称为绝对高程或海拔，本技术领域简称高程。

正则化影像坐标定义公式为

$X=\frac{\mathrm{Sample}-\mathrm{SAMP}\_\mathrm{OFF}}{\mathrm{SAMP}\_\mathrm{SCALE}}$

$Y=\frac{\mathrm{Line}-\mathrm{LINE}\_\mathrm{OFF}}{\mathrm{LINE}\_\mathrm{SCALE}}$

其中，正则化因子SAMP_OFF、SAMP_SCALE、LINE_OFF、LINE_SCALE也是RPC文件中包含的影像坐标正则化模型参数，sample代表影像列坐标，其数值即为s，line代表影像行坐标，其数值即为l。

在基于RPC模型的几何应用中首先需要基于RPC模型的正反变换。RPC正变换是指根据RPC模型的关系式，从地面控制点坐标到遥感影像平面坐标的变换过程。RPC反变换则是根据RPC模型的关系式，从原始遥感影像和影像中每个像素对应的高程数据到一定地图投影的变换过程。

但是在现有技术中，基于RPC模型的正变换都是利用最小二乘平差理论实现正变换，首先假定在地面的初始位置，然后迭代求解，由于假定的初始位置不一定准确，并且收敛的控制条件不一定是是最优的，有可能导致收敛位置的不准确。因此在RPC模型的实际应用中急切需要应用先进的科技手段，提高RPC模型正变换的可靠性，并提供一种新的RPC模型的正变换的方法以使基于RPC模型的几何处理获得更高的精度和可靠性。

发明内容

本发明目的在于解决现有技术不足，提出了一种全新的基于RPC模型反变换的正变换方法，解决了目前RPC模型正变换迭代不一定收敛的RPC模型应用中的瓶颈问题。

本发明的技术方案是：将有理多项式成像模型反变换结果引入有理多项式成像模型正变换过程进行迭代变换，即采用计算机进行以下步骤，

步骤1，设置高程数值D_hei和有理多项式成像模型参数，设定中心点经纬度坐标为(D_{lat_off}，D_{lon_off})、边长为D_{lat_scale}和D_{lon_scale}的地面矩形区域，其中4个角点的经纬度坐标为

(D_{lat_off}+D_{lat_scale}/2.0，D_{lon_off}+D_{lon_scale}/2.0)；

(D_{lat_off}-D_{lat_scale}/2.0，D_{lon_off}+D_{lon_scale}/2.0)；

(D_{lat_off}+D_{lat_scale}/2.0，D_{lon_off}-D_{lon_scale}/2.0)；

(D_{lat_off}-D_{lat_scale}/2.0，D_{lon_off}-D_{lon_scale}/2.0)；

步骤2，利用有理多项式成像模型的反变换模型计算地面矩形区域4个角点和中心点对应的影像像点坐标；

步骤3，利用上述5个点的经纬度坐标和影像像点坐标，根据仿射变换关系式 $\begin{array}{c}{D}_{\mathrm{lat}}=f_0+f_{1}x+f_{2}y\\ D_{\mathrm{lon}}=g_0+g_{1}x+g_{2}y\end{array},$ 利用最小二乘原理计算影像的近似仿射变换参数f₀，f₁，f₂，g₀，g₁，g₂，x、y分别表示正则化的影像像点坐标的列数值和行数值，D_lat、D_lon分别表示正则化的地面点坐标的纬度数值和经度数值；

步骤4，根据影像的近似仿射变换参数f₀，f₁，f₂，g₀，g₁，g₂求解量测给定的量测像点坐标(x_a，y_a)对应的近似经纬度坐标

步骤5，根据有理多项式成像模型的反变换模型计算地面点坐标

对应的理论像点坐标(x_b，y_b)；

步骤6，计算残差e，计算公式为e＝(|x_a-x_b|)²+(|y_a-y_b|)²，如果残差e满足设定的差值范围，直接输出

作为量测像点坐标(x_a，y_a)对应的经纬度坐标

否则以为中心点经纬度坐标缩小地面矩形区域范围，然后返回步骤2迭代运行，直到残差e满足设定的差值范围。

而且，步骤6.1，根据公式PS＝θ×D_{lat_scale}/sqrt((x₁-x₂)*(x₁-x₂)+(y₁-y₂)*(y₁-y₂))计算(

D_hei)位置附近的影像的像素大小PS，

其中像点坐标(x₁，y₁)为地面点(D_hei)所对应的像点坐标，像点坐标(x₂，y₂)为地面点 $(D_{{\mathrm{lat}}_a}+\mathrm{\θ}\×D_{\mathrm{lat}\_\mathrm{scale}},D_{{\mathrm{lon}}_a},D_{\mathrm{hei}})$ 所对应的像点坐标，θ为缩小取值系数；

步骤6.2，根据像素大小PS调整地面矩形区域范围，地面矩形区域4个角点的经纬度坐标为 $(D_{{\mathrm{lat}}_a}+e_s\×{\mathrm{PS},D}_{{\mathrm{lon}}_a}+e_l\×\mathrm{PS}),$ $(D_{{\mathrm{lat}}_a}-e_s\×{\mathrm{PS},D}_{{\mathrm{lon}}_a}+e_l\×\mathrm{PS}),$ $(D_{{\mathrm{lat}}_a}+e_s\×{\mathrm{PS},D}_{{\mathrm{lon}}_a}-e_l\×\mathrm{PS}),$ $(D_{{\mathrm{lat}}_a}-e_s\×{\mathrm{PS},D}_{{\mathrm{lon}}_a}-e_l\×\mathrm{PS}),$

其中e_s＝|x_a-x_b|，表示列方向的残差，e_l＝|y_a-y_b|，表示行方向的残差。

而且，缩小取值系数θ的取值为0.001～0.05。

本发明提出了一种全新的基于RPC模型反变换的正变换方法，将反变换结果引入进行迭代，克服了现有技术收敛难度，提高RPC模型正变换的可靠性，使基于RPC模型的几何处理获得更高的精度和可靠性。并且本发明的实现可在提供的现有模型参数下完成，实现方便。本发明提供的方案将可在卫星遥感影像处理的各方面工作中均到重要作用，例如影像的几何纠正、影像匹配和核线排列等方面。

附图说明

图1是本发明流程图。

具体实施方式

本发明将RPC模型反变换结果引入RPC模型正变换过程进行迭代变换，具体实施时可采用计算机程序完成，流程可参见附图。本发明的RPC模型正变化过程包括以下步骤，

步骤1，设置高程数值D_hei和RPC模型参数，设定中心点经纬度坐标为(D_{lat_off}，D_{lon_off})、边长为D_{lat_scale}和D_{lon_scale}的地面矩形区域，其中4个角点的经纬度坐标为

(D_{lat_off}+D_{lat_scale}/2.0，D_{lon_off}+D_{lon_scale}/2.0)；

(D_{lat_off}-D_{lat_scale}/2.0，D_{lon_off}+D_{lon_scale}/2.0)；

(D_{lat_off}+D_{lat_scale}/2.0，D_{lon_off}-D_{lon_scale}/2.0)；

(D_{lat_off}-D_{lat_scale}/2.0，D_{lon_off}-D_{lon_scale}/2.0)。

经纬度坐标的含义：在大地坐标系中，某点大地经度，就是通过该点的子午面与起始子午面所构成的二面角。由起始子午面起算，向东为正，叫东经，向西为负，叫西经。这点的法线与赤道面的夹角，叫做大地纬度，由赤道面起算，向北为正，叫北纬，向南为负，叫南纬。由大地经度和大地纬度所表示的地面点坐标就成为经纬度坐标。经纬度坐标加上高程即可标注地面点的坐标。为了满足下列步骤的计算条件，本发明先设置高程数值D_hei，RPC模型参数则由影像供应商提供，具体实施时输入即可。设定的地面矩形区域，需要迭代时会进行调整。

步骤2，利用RPC反变换模型计算地面矩形区域4个角点和中心点对应的影像像点坐标。

步骤3，利用上述5个点的经纬度坐标和影像像点坐标，根据仿射变换关系式 $\begin{array}{c}{D}_{\mathrm{lat}}=f_0+f_{1}x+f_{2}y\\ D_{\mathrm{lon}}=g_0+g_{1}x+g_{2}y\end{array},$ 利用最小二乘原理计算影像的近似仿射变换参数f₀，f₁，f₂，g₀，g₁，g₂，x、y分别表示正则化的影像像点坐标的列数值和行数值，D_lat、D_lon分别表示正则化的地面点坐标的纬度数值和经度数值；

仿射变换是平移、旋转、缩放等运算的组合，仿射变换关系式的系数即为仿射变换参数。

步骤4，根据影像的近似仿射变换参数f₀，f₁，f₂，g₀，g₁，g₂求解量测给定的量测像点坐标(x_a，y_a)对应的近似经纬度坐标

所谓量测像点坐标(x_a，y_a)，是具体实施时根据卫星图像量测某像点得到的，其对应的经纬度坐标

为待求解值，此步骤得到的结果为近似值，还需视能否满足精度要求决定直接输出该值为待求解值或迭代修正。

步骤5，根据RPC模型的反变换模型计算地面点坐标(D_hei)对应的理论像点坐标(x_b，y_b)；

计算理论像点坐标(x_a，y_a)，是用于将反变换结果引入正变换过程以便修正。

RPC模型的反变换模型即公式 $\begin{array}{c}X=\frac{\mathrm{Nu}{m}_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}$

在步骤4中获得的近似经纬度坐标

结合高程数值D_hei，即可分别作为P、L、H带入上式求解X、Y值，获得对应的理论像点坐标(x_b，y_b)。步骤6，计算残差e，计算公式为e＝(|x_a-x_b|)²+(|y_a-y_b|)²，如果残差e满足设定的差值范围，直接输出

作为量测像点坐标(x_a，y_a)对应的经纬度坐标

否则以

为中心点经纬度坐标缩小地面矩形区域范围，然后返回步骤2迭代运行，直到残差e满足设定的差值范围。

因为高程数值D_hei可给定，本发明通过迭代运算后获得给定高程数值对应的经纬度坐标

那么量测像点坐标(x_a，y_a)对应的地面点坐标即 $(D_{{\mathrm{lat}}_b},D_{{\mathrm{lon}}_b},D_{\mathrm{hei}}).$

本发明实施时，可根据精度需要设定差值范围，例如设定为残差e小于0.01像素，残差是指变量的真值与观测值之差。e_s＝|x_a-x_b|，表示列方向的残差，e_l＝|y_a-y_b|，表示行方向的残差，e＝(|x_a-x_b|)²+(|y_a-y_b|)²，表示距离上的残差。

为了提高每次迭代运算结果精度，保证收敛过程，本发明还提供了迭代时调整地面矩形区域的具体方法：

步骤6.1，根据公式PS＝θ×D_{lat_scale}/sqrt((x₁-x₂)*(x₁-x₂)+(y₁-y₂)*(y₁-y₂))计算(

D_hei)位置附近的影像的像素大小PS，

其中像点坐标(x₁，y₁)为地面点(

D_hei)所对应的像点坐标，并且数值为(x_b，y_b)；像点坐标(x₂，y₂)为地面点 $(D_{{\mathrm{lat}}_a}+\mathrm{\θ}\×D_{\mathrm{lat}\_\mathrm{scale}},D_{{\mathrm{lon}}_a},D_{\mathrm{hei}})$ 所对应的像点坐标，θ为缩小取值系数，sqrt表示开平方；缩小取值系数θ的取值范围可根据具体情况选取，例如0.001～0.05，本发明实施例选用0.01，则PS＝0.01×D_{lat_scale}/sqrt((x₁-x₂)*(x₁-x₂)+(y₁-y₂)*(y₁-y₂))

具体实施时，对影像的像素大小PS的计算并不限于此种方式，例如调整纬度数值，计算 $(D_{{\mathrm{lat}}_a}+0.01\×D_{\mathrm{lat}\_\mathrm{scale}},D_{{\mathrm{lon}}_a},D_{\mathrm{hei}})$ 对应的像点坐标；如果改为调整经度，相应可改成计算 $(D_{{\mathrm{lat}}_a},D_{{\mathrm{lon}}_a}+0.01\×D_{\mathrm{lat}\_\mathrm{scale}},D_{\mathrm{hei}})$ 对应的像点坐标；或者经度纬度同时调整也可。相应的，公式中分子也需要相应的改变。

步骤6.2，根据像素大小PS调整地面矩形区域范围，地面矩形区域中心点和4个角点的经纬度坐标为 $(D_{{\mathrm{lat}}_a}+e_s\×{\mathrm{PS},D}_{{\mathrm{lon}}_a}+e_l\×\mathrm{PS}),$ $(D_{{\mathrm{lat}}_a}-e_s\×{\mathrm{PS},D}_{{\mathrm{lon}}_a}+e_l\×\mathrm{PS}),$ $(D_{{\mathrm{lat}}_a}+e_s\×{\mathrm{PS},D}_{{\mathrm{lon}}_a}-e_l\×\mathrm{PS}),$ $(D_{{\mathrm{lat}}_a}-e_s\×{\mathrm{PS},D}_{{\mathrm{lon}}_a}-e_l\×\mathrm{PS}).$

Claims (2)

1.基于有理多项式成像模型的正变换方法，其特征在于：将有理多项式成像模型反变换结果引入有理多项式成像模型正变换过程进行迭代变换，即采用计算机进行以下步骤，步骤1，设置高程数值D_hei和有理多项式成像模型参数，设定中心点经纬度坐标为(D_{lat_off}，D_{lon_off})、边长为D_{lat_scale}和D_{lon_scale}的地面矩形区域，其中4个角点的经纬度坐标为

(D_{lat_off}+D_{lat_scale}/2.0，D_{lon_off}+D_{lon_scale}/2.0)；

(D_{lat_off}-D_{lat_scale}/2.0，D_{lon_off}+D_{lon_scale}/2.0)；

(D_{lat_off}+D_{lat_scale}/2.0，D_{lon_off}-D_{lon_scale}/2.0)；

(D_{lat_off}-D_{lat_scale}/2.0，D_{lon_off}-D_{lon_scale}/2.0)；

步骤2，利用有理多项式成像模型的反变换模型计算地面矩形区域4个角点和中心点对应的影像像点坐标；

步骤3，利用上述5个点的经纬度坐标和影像像点坐标，根据仿射变换关系式 $\begin{array}{c}{D}_{\mathrm{lat}}=f_0+f_{1}x+f_{2}y\\ D_{\mathrm{lon}}=g_0+g_{1}x+g_{2}y\end{array},$ 利用最小二乘原理计算影像的近似仿射变换参数f₀，f₁，f₂，g₀，g₁，g₂，x、y分别表示正则化的影像像点坐标的列数值和行数值，D_lat、D_lon分别表示正则化的地面点坐标的纬度数值和经度数值；

步骤4，根据影像的近似仿射变换参数f₀，f₁，f₂，g₀，g₁，g₂求解量测给定的量测像点坐标(x_a，y_a)对应的近似经纬度坐标

步骤5，根据有理多项式成像模型的反变换模型计算地面点坐标

对应的理论像点坐标(x_b，y_b)；

步骤6，计算残差e，计算公式为e＝(|x_a-x_b|)²+(|y_a-y_b|)²如果残差e满足设定的差值范围，直接输出作为量测像点坐标(x_a，y_a)对应的经纬度坐标

否则以为中心点经纬度坐标缩小地面矩形区域范围，然后返回步骤2迭代运行，直到残差e满足设定的差值范围，

以

为中心点经纬度坐标缩小地面矩形区域范围的具体步骤如下，步骤6.1，根据公式PS＝θ×D_{lat_scale}/sqrt((x₁-x₂)*(x₁-x₂)+(y₁-y₂)*(y₁-y₂))计算 $(D_{{\mathrm{lat}}_a},D_{{\mathrm{lon}}_a},D_{\mathrm{hei}})$ 位置附近的影像的像素大小PS，

其中像点坐标(x₁，y₁)为地面点

所对应的像点坐标，像点坐标(x₂，y₂)为地面点 $(D_{{\mathrm{lat}}_a}+\mathrm{\θ}\×D_{\mathrm{lat}\_\mathrm{scale}},D_{{\mathrm{lon}}_a},D_{\mathrm{hei}})$ 所对应的像点坐标，θ为缩小取值系数；

步骤6.2，根据像素大小PS调整地面矩形区域范围，地面矩形区域4个角点的经纬度坐标为 $(D_{{\mathrm{lat}}_a}+e_s\×\mathrm{PS},D_{{\mathrm{lon}}_a}+e_l\×\mathrm{PS}),$ $(D_{{\mathrm{lat}}_a}-e_s\×\mathrm{PS},D_{{\mathrm{lon}}_a}+e_l\×\mathrm{PS}),$ $(D_{{\mathrm{lat}}_a}+e_s\×\mathrm{PS},D_{{\mathrm{lon}}_a}-e_l\×\mathrm{PS}),$ $(D_{{\mathrm{lat}}_a}-e_s\×\mathrm{PS},D_{{\mathrm{lon}}_a}-e_l\×\mathrm{PS}),$

其中e_s＝|x_a-x_b|，表示列方向的残差，e_l＝|y_a-y_b|，表示行方向的残差。

2.根据权利要求1所述基于有理多项式成像模型的正变换方法，其特征在于：缩小取值系数θ的取值范围为0.001～0.05。

Images