CN103063200A - 高分辨率光学卫星正射纠正影像生成方法 - Google Patents

Abstract

本发明提供了一种基于RPC模型反变换的数字正射影像生产方法，其过程包括利用有理函数模型和少量地面控制点进行原始影像高精度定向，采用数字高程模型为辅助数据，使用三次卷积法或Raised Cosine插值算法，逐像素点纠正成正射投影，从而实现高通用性、高精度、高效率的正射影像生产，解决了目前正射纠正影像生产中无法在一种方法中同时顾及效率、质量及通用性的瓶颈问题。

Description

高分辨率光学卫星正射纠正影像生成方法

技术领域

本发明属于摄影成像技术领域，特别涉及高分辨率光学卫星正射纠正影像生成方法。

背景技术

为了满足不同用户对卫星遥感影像不同的需求，本领域有专家提出对卫星影像进行几何产品分级。正射纠正影像产品是一种遥感影像最广泛的应用形式，正射纠正产品是将传感器获取的原始影像经过辐射校正、传感器校正、在合适的DEM支持下，一定控制点条件下投影到用户选择的投影和参考大地基准下，其改正了由于地形起伏而造成的像点位移，整个影像具有统一的地面分辨率，并带有地理编码。数字正射纠正的实质就是将中心投影的影像通过数字元纠正形成正射投影的过程，其原理是将影像化为很多微小的区域，根据有关的参数利用相应的构像方程式或按一定的数学模型用控制点解算，求得解算模型，然后利用数字元高程模型对原始非正射影像进行纠正，使其转换为正射影像。正射纠正通常的方法有多项式纠正法、共线方程纠正法、有理函数模型纠正法、DELAUNEY三角形法等。

多项式纠正法是实践中经常使用的一种方法，其原理比较直观，计算较为简单，特别是对地面相对平坦的情况，具有足够的纠正精度。该法的基本思想是回避成像的空间几何过程，而直接对图像变形的本身进行数字模拟。该法对各种类型传感器的纠正都是普遍适用的，它认为遥感图像的整体变形可以看作是平移、缩放、旋转、仿射、偏扭、弯曲以及更高次的基本变形的综合作用结果，因而纠正前后图像相应点之间的坐标关系中可用一个相应的多项式来表达： $\left\{\begin{array}{c}x=\underset{i=0}{\overset{m}{\mathrm{\Σ}}}\underset{j=0}{\overset{n}{\mathrm{\Σ}}}\underset{k=0}{\overset{p}{\mathrm{\Σ}}}\mathrm{aijk}{X}^i{Y}^j{Z}^k\\ y=\underset{i=0}{\overset{m}{\mathrm{\Σ}}}\underset{j=0}{\overset{n}{\mathrm{\Σ}}}\underset{k=0}{\overset{p}{\mathrm{\Σ}}}\mathrm{bijk}{X}^i{Y}^j{Z}^k\end{array}\right.,$ 式中，x、y是像点坐标，X、Y、Z是地面点坐标，a_ijk和b_ijk是待求解的多项式系数。利用多项式纠正法进行正射校正，其定位精度与地面控制点的精度、分布和数量及实际地形有关。对于地形起伏较大的地区，该方法往往得不到满意的结果，特别是当倾斜角较大时，效果更差，以SPOT影像为例，当倾斜角大于100时，就不再适合于用多项式纠正了。采用此法时，在控制点上拟合很好，但在其它点的内插值可能有明显偏离，而与相邻控制点不协调，即在某些点处产生振荡现象。

共线方程纠正法是建立在对传感器成像时的位置和姿态进行模拟和解算的基础上，由于其严格给出了成像瞬间物方空间和像方空间的几何对应关系，是对成像空间几何形态的直接描述。因而从理论上来说，该法比多项式纠正法严密，特别是该法在纠正过程中还引入了地面高程的信息，因此在地形起伏较大情况下，它比多项式法更能体现出纠正精度上的优越性。其几何校正精度是目前认为最高的。目前大部分光学遥感均采用线阵CCD传感器推扫获取地面的图像，每行影像与地面符合严格的中心投影关系，并且都有各自的外方位元素，因此图像坐标(x_i，y_i)和地面坐标(X，Y，Z)符合共线方程，其表达式为： $\left[\begin{array}{c}{x}_i\\ y_i=0\\ -f\end{array}\right]=\frac{1}{\mathrm{\λ}}{M}_i^T\left[\begin{array}{c}X-X_{\mathrm{si}}\\ Y-Y_{\mathrm{si}}\\ Z-Z_{\mathrm{si}}\end{array}\right],$ 式中，(X，Y，Z)为地面点的物方空间坐标，(X_si，Y_si，Z_si)为第i行投影中心的物方空间坐标，λ为比例因子，M_i是由第i扫描行外方位角元素毋

构成的旋转矩阵。鉴于星载CCD传感器姿态平稳，因此在一定的范围内，可以近似认为外方位元素随时间线性变化。共线方程算法理论严密，充分考虑了星载光学影像的成像几何关系，因而模型的定位精度较高。但模型的定向参数之间存在很强的相关性，导致参数求解有时不稳定。

有理函数模型(RFM)纠正法在近年来才受到普遍关注，特别是IKONOS卫星的成功发射推动了对有理函数的全面研究，国际摄影测量与遥感协会成立了专门工作组研究有关RFM的校正精度、稳定性等各方面问题。有理函数模型是各种传感器几何模型的一种更广义的表达形式。是对不同的传感器模型更为精确的表达形式。它能适用于各类传感器，包括最新的航空和航天传感器。它的缺点是模型解算复杂，运算量大，并且要求控制点数目相对较多；但其优点是由于引入较多定向参数，模拟精度很高。有理函数模型将像点坐标(r，c)表示为以相应地面点空间坐标(X、Y、Z)为自变量的多项式的比值： $\left\{\begin{array}{c}{r}_n=\frac{p_1(X_n,Y_n,Z_n)}{p_2(X_n,Y_n,Z_n)}\\ c_n=\frac{p_3(X_n,Y_n,Z_n)}{p_4(X_n,Y_n,Z_n)}\end{array}\right.,$ 式中，(r_n，c_n)和(X_n，Y_n，Z_n)分别表示像素坐标(r，c)和地面点坐标(X，Y，Z)经平移和缩放后的归一化坐标，取值位于(-1.0～1.0)之间，采用归一化坐标的目的是减少计算过程中由于数据数量级差别过大引人的舍人误差。多项式中每一项的各个坐标分量X、Y、Z的幂最大不超过3，每一项各个地面坐标分量的幂的总和也不超过3。分母项p₂，p₄的取值可以有两种情况：p₂＝p₄，p₂≠p₄，每个多项式的形式为：

$p=\underset{i=0}{\overset{m_1}{\mathrm{\Σ}}}\underset{j=0}{\overset{m_2}{\mathrm{\Σ}}}\underset{k=0}{\overset{m_3}{\mathrm{\Σ}}}{a}_{\mathrm{ijk}}{X}^i{Y}^i{Z}^i=a_0+a_{1}Z+a_{2}Y+a_{3}X+a_4\mathrm{ZY}+a_5\mathrm{ZX}+$

$a_6\mathrm{YX}+a_7{Z}^2+a_8{Y}^2+a_9{X}^2+a_{10}\mathrm{ZYX}+a_{11}{Z}^{2}Y+a_{12}{Z}^{2}X+$

$a_{13}{\mathrm{ZY}}^2+a_{14}{Y}^{2}X+a_{15}{\mathrm{ZX}}^2+a_{16}{\mathrm{YX}}^2+a_{17}{Z}^3+a_{18}{Y}^3+a_{19}{X}^3,$ 式中，a_ijk是待求解的多项式系数。有理函数模型不需要知道卫星轨道信息及成像参数，且其系数包含了各种因素的影响(传感器构造、地球曲率、大气折光等)，因此适用于大多数传感器；但高阶有理函数模型计算复杂，通常由于参数过多会导致解的不稳定性；而且模型系数没有具体物理含义，因此校正精度存在一定局限性。

此外还有直接线性变换纠正、多种改良的线性变换模型纠正、多种仿射变换纠正、不规则三角网(TIN)纠正等正射纠正方法。由于有理函数模型是一种能获得和严格成像模型近似一致精度的形式简单的广义成像模型，基于RFM模型的正射纠正方法具有模拟精度高，通用性好，应用方便，计算量小等等优点，已成为当前最重要的一种高精度正射纠正影像生成方法。

但是现有的各种纠正方法在运算速度、算法可操作性、纠正精度、不同影像覆盖地理区域类型适用性等方面各有优缺点，无法做到能同时满足各项评价指标的最优，因此在正射纠正影像的实际生产中急切需要应用先进的科技手段，提高正射纠正算法的总体性能和可靠性。

发明内容

本发明目的在于解决现有技术不足，提出了一种全新的基于RPC模型反变换的数字正射影像生产方法，实现高通用性、高精度、高效率的正射影像生产。解决了目前正射纠正影像生产中无法在一种方法中同时顾及效率、质量及通用性的瓶颈问题。

本发明的技术方案是：利用有理函数模型和少量地面控制点进行原始影像高精度定向，采用数字高程模型为辅助数据，使用三次卷积法或Raised Cosine插值算法，逐像素点纠正成正射投影，其过程包括以下步骤：

步骤1，利用有理函数模型和少量地面控制点进行原始影像高精度定向；

步骤2，根据原始影像的覆盖范围，求出正射纠正后影像的覆盖范围；

步骤3，根据正射纠正影像像元阵列数目及分布情况，建立数字高程模型数据(DEM)与物方空间的坐标映射关系，并且在数字高程模型数据(DEM)上采用双线性内插算法，内插出各像元点对应的高程；

步骤4，基于正射纠正影像与物方空间的坐标映射关系以及物方空间与原始影像的坐标映射关系，构建正射纠正影像到原始影像的坐标映射关系；

步骤5，根据正射纠正影像到原始影像的坐标映射关系，进行原始影像灰度重采样及正射纠正影像灰度赋值；

步骤6，经过逐个像元的几何位置变换和灰度重采样后，可将校正后的具有较精确地理编码的影像以规定的格式输出到文件中。

优选地，所述步骤1具体包括：

利用少量合理分布的地面控制点，采用定义在影像面的仿射变换 $\begin{array}{c}y=e_0+e_1\·\mathrm{sample}+e_2\·\mathrm{line}\\ x=f_0+f_1\·\mathrm{sample}+f_2\·\mathrm{line}\end{array},$ 式中，x、y分别表示控制点在影像上的影像坐标，sample、line分别表示控制点大地坐标经过基于有理函数模型反算后在影像上的影像坐标，以此校正隐含的系统误差。

优选地，所述步骤3具体包括以下步骤：

步骤3.1，根据公式 $\begin{array}{c}{D}_{\mathrm{east}}=a+x\×d_x\\ D_{\mathrm{north}}=b+y\×d_y\end{array},$ 式中，a、b为偏移参数，d_x、d_y为影像分辨率，将正射纠正影像上的像素坐标(x，y)变换为投影面坐标(D_east，D_north)；

步骤3.2，根据投影反变换公式

将投影面坐标(D_east D_north)反投影为大地坐标(D_lat1，D_lon1)；

步骤3.3，通过数字高程模型数据(DEM)的投影模型建立的物方空间与DEM数据的坐标映射关系(x_dem，y_dem)＝T_dem ^-1(D_lat，D_lon)，式中，T_dem ^-1表示由地面点的大地坐标(D_lat，D_lon)反算到DEM数据上像素坐标(x_dem，y_dem)的转换关系，将大地坐标(D_lat1，D_lon1)变换到DEM数据上的像素坐标(x_dem，y_dem)；

步骤3.4，如果获取的DEM数据上像点坐标(x_dem，y_dem)不位于整像素上时，采用双线性内插算进行高程值内插。

优选地，所述步骤4具体包括：

步骤4.1，根据公式 $\begin{array}{c}{D}_{\mathrm{east}}=a+x\×d_x\\ D_{\mathrm{north}}=b+y\×d_y\end{array},$ 其中，a、b为偏移参数，d_x、d_y为影像分辨率，将正射纠正影像上的像素坐标(x₀，y₀)变换为投影面坐标(D_east，D_north)；

步骤4.2，根据投影反变换公式

将投影面坐标(D_east，D_north)反投影为大地坐标(D_lat1，D_lon1)；

步骤4.3，利用步骤3中取得的数字高程模型数据(DEM)，获取正射纠正影像上的像素坐标(x，y)处对应的实际地面高程值h_fact，通过原始影像有理函数成像模型反变换公式(x，y)＝T^-1(D_lat，D_lon，h)，式中，T^-1表示由地面点地理标(D_lat，D_lon，h)反算到原始影像上像素坐标(x₁，y₁)的转换关系，将大地坐标(D_lat1，D_lon1，h_fact)变换到原始影像的像素坐标(s ample₁，line₁，h_fact)；

步骤4.4，根据步骤1建立的仿射变换公式 $\begin{array}{c}y=e_0+e_1\·\mathrm{sample}+e_2\·\mathrm{line}\\ x=f_0+f_1\·\mathrm{sample}+f_2\·\mathrm{line}\end{array},$ 将像素坐标(s ample₁，line₁，h_fact)变换为控制点纠正后的像面坐标(x₁，y₁，h_fact)。

优选地，所述步骤5中采用Raised Cosine六点法或双三次卷积插值法进行原始影像灰度重采样及正射纠正影像灰度赋值。

可见，本发明实现了高通用性、高精度、高效率的正射影像生产；利用有理函数模型和少量地面控制点进行原始影像高精度定向，采用数字高程模型为辅助数据，使用三次卷积法或Raised Cosine插值算法，逐像素点纠正成正射投影。解决了目前正射纠正影像生产中无法在一种方法中同时顾及效率、质量及通用性的瓶颈问题。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明：

图1是本发明实施例的正射纠正影像影像生产方法示意图。

具体实施方式

为了使本技术领域的人员更好地理解本发明的技术方案，并使本发明的上述目的、特征和优点能够更加明显易懂，下面结合实施例及实施例附图对本发明作进一步详细的说明。

本发明通过有理函数模型和少量地面控制点进行原始影像高精度定向，采用数字高程模型为辅助数据，使用三次卷积法或Raised Cosine插值算法，逐像素点纠正成正射投影。图1是本发明所述的正射纠正影像影像生产方法流程图，其过程包括以下步骤：

步骤1：利用少量合理分布的地面控制点，采用定义在影像面的仿射变换 $\begin{array}{c}y=e_0+e_1\·\mathrm{sample}+e_2\·\mathrm{line}\\ x=f_0+f_1\·\mathrm{sample}+f_2\·\mathrm{line}\end{array},$ 式中，x、y分别表示控制点在影像上的影像坐标，sample、line分别表示控制点大地坐标经过基于有理函数模型反算后在影像上的影像坐标，以此校正隐含的卫星遥感平台的系统误差，减少像方畸变，提高原始影像的RPC模型精度，改善原始影像定向精度。

对于带有有理多项式模型参数的原始遥感影像，可以根据RPC模型参数，假设一定的高程数值，便可求得该点的三维坐标，将其投影到某投影系中，获得其平面坐标。而根据严格成像模型求解的RPC参数的定位精度与严格成像模型的精度一致，两者都可以视为一个3D坐标到2D坐标的转换关系，里面隐含了卫星遥感平台的系统误差，产生的像方畸变基本上可以看作仿射变换。因此可以利用控制点来改善基于RPC定位精度。

可采用两种方式来利用控制点提高RPC模型的精度，一种方式是利用控制点直接对RPC参数进行校正，该方法需要使用大量的控制点来解求RPC模型中的80个参数，且参数间可能存在相关性，使求解比较困难；另外一种方式是通过少量的地面控制点来计算影像的变换参数，而不校正RPC参数，该方法需要求解像点的量测坐标与利用RPC模型计算影像坐标之间的变换关系，使用少量的控制点就足够了，理论上一般采样1～6个控制点(1个控制点可以进行像方平移变换，3个控制点可以进行像方仿射变换，6个控制点可以进行像方二次变换)。

本实施例将使用第二种方式实现原始影像定向精度的改进，分析卫星系统参数对影像几何纠正精度的影响，需要改正两类误差，一类参数纠正行方向的误差，一类参数纠正列方向的误差。其中行参数吸收轨道、姿态在行方向上的影响，列参数吸收轨道、姿态在列方向上的影响，因此采用定义在影像面的仿射变换来校正此类误差： $\begin{array}{c}y=e_0+e_1\·\mathrm{sample}+e_2\·\mathrm{line}\\ x=f_0+f_1\·\mathrm{sample}+f_2\·\mathrm{line}\end{array},$ 式中，x、y分别表示控制点在影像上的影像坐标sample、line分别表示控制点大地坐标经过基于RPC模型反算后在影像上的影像坐标。根据最小二乘平差求解影像面的仿射变换参数f₀，f₁，f₂，e₀，e₁，e₂，完成利用控制点提高RPC模型的精度。因为仿射变换为线性模型，求解参数不需要初值。

在缺少控制点(控制点小于三个)条件下，为了获得较好的精度，待求解的仿射变换参数需要分析，如果仅有一个控制点，求解偏移参数e₀和f₀来消除平移误差；当有两个控制点，求解平移参数(e₀和f₀)和line方向的系数(e₂和f₂)可以获得较高的精度。

步骤2：根据原始影像的覆盖范围，求出正射纠正后影像的覆盖范围。

若输出正射影像的范围定义不恰当，会造成影像未校正完全，同时也包含了过多的空白影像。因此，在输出影像地面边界坐标值时，应当按这样的原则来确定，即由它们确定的范围，是包括校正后影像在内的最小长方形的范围。具体而言，为了确定输出影像的边界范围，需对待校正原始影像的4个角点a，b，c，d按照校正模型进行坐标变换得到4对地面坐标，从而分别找出X轴的最大值x₂和最小值x₁，Y轴的最大值y₂和最小值y₁，于是可确定输出影像的范围。鉴于数字影像是以阵列形式表示，我们可令地面坐标为(X1，Y1)的点为输出影像的第一行第一列，ΔX，ΔY为相应的地面分辨率，其反映了每个输出影像的像元所对应的地面尺寸，故可得到影像阵列的总行数M和总列数N：

$\left\{\begin{array}{c}M=\frac{(Y_2-Y_1)}{\mathrm{\ΔY}}+1\\ N=\frac{(X_2-X_1)}{\mathrm{\ΔX}}+1\end{array}\right.$

当输出影像阵列某点的地面坐标为(X，Y)，则该点在输出影像的行号I和列号J可通过下式求得

$\left\{\begin{array}{c}I=\frac{(Y_2-Y)}{\mathrm{\ΔY}}+1\\ J=\frac{(X-X_1)}{\mathrm{\ΔX}}+1\end{array}\right.$

步骤3：根据正射纠正影像像元阵列数目及分布情况，在数字高程模型数据(DEM)上采用双线性内插算法内插出各像元点对应的高程。使得正射纠正影像像元阵列和DEM像元阵列一一对应，便于后续计算，包括以下步骤：

步骤3.1，根据公式 $\begin{array}{c}{D}_{\mathrm{east}}=a+x\×d_x\\ D_{\mathrm{north}}=b+y\×d_y\end{array},$ 式中，a、b为偏移参数，d_x、d_y为影像分辨率，将正射纠正影像上的像素坐标(x，y)变换为投影面坐标(D_east，D_north)。

由于影像是以一定的地面分辨率进行排列的，所以正射纠正影像的像面坐标与其投影坐标间也具有线性的关系，该关系可以表达为下面的形式： $\begin{array}{c}{D}_{\mathrm{east}}=a+x\×d_x\\ D_{\mathrm{north}}=b+y\×d_y\end{array},$ 式中，a，b，d_x，d_y均可从正射纠正影像的影像信息中获取。a，b为正射纠正影像左上点的大地坐标，d_x，d_y为X、Y方向的分辨率。由于生成正射纠正影像时，记录了a，b，d_x，d_y，因此可以直接利用这些几何参数进行坐标转换。根据该公式即可通过像面坐标求出投影坐标。

步骤3.2，根据投影反变换公式

将投影面坐标(D_east，D_north)反投影为大地坐标(D_lat1，D_lon1)。

大地坐标也叫大地经纬度坐标，其含义：在大地坐标系中，某点大地经度，就是通过该点的子午面与起始子午面所构成的二面角。由起始子午面起算，向东为正，叫东经，向西为负，叫西经。这点的法线与赤道面的夹角，叫做大地纬度，由赤道面起算，向北为正，叫北纬，向南为负，叫南纬。由大地经度和大地纬度所表示的地面点坐标就成为经纬度坐标。经纬度坐标加上高程即可标注地面点的坐标。

步骤3.3，通过数字高程模型数据(DEM)的投影模型建立的物方空间与DEM数据的坐标映射关系(x_dem，y_dem)＝T_dem ^-1(D_lat，D_lon)，式中，T_dem ^-1表示由地面点的地理标(D_lat，D_lon)反算到DEM数据上像素坐标(x_dem，y_dem)的转换关系，将大地坐标(D_lat1，D_lon1)变换到DEM数据上的像素坐标(x_dem，y_dem)。

步骤3.4，如果获取的DEM数据上像点坐标(x_dem，y_dem)不位于整像素上时，采用双线性内插算进行高程值内插。

双线性内插法最适合基于格网的内插，其函数公式为：

H(P)＝(1-Δx)(1-Δy)H₁₁+(1-Δx)ΔyH₁₂+Δx(1-Δy)H₂₁+ΔxΔyH₂₂

式中：H(P)为P点的灰度值，H₁₁，H₁₂，H₂₁，H₂₂分别为最靠近待插值点的4个已知点的高程值，Δx＝x-INT(x)，Δy＝y-INT(y)，x、y为待插值点平面坐标。

步骤4：基于正射纠正影像的投影模型建立的正射纠正影像与物方空间的坐标映射关系，以及原始影像的有理多项式模型反变换和步骤1中的仿射变换模型共同建立的物方空间与原始影像的坐标映射关系，构建正射纠正影像到原始影像的坐标映射关系。

构建正射纠正影像和原始影像间的像元几何位置关系，常用的有“直接法”和“间接法”，直接变换法是从原始影像出发，依次对其中的每一像元(x，y)依据选定的校正模型F(x，y)求其在新影像中的地面坐标(X，Y)，同时把(x，y)上的灰度值送到地面坐标为(X，Y)所对应的输出影像像点。F(x，y)是以“直接法”形式表现的校正模型。间接变换法是从校正后的影像阵列出发，依次计算每个像元在原始影像中的像点位置(x，y)，最后把(x，y)处灰度值赋给地面坐标为(X，Y)的输出影像像点，以“间接法”形式表现的校正模型G(x，y)是F(x，y)的逆函数。相比较而言，直接法在某些校正变换中解算复杂，同时校正后的影像，有些像素可能重叠，有些像素之间可能出现空白，因此，一般很少被采用；而间接法通过灰度重采样获取校正后的像元灰度值，在实际中被广泛应用，本实施例将采用间接法。包括以下步骤：

步骤4.1，根据公式 $\begin{array}{c}{D}_{\mathrm{east}}=a+x\×d_x\\ D_{\mathrm{north}}=b+y\×d_y\end{array},$ 其中，a b为偏移参数，d_x d_y为影像分辨率，将正射纠正影像上的像素坐标(x₀，y₀)变换为投影面坐标(D_east，D_north)。该步骤在步骤3.1中已经有详细的描述，在此不再重复叙述。

步骤4.2，根据投影反变换公式

将投影面坐标(D_east，D_north)反投影为大地坐标(D_lat1，D_lon1)。该步骤在步骤3.2中已经有详细的描述，在此不再重复叙述。

步骤4.3，获取正射纠正影像上的像素坐标(x，y)处对应的实际地面高程值h_fact(利用步骤3处理后的DEM数据获取)，通过原始影像有理函数成像模型反变换公式(x，y)＝T^-1(D_lat，D_lon，h)，式中，T^-1表示由地面点地理标(D_lat，D_lon，h)反算到原始影像上像素坐标(x₁，y₁)的转换关系，将大地坐标(D_lat1，D_lon1，h_fact)变换到原始影像的像素坐标(s ample₁，line₁，h_fact)。

有理函数模型将地面点大地坐标与其对应的像点坐标用比值多项式关联起来。为了增强参数求解的稳定性，将地面坐标和影像坐标正则化到-1和1之间。

其比值多项式定义为： $\begin{array}{c}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)}\end{array},$ 式中，(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)为三次多项式。基于有理函数成像模型的反变换的具体实施步骤如下：

1)根据公式 $\begin{array}{c}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}}\end{array},$ 式中，正则化因子LAT_OFF、LAT_SCALE、LONG_OFF、LONG_SCALE、HEIGHT_OFF、HEIGHT_SCALE是RPC文件中包含的地面点坐标正则化模型参数。Latitude表示经度、Longitude表示纬度、Height表示高程，将地面坐标(D_lat1，D_lon1，h_fact)正则化为(P₁，L₁，H₁)；

2)由公式 $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)},$ 计算地面坐标(P₁，L₁，H₁)的像点正则化坐标(X，Y)；

3)由公式 $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，将像点正则化坐标(X，Y)转换成(s ample₁，line₁，h_fact)坐标。

步骤4.4，根据步骤1建立的仿射变换公式 $\begin{array}{c}y=e_0+e_1\·\mathrm{sample}+e_2\·\mathrm{line}\\ x=f_0+f_1\·\mathrm{sample}+f_2\·\mathrm{line}\end{array},$ 将像素坐标(s ample₁，line₁，h_fact)变换为控制点纠正后的像面坐标(x₁，y₁，h_fact)。

步骤5：根据正射纠正影像到原始影像的坐标映射关系，采用RaisedCosine六点法或双三次卷积插值法进行原始影像灰度重采样及正射纠正影像灰度赋值。

如果正射纠正影像上的像元坐标(x₀，y₀，h_fact)在原始影像中的相应位置的坐标值(x₁，y₁，h_fact)不为整数时，则必须利用在原始影像中该点附近的若干像元的灰度，并考虑这些像元对其所作的贡献和影响，进行灰度重采样。

Raised Cosine六点算法：没有码元窜扰的脉冲形状对于最佳传送是十分重要的，Raised Cosine函数就是满足Nyquist频率标准的。其插值公式：

$i\left(x\right)=\sin c\left(x\right)\frac{\cos \left(\mathrm{\πvx}\right)}{1-4v^2{x}^2}\mathrm{rect}\left(\frac{x}{L}\right)$

式中，L表示窗口大小，

r表示采样频率与Nyquist频率之比，x表示采样点与格网点的水平或垂直距离，sin c函数形式为

最临近像素法rect函数形式为 $\mathrm{rect}\left(x\right)=\left\{\begin{array}{c}0,\left|x\right|>\frac{1}{2}\\ \frac{1}{2},\left|x\right|=\frac{1}{2}\\ 1,\left|x\right|>\frac{1}{2}\end{array}\right..$ 实验证明该算法，与其它的算法相比，即便是使用较小的窗口，该算法的误差也比其他算法的误差小。理论上，Raised Cosine六点算法的相位误差比十二点的Knab SW算法的相位误差的四分之一还小。

三次卷积法：其基本思想是增加邻点来获得最佳插值函数。取与计算点(x，y)周围相邻的16个点，先在某一个方向上内插，如先在x方向上，每4个值依次内插4次，求出f(i，j-1)，f(i，j)，f(i，j+1)，f(i，j+2)，再根据这四个计算结果在y方向上内插，得到f(x，y)。每一组4个样点组成一个连续内插函数。可以证明，这种三次多项式内插过程实际上是一种卷积运算，故称为三次卷积内插。

在x方向上，设m＝j-1，j，j+1，j+2。计算公式为：

f(x，m)＝α²(α-1)f(i+2，m)+α(1+α-α²)f(i+1，m)+

(1-2α²+α³)f(i，m)-α(1-α²)f(i-1，m)

令α＝x-i，β＝y-i，i＝Integer(x)，j＝Integer(y)，共计算四次，再用同样方法计算y方向的值：

f(x，y)＝β²(β-1)f(x，j+2)+β(1+β-β²)f(x，j+1)+

(1-2β²+β³)f(x，j)-β(1-β²)f(x，j-1)

需注意的是，欲以三次卷积内插获得好的图像效果，就要求位置校正过程更准确，即对控制点选取的均匀度要求更高。

步骤6：经过逐个像元的几何位置变换和灰度重采样后，可将校正后的具有较精确地理编码的影像以规定的格式输出到文件中，完成原始影像正射纠正处理。

以上所述，仅为本发明的具体实施方式，本发明还可以应用在其它设备中；以上描述中的尺寸和数量均仅为参考性的，本领域技术人员可根据实际需要选择适当的应用尺寸，而不脱离本发明的范围。本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求所界定的保护范围为准。

Claims (5)

1.一种高分辨率光学卫星正射纠正影像生成方法，其特征在于，包括以下步骤：

步骤1，利用有理函数模型和少量地面控制点进行原始影像高精度定向；

步骤2，根据原始影像的覆盖范围，求出正射纠正后影像的覆盖范围；

步骤3，根据正射纠正影像像元阵列数目及分布情况，建立数字高程模型数据(DEM)与物方空间的坐标映射关系，并且在数字高程模型数据(DEM)上采用双线性内插算法，内插出各像元点对应的高程；

步骤4，基于正射纠正影像与物方空间的坐标映射关系以及物方空间与原始影像的坐标映射关系，构建正射纠正影像到原始影像的坐标映射关系；

步骤5，根据正射纠正影像到原始影像的坐标映射关系，进行原始影像灰度重采样及正射纠正影像灰度赋值；

步骤6，经过逐个像元的几何位置变换和灰度重采样后，可将校正后的具有较精确地理编码的影像以规定的格式输出到文件中。

2.根据权利要求1所述的生成方法，其特征在于，所述步骤1具体包括：

利用少量合理分布的地面控制点，采用定义在影像面的仿射变换 $\begin{array}{c}y=e_0+e_1\·\mathrm{sample}+e_2\·\mathrm{line}\\ x=f_0+f_1\·\mathrm{sample}+f_2\·\mathrm{line}\end{array},$ 式中，x、y分别表示控制点在影像上的影像坐标，sample、line分别表示控制点大地坐标经过基于有理函数模型反算后在影像上的影像坐标，以此校正隐含的系统误差。

3.根据权利要求2所述的生成方法，其特征在于，所述步骤3具体包括以下步骤：

步骤3.1，根据公式 $\begin{array}{c}{D}_{\mathrm{east}}=a+x\×d_x\\ D_{\mathrm{north}}=b+y\×d_y\end{array},$ 式中，a b为偏移参数，d_x d_y为影像分辨率，将正射纠正影像上的像素坐标(x，y)变换为投影面坐标(D_east，D_north)；

步骤3.2，根据投影反变换公式

将投影面坐标(D_east，D_north)反投影为大地坐标(D_lat1，D_lon1)；

步骤3.3，通过数字高程模型数据(DEM)的投影模型建立的物方空间与DEM数据的坐标映射关系(x_dem，y_dem)＝T_dem ^-1(D_lat，D_lon)，式中，T_dem ^-1表示由地面点的大地坐标(D_lat D_lon)反算到DEM数据上像素坐标(x_dem，y_dem)的转换关系，将大地坐标(D_lat1，D_lon1)变换到DEM数据上的像素坐标(x_dem，y_dem)；

步骤3.4，如果获取的DEM数据上像点坐标(x_dem，y_dem)不位于整像素上时，采用双线性内插算进行高程值内插。

4.根据权利要求3所述的生成方法，其特征在于，所述步骤4具体包括：

步骤4.1，根据公式 $\begin{array}{c}{D}_{\mathrm{east}}=a+x\×d_x\\ D_{\mathrm{north}}=b+y\×d_y\end{array},$ 其中，a b为偏移参数，d_x d_y为影像分辨率，将正射纠正影像上的像素坐标(x₀，y₀)变换为投影面坐标(D_east，D_north)；

步骤4.2，根据投影反变换公式将投影面坐标(D_east，D_north)反投影为大地坐标(D_lat1，D_lon1)；

步骤4.3，利用步骤3中取得的数字高程模型数据(DEM)，获取正射纠正影像上的像素坐标(x，y)处对应的实际地面高程值h_fact，通过原始影像有理函数成像模型反变换公式(x，y)＝T^-1(D_lat D_lon，h)，式中，T^-1表示由地面点地理标(D_lat，D_lon，h)反算到原始影像上像素坐标(x₁，y₁)的转换关系，将大地坐标(D_lat1，D_lon1，h_fact)变换到原始影像的像素坐标(s ample₁，line₁，h_fact)；

步骤4.4，根据步骤1建立的仿射变换公式 $\begin{array}{c}y=e_0+e_1\·\mathrm{sample}+e_2\·\mathrm{line}\\ x=f_0+f_1\·\mathrm{sample}+f_2\·\mathrm{line}\end{array},$ 将像素坐标(s ample₁，line₁，h_fact)变换为控制点纠正后的像面坐标(x₁，y₁，h_fact)。

5.根据权利要求4所述的生成方法，其特征在于，所述步骤5中采用RaisedCosine六点法或双三次卷积插值法进行原始影像灰度重采样及正射纠正影像灰度赋值。

Images