CN105387846A - 卫星影像的正射校正方法和系统 - Google Patents

Abstract

本发明涉及一种卫星影像的正射校正方法，包括：划定正射校正变换范围；获取变换范围内的DEM数据，将所述DEM数据转换到地理坐标系下；获取待校正卫星影像的RPC系数，通过RPC模型估算影像的顶点像元物方坐标；提取影像对应的DEM数据；将DEM数据通过RPC模型重新计算顶点像元物方坐标；建立地理坐标系与目标平面坐标系之间的投影变换关系；建立所述目标平面坐标系与所述待校正卫星影像的像方坐标系之间的对应关系。通过采用本发明的卫星影像的正射校正方法可以减少中间的读写，提高正射校正的计算速度，减少硬盘存储空间的压力，减少硬盘损耗，适应计算区域范围的变化，满足工程化应用要求，输出平面坐标参考下卫星影像正射校正产品。

Description

卫星影像的正射校正方法和系统

技术领域

本发明涉及图像处理技术领域，特别涉及一种卫星影像的正射校正方法和系统。

背景技术

高分一号卫星(简称GF-1卫星)于2013年4月在酒泉卫星发射中心发射升空。卫星运行于太阳同步回归轨道，轨道高度645km，搭载两台2米分辨率全色/8米分辨率多光谱相机(简称PMS，组合幅宽60km)，四台16米分辨率多光谱相机(简称WFV，组合幅宽800km)，其中WFV影像实现四天的覆盖周期。其中，GF-1LV1A影像产品具备详细的辐射定标系数、光谱响应参数、有理多项式系数，以其良好的时空分辨率特性在农情资源与灾害监测等领域具有广泛的应用前景。GF-1LV1A卫星影像产品不具备空间参考，且幅宽较大，其几何畸变随着像元位置与星下点距离增大而逐渐严重，同时考虑到影像更新快速、数据量大等特征，如何快速执行正射校正，改正地形起伏和传感器误差引起的卫星影像几何畸变，成为GF-1LV1A数据处理的基础性工作之一。

正射校正模型通常分为传感器物理模型(即严格成像模型)和通用传感器模型(包括RPC模型，格网插值模型，通用实时传感器模型)两类。其中传感器物理模型以传统摄影测量中的共线条件方程为代表，综合考虑传感器成像时引起几何畸变的地形起伏、大气折射、相机透镜畸变及传感器平台的位置、姿态变化等因素，建立像方和物方的数学关系，成像理论严密，定位精度较高，但解算复杂，且因为某些商业原因传感器的核心信息和卫星轨道参数并未公开而难以广泛适用。通用传感器模型是一种依赖于传感器物理模型或地面控制点建立的独立于不同传感器的通用数学模型，直接采用数学函数形式描述地面点物方坐标和相应像方坐标间的几何关系，其中发展最快、应用最广和最具代表性的便是有理函数模型(rationalpolynomialcoefficient，RPC)。国内外许多学者采用RPC模型进行了SPOT-5卫星、QuickBird卫星、IKONOS卫星等高分辨率卫星影像的正射校正实验，得到了较高精度的实验结果。在此基础上，有学者对比了多种形式下的RPC模型参数求解精度，证明了RPC模型对SAR卫星传感器成像几何拟合的一致性；还有学者提出了“地形相关”和“地形无关”的RPC模型参数求解方式，研究了用最小二乘方法求解RPC模型参数的算法；相关的研究还包括RPC模型精度与误差补偿，控制点动态提取，基于GPU(通用计算编程模型)提高校正效率等等。

现有技术的RPC模型正射变换方法虽然能在校正精度上满足实际应用的需求，但在应用过程中需要反复的文件读、写过程，中间文件的读写过程往往会消耗一定的计算时间，这无疑降低了数据处理效率，所以提供一种快速的卫星影像的正射校正方法具有深远的意义。

发明内容

本发明所要解决的技术问题是如何提高卫星影像的正射校正效率。

为此目的，本发明提出了一种卫星影像的正射校正方法，包括：

S1：在地理坐标系下划定正射校正变换范围；获取所述变换范围内的数字高程模型数据，将所述数字高程模型数据转换到所述地理坐标系下；

S2：获取待校正卫星影像的有理函数系数，根据所述有理函数系数建立有理函数模型，通过所述有理函数模型正变换估算所述待校正卫星影像的顶点像元物方坐标；

S3：根据所述待校正卫星影像的顶点像元物方坐标，提取与所述待校正卫星影像对应的所述转换到地理坐标系下的数字高程模型数据；

S4：将与所述待校正卫星影像对应的所述数字高程模型数据通过有理函数模型正变换，重新计算所述待校正卫星影像的顶点像元物方坐标，从中提取物方坐标的最大值与最小值；

S5：根据投影变换关系，计算所述物方坐标的最大值与最小值对应的目标平面坐标；

S6：根据所述物方坐标最大值与最小值对应的目标平面坐标与所述待校正卫星影像的像元空间分辨率，建立所述目标平面坐标系与所述待校正卫星影像的像方坐标系之间的对应关系。

优选的，所述步骤S1还包括：

在所述变换范围内建立地理格网,基于所述地理格网对所述数字高程模型数据执行重采样、裁剪、镶嵌，将所述数字高程模型数据分成若干分块栅格文件，使所述地理格网的每个网格对应唯一的一个分块栅格文件。

优选的，所述步骤S3包括：

将所述顶点像元物方坐标的最大值增加0.01°，将所述顶点像元物方坐标的最小值减小0.01°，作为提取所述数字高程模型数据的边界范围，并根据所述边界范围建立高程二维数组，建立所述高程二维数组的数组元素与所述数字高程模型数据的对应关系，将所述数字高程模型数据赋值到数组。

优选的，所述高程二维数组为[xls,yls]：

xls＝[((lon_max+0.01)-(lon_min-0.01))/0.00025]+1

yls＝[((lat_max+0.01)-(lat_min-0.01))/0.00025]+1

其中，“[]”为向下取整符号，(lon_max，lat_max)为顶点像元物方坐标的最大值，(lon_min，lat_min)为顶点像元物方坐标的最小值。

优选的，所述步骤S5还包括建立地理坐标系与目标平面坐标系之间的投影变换关系，具体包括如下步骤：

在所述目标平面坐标系下针对数据处理区域建立格网；

取所述格网的至少三个顶点的坐标；

计算所述三个顶点的坐标在地理坐标系下的同名点坐标；

根据所述三个顶点的坐标以及相应的同名点坐标，计算变换系数；

根据所述变换系数建立投影逆变换线性方程组；

优选的，所述投影逆变换线性方程组为：

$L=L_a+{\∂}_x(x-x_a)+{\∂}_y(y-y_a)$

B＝B_a+θ_y(y-y_a)+θ_x(x-x_a)

其中，(B，L)，(x，y)为待投影要素分别在地理坐标系与目标平面坐标系下的坐标；θ_x，θ_y分别为不同的变换系数，所述变换系数的计算公式为：

${\∂}_x=(L_b-L_a)/(x_b-x_a);{\∂}_y=(L_d-L_a)/(y_d-y_a);$

θ_x＝(B_b-B_a)/(x_b-x_a)；θ_y＝(B_d-B_a)/(y_d-y_a)；

其中，取所述至少三个顶点的坐标A(x_a，y_a)、B(x_b，y_b)、D(x_d，y_d),且x_a＝x_d，y_a＝y_b，所述顶点坐标A(x_a，y_a)、B(x_b，y_b)、D(x_d，y_d)在地理坐标系下的同名点坐标为A’(B_a，L_a)，B’(B_b，L_b)，D’(B_d，L_d)。

优选的，所述步骤S6具体包括：

根据所述物方坐标最大值与最小值对应的目标平面坐标与所述待校正卫星影像的像元空间分辨率建立校正结果图像二维数组，每个所述校正结果图像二维数组对应所述待校正卫星影像的一个光谱波段；

针对所述校正结果图像二维数组逐元素提取元素中心点的目标平面坐标；

根据所述投影变换关系，计算所述中心点的目标平面坐标对应的地理坐标系下的物方坐标，并从所述数字高程模型数据中提取所述物方坐标位置的高程值，通过所述有理函数模型计算所述物方坐标对应的待校正卫星影像的像方坐标。

优选的，在步骤S6之后还包括：

提取与所述像方坐标相邻的若干像元值；

通过双线性内插计算所述像方坐标对应的物方坐标的像元值；

将所述像元值作为对应的所述校准结果图像二维数组的数组元素；

在目标平面坐标系下，根据所述校准结果图像二维数组创建多光谱栅格文件。

另一方面，本发明还提供了一种卫星影像的正射校正系统，包括：

重组单元，用于在地理坐标系下划定正射校正变换范围；获取所述变换范围内的数字高程模型数据，将所述数字高程模型数据转换到所述地理坐标系下；

估算单元，用于获取待校正卫星影像的有理函数系数，根据所述有理函数系数建立有理函数模型，通过所述有理函数模型正变换估算所述待校正卫星影像的顶点像元物方坐标；

提取单元，用于根据所述待校正卫星影像的顶点像元物方坐标，提取与所述待校正卫星影像对应的所述数字高程模型数据；

重新计算单元，用于将与所述待校正卫星影像对应的所述数字高程模型数据通过有理函数模型正变换，重新计算所述待校正卫星影像的顶点像元物方坐标，从中提取物方坐标最大值与最小值；

变换单元，用于建立所述地理坐标系与目标平面坐标系之间的投影变换关系，根据所述投影变换关系，计算所述物方坐标最大值与最小值对应的目标平面坐标；

校正单元，用于根据所述物方坐标最大值与最小值对应的目标平面坐标与所述待校正卫星影像的像元空间分辨率，建立所述目标平面坐标系与所述待校正卫星影像的像方坐标系之间的对应关系。

本发明针对正射校正结果需要与平面坐标系下空间数据集成应用的情况，集成DEM数据读取、RPC变换、投影变换过程于一体，一方面基于地理格网组织DEM数据并根据待校正影像范围快速读取DEM数据进入内存数组以辅助正射校正计算；另一方面将基于数值法的快速投影变换过程集成进入正射校正计算过程，通过减少中间文件的读写提高计算效率，并可以适应计算区域范围的变化，满足工程化应用要求，输出平面坐标参考下卫星影像正射校正产品；此外还可以减少硬盘存储空间的压力，减少硬盘损耗。

附图说明

通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，在附图中：

图1示出了本发明提供的卫星影像的正射校正方法的流程示意图；

图2示出了本发明正射校正变换范围内数字高程模型数据编码方案示意图；

图3示出了本发明卫星影像正射校正方法的高程二维数组提取方案示意图；

图4示出了本发明卫星影像正射校正方法的线性投影变换示意图

具体实施方式

下面将结合附图对本发明的实施例进行详细描述。

本发明针对正射校正结果需要与平面坐标参考下空间数据集成应用的情况，提供了一种可以适用于高分一号卫星上的卫星影像的正射校正方法，如图1所示，本发明包括如下步骤：S1：在地理坐标系下划定正射校正变换范围；获取所述变换范围内的数字高程模型数据，将所述数字高程模型数据转换到所述地理坐标系下；

S2：获取待校正卫星影像的有理函数系数，根据所述有理函数系数建立有理函数模型，通过所述有理函数模型正变换估算所述待校正卫星影像的顶点像元物方坐标；

S3：根据所述待校正卫星影像的顶点像元物方坐标，提取与所述待校正卫星影像对应的所述转换到地理坐标系下的数字高程模型数据；

S4：将与所述待校正卫星影像对应的所述数字高程模型数据通过有理函数模型正变换，重新计算所述待校正卫星影像的顶点像元物方坐标，从中提取物方坐标的最大值与最小值；

S5：根据投影变换关系，计算所述物方坐标的最大值与最小值对应的目标平面坐标；

S6：根据所述物方坐标最大值与最小值对应的目标平面坐标与所述待校正卫星影像的像元空间分辨率，建立所述目标平面坐标系与所述待校正卫星影像的像方坐标系之间的对应关系。

下面结合具体的实施例对本发明的技术方案展开详细的描述。

本发明以局部区域内单景GF-1LV1A多光谱影像数据为例，对卫星影像的正射校正方法进行说明。

S1：在地理坐标系下划定正射校正变换范围；获取所述变换范围内的数字高程模型(DEM)数据，将所述数字高程模型数据转换到所述地理坐标系下；优选地，S1还包括：在所述变换范围内建立地理格网,基于所述地理格网对所述数字高程模型数据执行重采样、裁剪、镶嵌，将所述数字高程模型数据分成若干分块栅格文件，使所述地理格网的每个网格对应唯一的一个分块栅格文件。

具体的，在WGS1984地理坐标系下以1°(度)为最小单元划定正射校正变换范围。如图2所示为正射校正变换范围内数字高程模型数据编码方案示意图，变换范围zone以“度”单位，其经度跨度为△lon，纬度跨度为△lat；最小经度为lon0，最大经度为lon0+△lon，最小纬度为lat0，最大纬度为lat0+△lat，均可以被1°整除。在变换范围内建立0.5°*0.5°的地理格网，将地理格网视为相对独立的DEM数据分块存储单元；针对每个地理格网边界坐标，其经度最大值、最小值必有一个可以被1°(度)整除，纬度最大值、最小值必有一个可以被1°(度)整除；依据地理格网边界坐标最小值对地理格网编码，编码方式采用“(纬度最小值*2)-(经度最小值*2)”。将DEM数据集转换到地理坐标系下，将DEM数据的栅格分辨率统一重采样到0.00025°，并应用所述地理格网划分DEM数据集，通过裁剪、镶嵌处理使每个地理格网对应唯一的、边界与地理格网相同的DEM分块栅格文件。所述DEM分块栅格文件相互独立，行、列数均为2000并且文件名称与格网编码一致。在此基础上，在执行正射校正过程中可以根据GF-1LV1A多光谱影像数据边界坐标迅速提取相交的DEM分块栅格文件。

S2：获取待校正卫星影像的有理函数(RPC)系数，根据所述有理函数系数建立有理函数模型，通过所述有理函数模型正变换估算所述待校正卫星影像的顶点像元物方坐标；具体地，针对单景GF-1LV1A多光谱影像数据提取RPC系数，建立RPC模型。通过RPC模型正变换估算影像顶点像元物方坐标E(lon_e，lat_e)，F(lon_f，lat_f)，G(lon_g，lat_g)，H(lon_h，lat_h)，从中可以提取物方坐标最大值(lon_max，lat_max)与最小值(lon_min，lat_min)。

S3：根据所述待校正卫星影像的顶点像元物方坐标，提取与所述待校正卫星影像对应的所述数字高程模型数据；其中较优的，将所述顶点像元物方坐标的最大值增加0.01°，将所述顶点像元物方坐标的最小值减小0.01°，作为提取所述数字高程模型数据的边界范围，计算与该边界范围相交的地理格网在经度方向的数量Gxn以及在纬度方向的数量Gyn，其中，Gxn与Gyn的计算公式如下：

Gyn＝[(lat_max+0.01)/0.5]-［(lat_min-0.01)/0.5]+1(1)

Gxn＝[(lon_max+0.01)/0.5]-[(lon_min-0.01)）/0.5]+1(2)

所述相交的地理格网编码

，“[]”为向下取整符号。提取文件名称与Gcode集合对应的DEM数据分块栅格文件。

根据所述待提取DEM数据的边界范围建立高程二维数组(DemArray)，如图3所示为高程二维数组提取方案示意图，所述高程二维数组对应经度方向像元数量为xls，对应纬度方向的像元数量为yls，所述高程二维数组为[xls,yls]：

xls＝[((lon_max+0.01)-(lon_min-0.01))/0.00025]+1(3)

yls＝[((lat_max+0.01)-(lat_min-0.01))/0.00025]+1(4)

其中，“[]”为向下取整符号，(lon_max，lat_max)为顶点像元物方坐标的最大值，(lon_min，lat_min)为顶点像元物方坐标的最小值。

在统一栅格分辨率的基础上，所述高程二维数组中的元素分布与所述DEM分块栅格文件中的像元分布完全重合，可以通过像方坐标建立所述高程二维数组的数组元素与所述数字高程模型数据的分块栅格文件像元的对应关系，将所述数字高程模型数据赋值到数组。如图3所示，所述DEM分块栅格文件在空间分布上无缝邻接，可以根据其在空间上的排列顺序，以西北顶点为像元起点读取DEM分块栅格文件，建立初始高程二维数组OrigDemArray[Gox,Goy]。其中，Gox＝Gxn*2500；Goy＝Gyn*2500。

OrigDemArray[j*2500+k,i*2500+w]＝FN_i,j[k,v]；

i＝{0，1，…Gyn-1}，j＝{0，1，…Gxn-1}，k＝{0，1，...2499}v＝{0，1，…2499}；其中Gxn，Gyn分别为所述待提取DEM数据边界相交的地理格网在经度方向数量、纬度方向数量；FN_i,j为所述地理格网对应的像元数组，FN_i,j对应DEM分块栅格文件名称为

在此基础上建立所述高程二维数组DemArray与所述初始高程二维数组OrigDemArray间的像元关系并给DemArray赋值。DemArray[u,w]＝OrigDemArray[iniX+u,iniY+w]；

其中 $iniX=\[\left(\right({\mathrm{lon}}_{min}-0.01)-\[\frac{1{\mathrm{on}}_{\min }-0.01}{0.5}\]*0.5)/0.00025\];$

$iniY=\[(\[\frac{{\mathrm{lon}}_{max}+0.01}{0.5}\]*0.5+1-({\mathrm{lat}}_{max}+0.01))/0.00025\].$

S4：将与所述待校正卫星影像对应的所述数字高程模型数据通过有理函数模型正变换，重新计算所述待校正卫星影像的顶点像元物方坐标，从中提取物方坐标最大值与最小值；具体的，将所述高程二维数组DemArray应用于RPC模型正变换，重新计算待校正单景GF-1LV1A多光谱影像顶点像元物方坐标K(lon_k，lat_k)，M(lon_m，lat_m)，N(lon_n，lat_n)，P(lon_p，lat_p),从中提取地理坐标最大值(lon_max’，lat_max’)与最小值(lon_min’，lat_min’)。

S5：建立所述地理坐标系与目标平面坐标系之间的投影变换关系，根据所述投影变换关系，计算所述物方坐标最大值与最小值对应的目标平面坐标；其中，建立地理坐标系与目标平面坐标系之间的投影变换关系，具体包括如下步骤：在所述目标平面坐标系下针对数据处理区域建立等公里跨度格网；每个格网经度跨度△x，纬度跨度△y，如图4所示，将每个格网视为相对独立的处理单元取所述格网的至少三个顶点的坐标；计算所述三个顶点的坐标在地理坐标系下的同名点坐标；根据所述三个顶点的坐标以及相应的同名点坐标，计算变换系数；根据所述变换系数建立投影逆变换线性方程组；较优的，可以取顶点坐标A(x_a，y_a)、B(x_b，y_b)、C(x_c，y_c)D(x_d，y_d),且x_a＝x_d，x_b＝x_c，y_a＝y_b，y_c＝y_d，所述顶点A、B、D在地理坐标系下的同名点坐标为A’(B_a，L_a)，B’(B_b，L_b)，D’(B_d，L_d)，则所述变换系数的计算公式为：

${\∂}_x=(L_b-L_a)/(x_b-x_a);{\∂}_y=(L_d-L_a)/(y_d-y_a);$

θ_x＝(B_b-B_a)/(x_b-x_a)；θ_y＝(B_d-B_a)/(y_d-y_a)；

所述投影逆变换线性方程组为：

$L=L_a+{\∂}_x(x-x_a)+{\∂}_y(y-y_a)---\left(5\right)$

B＝B_a+θ_y(y-y_a)+θ_x(x-x_a)(6)

其中，(B，L)，(x，y)为待投影要素分别在地理坐标系与目标平面坐标系下的坐标；θ_x，θ_y分别为不同的变换系数。

在此基础上，对于格网zone2内部任意点T(x_t，y_t)，其对应的地理坐标T’(B_t，L_t)的计算公式为：

$L_t=L_a+{\∂}_x(x_t-x_a)+{\∂}_y(y_t-y_a);---\left(7\right)$

B_t＝B_a+θ_y(y_t-y_a)+θ_x(x_t-x_a)(8)

S6：根据所述物方坐标最大值与最小值对应的目标平面坐标与所述待校正卫星影像的像元空间分辨率，建立所述目标平面坐标系与所述待校正卫星影像的像方坐标系之间的对应关系。其中较优的，步骤S6具体包括：根据所述物方坐标最大值与最小值对应的目标平面坐标与所述待校正卫星影像的像元空间分辨率建立n个校正结果图像二维数组Sarray[zxn,zyn]，每个所述校正结果图像二维数组对应所述待校正卫星影像的一个光谱波段；zxn＝[xmax-xmin/Sxy]+1；zyn＝[ymax-ymin/Sxy]+1；其中“[]”为向下取整符号。针对所述校正结果图像二维数组逐元素提取元素中心点的目标平面坐标；根据所述投影变换关系，计算所述中心点的目标平面坐标对应的地理坐标系下的物方坐标，并从所述数字高程模型数据中提取所述物方坐标位置的高程值，通过所述有理函数模型计算所述物方坐标对应的待校正卫星影像的像方坐标。以一个校正结果图像二维数组Sarray₀为例，设Sarray₀与所述待校正GF-1LV1A影像的第1波段对应，遍历全部数组元素：对于某一数组元素Sarray₀[x₀,y₀]，x₀＝{0，1，...zxn-1}，y₀＝{0，1，…zyn-1｝，提取元素中心点平面坐标(Xg,Yg)，应用所述基于格网的线性投影逆变换体系计算对应的地理坐标(Bg,Lg)，在此基础上从所述高程二维数组DemArray提取所述地理坐标位置的高程Hg。将(Bg,Lg,Hg)应用于RPC模型计算物方坐标对应的待校正单景GF-1LV1A影像的像方坐标(s,l)。具体的，s＝Tx*Scale+Soff；l＝Ty*Lcale+Loff；

$Tx=\frac{N_3(P,L,H)}{D_3(P.L.H)};Ty=\frac{N_1(P,L,H)}{D_1(P.L.H)};$

N₁(P,L,H)＝a₁+a₂L+a₃P+a₄H+a₅LP+a₆LH+a₇PH+a₈L²+a₉P²+a₁₀H²+a_{1 1}PLH+a₁₂L³+a₁₃LP²+a₁₄LH²+a₁₅L²P+a₁₆P³+a₁₇PH²+a₁₈L²H+a₁₉P²H+a₂₀H³

D₁(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³

N₃(P,L,H)＝c₁+c₂L+c₃P+c₄H+c₅LP+c₆LH+c₇PH+c₈L²+c₉P²+c₁₀H²+c_{1 1}PLH+c₁₂L³+c₁₃LP²+c₁₄LH²+c₁₅L²P+c₁₆P³+c₁₇PH²+c₁₈L²H+c₁₉P²H+c₂₀H³

D₃(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³

$P=\frac{Bg-D_{lat\_off}}{D_{lat\_scale}};L=\frac{Lg-D_{lon\_off}}{D_{lon\_scale}},H=\frac{Hg-D_{hei\_off}}{D_{hei\_scale}};$

其中D_{lat_scale}，D_{lat_off}，D_{lon_off}，D_{lon_scale}，D_{hei_off}，D_{hei_scale}，a₁，…，a₂₀，b₁，…，b₂₀，c₁，…，c₂₀，d₁，…，d₂₀，Scale，Soff，Lcale，Loff均为已知的所述单景GF-1LV1A多光谱影像数据RPC模型系数。N₁(P,L,H)、D₁(P,L,H)、N₃(P,L,H)、D₃(P,L,H)为RPC模型的多项式，P代表纬度值，L代表经度值，H代表高程值。

其中较优的，在步骤S6之后还可以包括：提取与所述像方坐标相邻的若干像元值；通过双线性内插计算所述像方坐标对应的物方坐标的像元值；将所述像元值作为对应的所述校准结果图像二维数组的数组元素；在目标平面坐标系下，根据所述校准结果图像二维数组创建多光谱栅格文件。具体的，从所述待校正GF-1LV1A影像的第1波段数据中，提取与所述像方坐标(s,l)相邻的若干像元值，应用双线性内插计算所述地理坐标(Bg,Lg)对应的像元值，并赋值给所述数组元素Sarray₀[x₀,y₀]。对n个校正结果图像二维数组Sarray执行上述赋值过程。在目标平面坐标参考系下，根据所述n个校正结果图像二维数组Sarray创建多光谱栅格文件，完成单景GF-1LV1A影像正射校正过程。

另一方面，采用上述的卫星影像正射校正方法，本发明还提供了一种卫星影像的正射校正系统，包括：重组单元，用于在地理坐标系下划定正射校正变换范围；获取所述变换范围内的数字高程模型数据，将所述数字高程模型数据转换到所述地理坐标系下；估算单元，用于获取待校正卫星影像的有理函数系数，根据所述有理函数系数建立有理函数模型，通过所述有理函数模型正变换估算所述待校正卫星影像的顶点像元物方坐标；提取单元，用于根据所述待校正卫星影像的顶点像元物方坐标，提取与所述待校正卫星影像对应的所述数字高程模型数据；重新计算单元，用于将与所述待校正卫星影像对应的所述数字高程模型数据通过有理函数模型正变换，重新计算所述待校正卫星影像的顶点像元物方坐标，从中提取物方坐标最大值与最小值；变换单元，用于建立所述地理坐标系与目标平面坐标系之间的投影变换关系，根据所述投影变换关系，计算所述物方坐标最大值与最小值对应的目标平面坐标；校正单元，用于根据所述物方坐标最大值与最小值对应的目标平面坐标与所述待校正卫星影像的像元空间分辨率，建立所述目标平面坐标系与所述待校正卫星影像的像方坐标系之间的对应关系。

本发明针对正射校正结果需要与平面坐标系下空间数据集成应用的情况，集成DEM数据读取、RPC变换、投影变换过程于一体，一方面基于地理格网组织DEM数据并根据待校正影像范围快速读取DEM数据进入内存数组以辅助正射校正计算；另一方面将基于数值法的快速投影变换过程集成进入正射校正计算过程，通过减少中间文件的读写提高计算效率，并可以适应计算区域范围的变化，满足工程化应用要求，输出平面坐标参考下GF-1LV1A影像正射校正产品；此外还可以减少硬盘存储空间的压力，减少硬盘损耗。

虽然结合附图描述了本发明的实施方式，但是本领域技术人员可以在不脱离本发明的精神和范围的情况下做出各种修改和变型，这样的修改和变型均落入由所附权利要求所限定的范围之内。

Claims (9)

1.一种卫星影像的正射校正方法，其特征在于，包括：

S1：在地理坐标系下划定正射校正变换范围；获取所述变换范围内的数字高程模型数据，将所述数字高程模型数据转换到所述地理坐标系下；

S2：获取待校正卫星影像的有理函数系数，根据所述有理函数系数建立有理函数模型，通过所述有理函数模型正变换估算所述待校正卫星影像的顶点像元物方坐标；

S3：根据所述待校正卫星影像的顶点像元物方坐标，提取与所述待校正卫星影像对应的所述转换到地理坐标系下的数字高程模型数据；

S4：将与所述待校正卫星影像对应的所述数字高程模型数据通过有理函数模型正变换，重新计算所述待校正卫星影像的顶点像元物方坐标，从中提取物方坐标的最大值与最小值；

S5：根据投影变换关系，计算所述物方坐标的最大值与最小值对应的目标平面坐标；

S6：根据所述物方坐标最大值与最小值对应的目标平面坐标与所述待校正卫星影像的像元空间分辨率，建立所述目标平面坐标系与所述待校正卫星影像的像方坐标系之间的对应关系。

2.根据权利要求1所述的卫星影像的正射校正方法，其特征在于，所述步骤S1还包括：

在所述变换范围内建立地理格网，基于所述地理格网对所述数字高程模型数据执行重采样、裁剪、镶嵌，将所述数字高程模型数据分成若干分块栅格文件，使所述地理格网的每个网格对应唯一的一个分块栅格文件。

3.根据权利要求1或2所述的卫星影像的正射校正方法，其特征在于，所述步骤S3包括：

将所述顶点像元物方坐标的最大值增加0.01°，将所述顶点像元物方坐标的最小值减小0.01°，作为提取所述数字高程模型数据的边界范围，并根据所述边界范围建立高程二维数组，建立所述高程二维数组的数组元素与所述数字高程模型数据的对应关系，将所述数字高程模型数据赋值到数组。

4.根据权利要求3所述的卫星影像的正射校正方法，其特征在于，所述高程二维数组为[xls,yls]：

xls＝[((lon_max+0.01)-(lon_min-0.01))/0.00025]+1

yls＝[((lat_max+0.01)-(lat_min-0.01))/0.00025]+1

其中，“[]”为向下取整符号，(lon_max，lat_max)为顶点像元物方坐标的最大值，(lon_min，lat_min)为顶点像元物方坐标的最小值。

5.根据权利要求1或2所述的卫星影像的正射校正方法，其特征在于，步骤S5还包括：建立地理坐标系与目标平面坐标系之间的投影变换关系，具体包括如下步骤：

在所述目标平面坐标系下针对数据处理区域建立格网；

取所述格网的至少三个顶点的坐标；

计算所述三个顶点的坐标在地理坐标系下的同名点坐标；

根据所述三个顶点的坐标以及相应的同名点坐标，计算变换系数；

根据所述变换系数建立投影逆变换线性方程组。

6.根据权利要求5所述的卫星影像的正射校正方法，其特征在于，所述投影逆变换线性方程组为：

$L=L_a+{\∂}_x(x-x_a)+{\∂}_y(y-y_a)$

B＝B_a+θ_y(y-y_a)+θ_x(x-x_a)

其中，(B，L)，(x，y)为待投影要素分别在地理坐标系与目标平面坐标系下的坐标； θ_x，θ_y分别为不同的变换系数，所述变换系数的计算公式为：

${\∂}_x=(L_b-L_a)/(x_b-x_a);$ ${\∂}_y=(L_d-L_a)/(y_d-y_a);$

θ_x＝(B_b-B_a)/(x_b-x_a)；θ_y＝(B_d-B_a)/(y_d-y_a)；

其中，取所述至少三个顶点的坐标A(x_a，y_a)、B(x_b，y_b)、D(x_d，y_d),且x_a＝x_d，y_a＝y_b，所述顶点坐标A(x_a，y_a)、B(x_b，y_b)、D(x_d，y_d)在地理坐标系下的同名点坐标为A’(B_a，L_a)，B’(B_b，L_b)，D’(B_d，L_d)。

7.根据权利要求1或2所述的卫星影像的正射校正方法，其特征在于，所述步骤S6具体包括：

根据所述物方坐标最大值与最小值对应的目标平面坐标与所述待校正卫星影像的像元空间分辨率建立校正结果图像二维数组，每个所述校正结果图像二维数组对应所述待校正卫星影像的一个光谱波段；

针对所述校正结果图像二维数组逐元素提取元素中心点的目标平面坐标；

根据所述投影变换关系，计算所述中心点的目标平面坐标对应的地理坐标系下的物方坐标，并从所述数字高程模型数据中提取所述物方坐标位置的高程值，通过所述有理函数模型计算所述物方坐标对应的待校正卫星影像的像方坐标。

8.根据权利要求7所述的卫星影像的正射校正方法，其特征在于，在步骤S6之后还包括：

提取与所述像方坐标相邻的若干像元值；

通过双线性内插计算所述像方坐标对应的物方坐标的像元值；

将所述像元值作为对应的所述校准结果图像二维数组的数组元素；

在目标平面坐标系下，根据所述校准结果图像二维数组创建多光谱栅格文件。

9.一种卫星影像的正射校正系统，其特征在于，包括：

重组单元，用于在地理坐标系下划定正射校正变换范围；获取所述变换范围内的数字高程模型数据，将所述数字高程模型数据转换到所述地理坐标系下；

估算单元，用于获取待校正卫星影像的有理函数系数，根据所述有理函数系数建立有理函数模型，通过所述有理函数模型正变换估算所述待校正卫星影像的顶点像元物方坐标；

提取单元，用于根据所述待校正卫星影像的顶点像元物方坐标，提取与所述待校正卫星影像对应的所述数字高程模型数据；

重新计算单元，用于将与所述待校正卫星影像对应的所述数字高程模型数据通过有理函数模型正变换，重新计算所述待校正卫星影像的顶点像元物方坐标，从中提取物方坐标最大值与最小值；

变换单元，用于建立所述地理坐标系与目标平面坐标系之间的投影变换关系，根据所述投影变换关系，计算所述物方坐标最大值与最小值对应的目标平面坐标；

校正单元，用于根据所述物方坐标最大值与最小值对应的目标平面坐标与所述待校正卫星影像的像元空间分辨率，建立所述目标平面坐标系与所述待校正卫星影像的像方坐标系之间的对应关系。