CN109903352B - 一种卫星遥感影像大区域无缝正射影像制作方法 - Google Patents

Abstract

本发明涉及一种卫星遥感影像大区域无缝正射影像制作方法，包括以下内容：1)构建遥感影像成像几何模型；2)读取待匹配遥感影像和DEM数据，对DEM和待配准遥感影像数据进行匹配；3)遥感影像正射校正，得到对应的DOM产品；4)根据步骤3)中获得正射影像间的位置关系，确定影像间的重叠区域，在重叠区域匹配得到同名点，对每组同名点建立基于TPS模型的误差方程，通过迭代求解得到影像间几何偏差的改正量；5)对所有校正后正射影像进行匀光和镶嵌，得到大区域的正射影像地图产品。本发明在正射影像平差过程中，引入薄板样条模型，解决了影像间重叠范围较小导致的传统几何纠正模型精度不高的问题。

Description

一种卫星遥感影像大区域无缝正射影像制作方法

技术领域

本发明是关于一种卫星遥感影像大区域无缝正射影像制作方法，涉及摄影测量中遥感影像技术领域。

背景技术

数字正射影像图(Digital Orthophoto Map,DOM)是传统测绘4D产品的重要组成，是反映地形地貌以及地物特征的重要载体，也是遥感数据服务于各行各业的一种基础产品。生成数字正射影像的过程一般是利用数字高程模型(Digital Elevation Model,DEM)对遥感影像进行逐像元的投影差改正，消除各种畸变和位移误差，最终得到包含地理信息和各种专题的卫星遥感影像产品。原始获取的遥感影像由于成像时受传感器内部光学***畸变、扫描***非线性、外部姿态变化、地球曲率和地形起伏等的影响，均有不同程度的畸变和失真，从而导致具有一定重叠度的影像直接拼接生成的投影影像产品，在接边位置处出现几何位置不一致的情况。此外，由于成像过程中引入的各种误差，还会引起影像与已有高程数据之间的位置偏差，从而对投影差和地形构成的精度造成影响。在月球和火星轨道器摄影测量中，高精度的基准和地面控制的缺乏，造成影像之间和影像与高程数据间的不一致性问题更加突出，最终影响地图产品的工程和科研应用研究。

现有方法以已有的DOM产品为媒介，在影像和DEM之间获取控制点，通过将控制点和所有影像间的连接点同时纳入区域网平差中，利用平差解算对影像的成像几何模型进行精化，最终提高影像间和影像与DEM间的几何一致性；然后利用精化后的几何模型和DEM对影像进行正射校正。但是在处理大区域影像数据时，由于涉及的制图范围更广、包含影像数量更多、控制数据精度不满足要求等，会导致区域网平差计算量过大、解算不稳定、迭代难以收敛、平差精度不满足要求等问题。另外，如果影像间的重叠范围相对于影像的整体面积而言过小的话，传统基于多项式模型的区域网平差方法难以得到满足要求的平差精度。

发明内容

针对上述问题，本发明的目的是提供一种使得平差迭代易于收敛且有利于提高平差精度的卫星遥感影像大区域无缝正射影像制作方法。

为实现上述目的，本发明采取以下技术方案：一种卫星遥感影像大区域无缝正射影像制作方法，包括以下内容：

1)构建遥感影像成像几何模型；

2)读取待匹配遥感影像和DEM数据，对DEM和待配准遥感影像数据进行匹配；

3)遥感影像正射校正，得到对应的DOM产品；

4)根据步骤3)中获得正射影像间的位置关系，确定影像间的重叠区域，在重叠区域匹配得到同名点，对每组同名点建立基于TPS模型的误差方程，通过迭代求解得到影像间几何偏差的改正量；

5)对所有校正后正射影像进行匀光和镶嵌，得到大区域的正射影像地图产品。

进一步地，对DEM和待配准遥感影像数据进行匹配采用几何与辐射信息共同约束的DEM与影像进行匹配，具体过程为：

2.1)将DEM生成模拟影像，作为基准影像；

2.2)利用仿射-尺度不变特征变换算法，获得基准影像与遥感影像初始匹配点；

2.3)利用初始匹配点，计算有理函数模型初始精化参数；

2.4)获取最小二乘精匹配的初值；

2.5)基于全局约束的最小二乘匹配，输出待匹配的遥感影像的辐射变形参数以及几何变形参数。

进一步地，上述步骤2.1)将DEM生成模拟影像，作为基准影像，具体过程为：

2.1.1)在DEM上，逐像素计算各像素在x与y方向的梯度p和q；

2.1.2)根据梯度值，逐像素计算坡度、坡向值：

A＝p/q

式中，S为坡度值，A为坡向值；

2.1.3)根据遥感影像数据的头文件信息，得到遥感影像获取时的太阳方位角与太阳高度角，创建与DEM相同大小的空影像，逐像素计算其对应的模拟辐射值信息，并存入空影像，构建得到基准影像。

进一步地，上述步骤2.2)利用仿射-尺度不变特征变换算法，获得基准影像与遥感影像初始匹配点，具体过程为：

首先获取摄像机经度角和纬度角的采样序列，用来模拟所有可能由摄像机光轴造成的仿射变形来实现图像变换；

然后将待匹配的遥感影像进行倾斜旋转变换生成模拟影像；

最后将生成的模拟影像根据尺度不变特征变换算法进行特征点检测和基准影像匹配，获取遥感影像与基准影像间的初始匹配点。

进一步地，上述步骤2.4)获取最小二乘精匹配的初值的具体过程为：

2.4.1)利用特征点提取算子在基准影像上提取特征点；

2.4.2)利用精化的有理函数模型，将基准影像上的特征点计算到遥感影像上，作为最小二乘匹配初值。

进一步地，上述步骤2.5)基于全局约束的最小二乘匹配，输出待匹配的遥感影像的辐射变形参数以及几何变形参数，具体过程为：

2.5.1)在求解遥感影像与基准影像配准参数的同时，解算有理函数模型精化参数，引入全局约束最小二乘匹配模型如下：

g₁(x,y)+n₁(x,y)＝h₀+h₁g₂(r,c)+n₂(x,y)

r＝a^r ₀+a^r ₁r′+a^r ₂c′+δ_r

c＝b^c ₀+b^c ₁r′+b^c ₂c′+δ_c

式中，g₁与g₂为基准影像与待配准遥感影像连接点对应坐标位置的辐射值函数，n₁与n₂为基准影像与待配准遥感影像噪声的函数，(x,y)为连接点在基准影像上的坐标，h₀，h₁为基准影像与遥感影像间的辐射变形参数，a^r ₀……a^r ₂,b^c ₀……b^c ₂为几何变形参数(精化参数)，(r′,c′)为地面点利用RFM反投影到像方的投影坐标，(r,c)为反投影到像方的投影坐标精化后得到的像点坐标，δ_r,δ_c为模型误差；

2.5.2)根据全局约束最小二乘匹配模型建立误差方程：

v_g＝h₀+h₁g₂(r,c)-g₁(x,y)

v_r＝a^r ₀+a^r ₁r′+a^r ₂c′+δ_r-r

v_c＝b^c ₀+b^c ₁r′+b^c ₂c′+δ_c-c

式中，v_g为方程在配准影像间辐射值的误差，v_r为方程在列方向的误差，v_c为方程在行方向的误差；

对误差方程进行一阶泰勒展开，得到线性化的误差方程：

式中，

为泰勒级数0次项，

分别表示列方向与行方向坐标的一阶偏导数；

2.5.3)对每个输入的匹配初值点按照误差方程，逐点计算法方程；

2.5.4)求解法方程，对未知数进行改正，并判断平差后计算得出的相关系数，如果所有点相关系数的和满足要求，则输出平差结果，否则进入步骤2.5.2迭代求解；

2.5.5)输出遥感影像的辐射变形参数以及几何变形参数，完成遥感影像与基准影像间的全局最小二乘匹配。

进一步地，上述步骤3)遥感影像正射校正采用正解法和反解法。

进一步地，采用反解法对遥感影像正射校正的具体过程为：

3.1)设置要生成正射影像的X方向和Y方向分辨率(R_xdom，R_ydom)，利用原始影像大小及影像的几何模型，计算影像的覆盖范围，根据影像的范围设置影像的左上角点坐标(X₀,Y₀)，设预生成DOM上任意一点像素中心P的像点坐标为(x_dom,y_dom),计算P点对应的地面坐标(X,Y)；

X＝X₀+R_xdom*x_dom

Y＝Y₀+R_ydom*y_dom

3.2)利用(X,Y)及DEM，内插出该点的高程Z，利用建立并精化后的影像成像几何模型，计算原始图像上相应像点坐标p(x,y)；

x＝f1(X,Y,Z)

y＝f2(X,Y,Z)

式中，f1和f2代表精化后的影像成像几何模型，Z表示P点的高程；

3.3)灰度内插，求得像点p的灰度值g(x,y)；

3.4)将像点p的灰度值赋给校正后的像元，也即DOM上的P点；

3.5)依次对DOM上每个像素完成上述运算，即能获得经过地形校正的DOM产品。

进一步地，上述步骤4)根据步骤3)中获得正射影像间的位置关系，确定影像间的重叠区域，在重叠区域匹配得到同名点，然后对每组同名点建立基于TPS模型的误差方法，通过迭代求解得到影像间几何偏差的改正量，具体过程为：

4.1)根据步骤3)中获得正射影像的范围信息，确定影像间的重叠区域，其中，影像的重叠区为地理坐标矩形范围的交集，在重叠区域得到同名点；

4.2)从DOM头文件中读取其左上角点坐标和分辨率，将匹配得到的同名点坐标解求得到地面坐标：

4.3)建立基于TPS模型的正射平差误差方程，并采用最小二乘原理进行求解，得到正射影像间的几何精化参数。

进一步地，上述步骤4.3)建立基于TPS模型的正射平差误差方程，并采用最小二乘原理进行求解，得到正射影像间的几何精化参数，具体过程为：

4.3.1)建立TPS的模型：

式中，a^x ₀,a^x ₁,a^x ₂,b^y ₀,b^y ₁,b^y ₂,δ_i(i＝1,2,…n)为TPS模型参数，x,y为原始坐标，x',y'为改正后的坐标，n为同名点的数量，ψ(r_i)为径向基函数：

式中，r为影像中任意一点到匹配同名点的欧式距离；

4.3.2)基于TPS模型的误差方程对于每个解求至物方的同名点，其中，基于TPS模型的误差方程：

式中，

分别为行列方向的坐标残差，

为原始的物方坐标，

为精化后的物方坐标，i表示点标号，j表示影像编号，n表示该影像上的同名点总数；

用矩阵形式表示如下：

v＝Ax-L

L＝[Dx₁ Dy₁]^T

4.3.3)对于每张影像上的每个同名点列立如上所示的误差方程，按照最小二乘原理求解得到TPS模型的系数；

4.4)根据4.3)中求解得到的TPS模型，对每张影像进行纠正和重采样，得到校正后的正射影像。

进一步地，上述步骤4.4)对每张影像进行校正和重采样采用正解法或反解法；其中，对每张影像进行纠正和重采样采用正解法的具体过程为：

4.4.1)确定纠正后正射影像的左上角点地面坐标以及分辨率：

式中，(X₀,Y₀)为原始正射影像左上角点地面坐标，(X₀',Y₀')为纠正后正射影像的左上角点地面坐标；

4.4.2)对于原始正射影像上每一点(x_dom,y_dom)，将其转换成地面坐标(X,Y)。

4.4.3)对X,Y进行坐标校正，得到其校正后的坐标X',Y'；

4.4.4)将校正后的地面坐标转换至校正后DOM的影像坐标(x_dom',y_dom')：

x_dom'＝(X'-X₀')/m_x'

y_dom'＝(Y'-Y₀')/m_y'

式中，(m_x',m_y')为校正后DOM的分辨率,(X₀’,Y₀’)为纠正后DOM的左上角点地面坐标；

4.4.5)将(x_dom,y_dom)的灰度值赋给(x_dom',y_dom')，逐像素完成上述校正过程，最后对DOM进行灰度重采样得到纠正后的正射影像。

本发明由于采取以上技术方案，其具有以下优点：

1、本发明的大区域无缝正射影像制首先将原始遥感影像与DEM进行了配准，然后解决由于控制数据精度不足带来的正射影像间几何精度不满足要求问题，因此将未知数的解算进行了分解，减少了单次求解的未知数数量，从而使得平差迭代易于收敛，且有利于提高平差精度；

2、本发明提出了基于正射影像的平差方法，原始影像利用成像几何模型进行投影达到降维处理的效果，解决了传统平差方法中因交会角小，几何强度不足引起的平差不稳定、求解病态的问题。

3、本发明在正射影像平差过程中，引入薄板样条模型，解决了影像间重叠范围较小导致的传统几何纠正模型精度不高的问题；

4、本发明提出DEM与影像直接匹配的方法，解决了传统DEM与影像匹配过程中需要已知DOM辅助的问题，并且在DEM与影像匹配过程中同时得到影像成像几何模型的精化参数，从而使得影像正射纠正的精度和质量得到提升。

附图说明

图1为本发明卫星遥感影像大区域无缝正射影像制作方法的流程示意图；

图2是本发明实施例3×3的窗口示意图。

具体实施方式

为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都属于本发明保护的范围。

如图1所示，本发明提出的卫星遥感影像大区域无缝正射影像制作方法，具体过程为：

1、构建遥感影像成像几何模型

目前常用的遥感影像成像几何模型主要分为严格成像几何模型和通用成像几何模型。严格成像几何模型是具有严密理论基础的数学模型，它主要以共线方程为基础，可以准确地表达影像坐标与地面点空间坐标之间的严格几何关系。通用成像几何模型则回避了成像过程的复杂关系，采用数学模型来拟合像点坐标与物方点三维坐标之间的相互关系，常用的拟合模型有一般多项式模型、直接线性变换模型和有理函数模型(RationalFunction Model,RFM)等，其中，有理函数模型由于其拟合精度高、通用性好、应用方便等优点，成为了遥感影像通用几何模型中应用最广泛的一种数学模型。本发明可以采用任何严格成像几何模型和通用成像几何模型。

本发明以月球目前最高分辨率的卫星影像——LRO NAC影像的成像几何模型构建为例，详细说明本发明的严格成像几何模型和通用几何模型的构建过程。严格成像几何模型的构建一般包括内定向和外定向两个过程，而通用几何模型的构建则需要以构建的严格成像几何模型为基础，具体过程为：

1.1)LRO NAC严格成像几何模型的构建

1.1.1)LRO NAC内定向

从LRO的IK辅助文件中获取NAC相机的内定向参数例如：焦距、行列方向中心坐标、像元尺寸和畸变参数等，然后根据LRO NAC的畸变模型，对NAC相机进行内定向，其中，LRONAC的畸变模型为：

xd＝(sample-BORESIGHT_SAMPLE)*PIXEL_PITCH

r＝xd

xc＝xd/(1+k1*r²) (1)

式中，sample为像点在NAC EDR原始数据上的列坐标，BORESIGHT_SAMPLE为列方向的中心坐标，PIXEL_PITCH为列方向的像元尺寸，xd为包含畸变差的坐标(量测坐标)，k1为径向畸变参数，r为像点到主点的距离，xc为纠正后像点在焦平面的坐标，单位为mm。由于NAC为CCD线阵扫描相机，因此行方向的同类参数yd＝0，yc＝0。

1.1.2)LRO NAC外定向

(a)建立共线方程

内定向完成后，可以得到每个像素在焦平面上经过畸变改正后的坐标，外定向即建立焦平面坐标系与星固坐标系的关系，其严格成像几何模型可以用共线方程来表达：

式中，(xc,yc)是像点的焦平面坐标，f是焦距，(X,Y,Z)是对应的物方点在星固坐标系的坐标，(Xs,Ys,Zs)是摄影中心在星固坐标系的坐标，称为外方位元素的线元素，λ是一个比例因子，R是像空间坐标系到星固坐标系的旋转矩阵，由三个外方位角元素

组成。

(b)初始外方位元素的读取

要对影像进行外定向，首先需要获取成像时刻的外方位元素。外方位元素从轨道测量得到的飞行器位置和姿态数据中获得，这些测得的数据作为辅助数据存储在LRO NAC影像的SPICE kernel文件中，所以每张影像的外方位元素可以从其对应SPICE kernel中读取。

(c)内插每条扫描线的外方位元素

对于推扫式成像的轨道器影像，每一条扫描线都有相应的外方位元素。但是卫星轨道测量时间间隔大于各行影像扫描成像时间间隔，要获取每条扫描线的外方位元素需要采用内插的方式。一般采用三阶多项式建立相对于成像时间t的外方位元素函数，根据记录的每行CCD成像时间，可以插值得到每条扫描线的外方位元素。

式中，X_s(t),Y_s(t),Z_s(t)表示t时刻摄影中心在星固坐标系中的坐标，即外方位线元素；

表示t时刻焦平面在星固坐标系中的姿态角，即外方位角元素；a^X ₀……a^κ ₃表示对应参数的多项式系数，这些系数可使用最小二乘法根据轨道测量数据求解。

(d)通过共线方程和求得的外方位元素，可以将经过畸变校正的焦平面坐标转换成物方坐标，完成传感器严格成像几何模型的建立。

1.2)LRO NAC有理函数模型的建立

LRO NAC有理函数模型的建立需要首先建立虚拟控制格网，再根据生成的虚拟控制点求解有理函数模型参数。

1.2.1)虚拟控制格网的建立

建立虚拟控制格网时，需要把影像区域的高程分成若干个高程面，在像方以一定的间距生成影像的格网点坐标作为像方虚拟控制点，然后根据严格几何模型将格网点影像坐标投影到各个高程面上得到物方虚拟控制点。

1.2.2)有理函数模型参数的求解

有理多项式模型通过比值多项式建立起任意地面点坐标(lat,lon,h)和对应影像坐标(sample,line)之间的一一对应关系，表达形式如下所示：

式中，

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_i,b_i,c_i,d_i(i＝1～20)为有理函数模型参数，b₁和d₁通常为1，(P,L,H)为归一化的地面坐标，(X,Y)为归一化的影像坐标，

归一化方式如下所示：

式中，LINE_SCALE，SAMP_SCALE，SAMP_OFF和LINE_OFF为像方坐标的归一化参数；LAT_OFF,LON_OFF,HEIGHT_OFF,LAT_SCALE,LON_SCALE,HEIGHT_SCALE为物方坐标归一化参数，lat为纬度，lon为经度，h为高程。

1.2.3)由1.2.1)中得到的虚拟控制点，通过最小二乘求解78个有理函数模型参数，可以建立每张影像的有理函数模型。

2、读取待匹配遥感影像和DEM数据，对DEM和待配准遥感影像数据进行匹配。DEM与影像的匹配方式有多种，此处采用利用几何与辐射信息共同约束的DEM与影像的匹配方法进行阐述，但是并不限于该方法，方法具体过程为：

2.1)将DEM生成模拟影像，作为基准影像。

2.1.1)在DEM上，逐像素计算各像素在x与y方向的梯度：

式中，e₁，e₂…e₉为以所计算的DEM像素为中心，3×3的窗口内所对应的各个像素值，如图2所示为3×3的窗口示意图。

2.1.2)根据梯度值，逐像素计算坡度、坡向值：

A＝p/q (8)

式中，S为坡度值，A为坡向值。

2.1.3)根据遥感影像数据的头文件信息，得到遥感影像获取时的太阳方位角(Azi)与太阳高度角(Alt)，创建与DEM相同大小的空影像，逐像素计算其对应的模拟辐射值信息，并存入空影像，构建得到基准影像。

I＝255×(cosAlt×cosS+sinAlt×sinS×cosAzi-A) (9)

式中，I为计算得到的模拟基准影像各像素的模拟辐射值，Alt为遥感影像获取时的太阳高度角，Azi为遥感影像获取时的太阳方位角。

2.2)利用仿射-尺度不变特征变换(ASIFT)算法，获得基准影像与遥感影像初始匹配点。

仿射-尺度不变特征变换(ASIFT)算法为遥感技术领域的常用算法，该算法可以获得遥感影像与基准影像之间的初始匹配点，具体为：

首先获取摄像机经度角和纬度角的采样序列，用来模拟所有可能由摄像机光轴造成的仿射变形来实现图像变换；然后将待匹配的遥感影像进行倾斜旋转变换生成模拟影像；最后将生成的模拟影像根据尺度不变特征变换(SIFT)算法进行特征点检测和基准影像匹配，最终获取遥感影像与基准影像间的初始匹配点。

2.3)利用初始匹配点，计算有理函数模型初始精化参数。

本发明采用像方的仿射变换来校正反投影误差，引入精化参数(列参数以及行参数)对上述误差进行改正：

式中，(r,c)为像点坐标，(r′,c′)为地面点利用RFM反投影到像方的投影坐标，

为RFM精化参数。

2.4)获取最小二乘精匹配的初值。

2.4.1)利用特征点提取算子，例如Harris、Forstner等，在基准影像上提取特征点。

2.4.2)利用精化的有理函数模型，将基准影像上的特征点计算到遥感影像上，作为最小二乘匹配初值。

将步骤2.4.1)中提取的特征点通过步骤1得出的有理函数模型以及步骤2.3)计算得出的精化参数反投影至像方，得到特征点在影像上的像点坐标，作为全局最小二乘匹配的初值：

(a)根据基准影像的头文件信息，将步骤2.4.1)提取的基准影像特征点计算至大地坐标：

Y＝Y₀+f_y×R_y

X＝X₀+f_x×R_x (11)

式中，Y与X分别为特征点的纬度与经度坐标,Y₀与X₀分别为基准影像左上角点的纬度与经度坐标，f_y，f_x为特征点在基准影像上的行列号，R_y，R_x为基准影像在y方向以及x方向的分辨率。

根据计算得出的经纬度坐标，在数字高程模型上提取对应经纬度的高程值Z。

(b)根据步骤1及步骤2.3)计算得出的RPC及精化参数以及公式(4)、公式(10)，通过(X，Y，Z)计算得出特征点在对应待配准遥感影像上的像点坐标(r,c)；

2.5)基于全局约束的最小二乘匹配，输出待匹配的遥感影像的辐射变形参数以及几何变形参数，具体过程为：

2.5.1)在求解遥感影像与基准影像配准参数的同时，解算有理函数模型精化参数，引入全局约束最小二乘匹配模型如下：

式中，g₁与g₂为基准影像与待配准遥感影像连接点对应坐标位置的辐射值函数，n₁与n₂为基准影像与待配准遥感影像噪声的函数，(x,y)为连接点在基准影像上的坐标，h₀，h₁为基准影像与遥感影像间的辐射变形参数，a^r ₀……a^r ₂,b^c ₀……b^c ₂为几何变形参数(精化参数)，(r′,c′)为地面点利用RFM反投影到像方的投影坐标，(r,c)为反投影到像方的投影坐标精化后得到的像点坐标，δ_r,δ_c为模型误差。

2.5.2)根据全局约束最小二乘匹配模型建立误差方程如下：

式中，v_g为方程在配准影像间辐射值的误差，v_r为方程在列方向的误差，v_c为方程在行方向的误差。

对公式(13)进行一阶泰勒展开，可以得到线性化的误差方程：

式中，

为泰勒级数0次项，

分别表示列方向与行方向坐标的一阶偏导数。

在求解匹配参数及精化参数过程中输入为步骤2.2)计算得出的最小二乘匹配初值，以及以初值点为中心n×n大小的窗口内的所有点，将这些点作为输入的点，此处n取值为9。对于所有输入的点，采用线性化误差方程式(14)建立误差方程，对于步骤2.2)计算得出的初值点，采用线性化误差方程式(15)建立误差方程，最后对建立的所有误差方程进行统一求解，得出匹配参数及有理函数模型精化参数。

2.5.3)对每个输入的匹配初值点按照公式(15)中的误差方程，逐点计算法方程。

2.5.4)对法方程求解，并对未知数进行改正，并判断平差后计算得出的相关系数，例如所有点相关系数的和满足要求(本次相关系数的和小于上一次迭代的相关系数的和)，则输出平差结果，如平差精度不满足要求，则重复步骤2.5.2),2.5.3)和2.5.4)，迭代求解，其中，相关系数ρ²：

2.5.5)输出遥感影像的辐射变形参数以及几何变形参数，完成遥感影像与基准影像间的全局最小二乘匹配。

3、遥感影像正射校正

利用步骤2中求解得到的几何变形参数以及步骤1中建立的影像成像几何模型，可以对每张影像进行正射校正，得到对应的DOM产品。影像正射校正一般分为正解法和反解法，本发明以反解法为例描述具体实施过程：

3.1)计算地面点坐标。

首先设置要生成正射影像的X方向和Y方向分辨率(R_xdom，R_ydom)，该分辨率可以根据需要自行设定，也可以根据遥感影像的原始分辨率进行设置。同时，利用原始影像大小及影像的几何模型，计算影像的覆盖范围。根据影像的范围设置影像的左上角点坐标(X₀,Y₀)。设预生成DOM上任意一点像素中心P的像点坐标为(x_dom,y_dom),计算P点对应的地面坐标(X,Y)；

X＝X₀+R_xdom*x_dom

Y＝Y₀+R_ydom*y_dom

3.2)计算像点坐标

利用(X,Y)及DEM，内插出该点的高程Z，利用建立并精化后的影像成像几何模型，计算原始图像上相应像点坐标p(x,y)；

x＝f1(X,Y,Z)

y＝f2(X,Y,Z)

式中，f1和f2代表精化后的影像成像几何模型，Z表示P点的高程。

3.3)灰度内插

由于所得到的像点坐标p(x,y)不一定落在原始影像的像元中心位置，为此必须进行灰度内插，一般可采用双线性内插，求得像点p的灰度值g(x,y)；

3.4)灰度赋值

将像点p的灰度值赋给纠正后的像元，也即DOM上的P点；

3.5)依次对DOM上每个像素完成上述运算，即能获得经过地形校正的DOM产品。

4、基于薄板样条(Thin Plate Spline,TPS)模型的正射平差：根据步骤3中获得正射影像间的位置关系，确定影像间的重叠区域，在重叠区域匹配得到同名点，然后对每组同名点建立基于TPS模型的误差方法，通过迭代求解得到影像间几何偏差的改正量，对正射影像进行进一步的几何位置精化，使得影像间达到无缝的效果，具体过程为：

4.1)根据步骤3中获得正射影像的范围信息，确定影像间的重叠区域，其中，影像的重叠区为地理坐标矩形范围的交集。在重叠区域利用最小二乘匹配、SIFT、SURF、A-SIFT等匹配方法得到同名点。

4.2)从DOM头文件中读取其左上角点坐标和分辨率，根据式(11)将匹配得到的同名点坐标解求得到地面坐标：

4.3)建立基于TPS模型的正射平差误差方程，并采用最小二乘原理进行求解，得到正射影像间的几何精化参数，具体为：

4.3.1)建立TPS的模型：

式中，a^x ₀,a^x ₁,a^x ₂,b^y ₀,b^y ₁,b^y ₂,δ_i(i＝1,2,…n)为TPS模型参数，x,y为原始坐标，x',y'为改正后的坐标，n为同名点的数量，ψ(r_i)为径向基函数：

式中，r为影像中任意一点到匹配同名点的欧式距离。

4.3.2)基于TPS模型的误差方程对于每个解求至物方的同名点，其中，基于TPS模型的误差方程如下：

式中，

分别为行列方向的坐标残差，

为原始的物方坐标，

为精化后的物方坐标，i表示点标号，j表示影像编号，n表示该影像上的同名点总数。

用矩阵形式表示如下：

v＝Ax-L

L＝[Dx₁ Dy₁]^T

4.3.3)对于每张影像上的每个同名点列立如上所示的误差方程，按照最小二乘原理求解得到TPS模型的系数。

4.4)根据4.3)中求解得到的TPS模型，对每张影像进行纠正和重采样，其求解过程与步骤3中生成正射影像的过程类似，也可分为正解法和反接法，这里以正解法为例进行说明：

4.4.1)确定纠正后正射影像的左上角点地面坐标以及分辨率，一般默认与纠正前影像的分辨率相同，左上角点坐标由原DOM左上角点坐标根据TPS模型参数变换得到：

式中，(X₀,Y₀)为原始正射影像左上角点地面坐标，可以从DOM的头文件中读取得到，(X₀',Y₀')为纠正后正射影像的左上角点地面坐标。

4.4.2)对于原始正射影像上每一点(x_dom,y_dom)，利用其头文件信息根据式(11)将其转换成地面坐标(X,Y)。

4.4.3)对X,Y进行坐标纠正，得到其纠正后的坐标X',Y'；

4.4.4)将纠正后的地面坐标转换至纠正后DOM的影像坐标(x_dom',y_dom')：

x_dom'＝(X'-X₀')/m_x'

y_dom'＝(Y'-Y₀')/m_y'

式中，(m_x',m_y')为纠正后DOM的分辨率,(X₀’,Y₀’)为纠正后DOM的左上角点地面坐标。

4.4.5)将(x_dom,y_dom)的灰度值赋给(x_dom',y_dom')，逐像素完成上述纠正过程，最后对DOM进行灰度重采样得到校正后的正射影像。

5、对所有校正后正射影像的进行匀光和镶嵌，得到大区域的正射影像地图产品。该地图产品可用于地形地貌分析、地物提取以及低分辨影像的控制基准等。

根据上述说明书的揭示和教导，本发明所属领域的技术人员还可以对上述实施方式进行适当的变更和修改。因此，本发明并不局限于上面揭示和描述的具体实施方式，对本发明的一些修改和变更也应当落入本发明的权利要求的保护范围内。此外，尽管本说明书中使用了一些特定的术语，但这些术语只是为了方便说明，并不对本发明构成任何限制。

Claims (8)

1.一种卫星遥感影像大区域无缝正射影像制作方法，其特征在于，包括以下内容：

1)构建遥感影像成像几何模型；

2)读取待匹配遥感影像和DEM数据，对DEM和待配准遥感影像数据进行匹配；

3)遥感影像正射校正，得到初步校正对应的DOM产品；

4)根据步骤3)中获得正射影像间的位置关系，确定影像间的重叠区域，在重叠区域匹配得到同名点，对每组同名点建立基于TPS模型的误差方程，通过迭代求解得到影像间几何偏差的改正量，并利用求解得到影像间几何偏差的改正量对步骤3)中初步校正的DOM进行进一步校正，得到校正后的DOM，具体过程为：

4.1)根据步骤3)中获得正射影像的范围信息，确定影像间的重叠区域，其中，影像的重叠区为地理坐标矩形范围的交集，在重叠区域得到同名点；

4.2)从DOM头文件中读取其左上角点坐标和分辨率，将匹配得到的同名点坐标解求得到地面坐标：

4.3)建立基于TPS模型的正射平差误差方程，并采用最小二乘原理进行求解，得到正射影像间的几何精化参数，具体过程为：

4.3.1)建立TPS的模型：

式中，a^x ₀,a^x ₁,a^x ₂,b^y ₀,b^y ₁,b^y ₂,δ_i(i＝1,2,…n)为TPS模型参数，x,y为原始坐标，x',y'为改正后的坐标，n为同名点的数量，ψ(r_i)为径向基函数：

式中，r为影像中任意一点到匹配同名点的欧式距离；

4.3.2)基于TPS模型的误差方程对于每个解求至物方的同名点，其中，基于TPS模型的误差方程：

式中，

分别为行列方向的坐标残差，

为原始的物方坐标，

为精化后的物方坐标，i表示点标号，j表示影像编号，n表示该影像上的同名点总数；

用矩阵形式表示如下：

v＝Ax-L

L＝[Dx₁ Dy₁]^T

4.3.3)对于每张影像上的每个同名点列立如上所示的误差方程，按照最小二乘原理求解得到TPS模型的系数；

4.4)根据4.3)中求解得到的TPS模型，对每张影像进行纠正和重采样，得到校正后的正射影像；

5)对所有校正后正射影像进行匀光和镶嵌，得到大区域的正射影像地图产品。

2.根据权利要求1所述的卫星遥感影像大区域无缝正射影像制作方法，其特征在于，对DEM和待配准遥感影像数据进行匹配采用几何与辐射信息共同约束的DEM与影像进行匹配，具体过程为：

2.1)将DEM生成模拟影像，作为基准影像；

2.2)利用仿射-尺度不变特征变换算法，获得基准影像与遥感影像初始匹配点；

2.3)利用初始匹配点，计算有理函数模型初始精化参数；

2.4)获取最小二乘精匹配的初值；

2.5)基于全局约束的最小二乘匹配，输出待匹配的遥感影像的辐射变形参数以及几何变形参数。

3.根据权利要求2所述的卫星遥感影像大区域无缝正射影像制作方法，其特征在于，上述步骤2.2)利用仿射-尺度不变特征变换算法，获得基准影像与遥感影像初始匹配点，具体过程为：

首先获取摄像机经度角和纬度角的采样序列，用来模拟所有可能由摄像机光轴造成的仿射变形来实现图像变换；

然后将待匹配的遥感影像进行倾斜旋转变换生成模拟影像；

最后将生成的模拟影像根据尺度不变特征变换算法进行特征点检测和基准影像匹配，获取遥感影像与基准影像间的初始匹配点。

4.根据权利要求2所述的卫星遥感影像大区域无缝正射影像制作方法，其特征在于，上述步骤2.4)获取最小二乘精匹配的初值的具体过程为：

2.4.1)利用特征点提取算子在基准影像上提取特征点；

2.4.2)利用精化的有理函数模型，将基准影像上的特征点计算到遥感影像上，作为最小二乘匹配初值。

5.根据权利要求2所述的卫星遥感影像大区域无缝正射影像制作方法，其特征在于，上述步骤2.5)基于全局约束的最小二乘匹配，输出待匹配的遥感影像的辐射变形参数以及几何变形参数，具体过程为：

2.5.1)在求解遥感影像与基准影像配准参数的同时，解算有理函数模型精化参数，引入全局约束最小二乘匹配模型如下：

g₁(x,y)+n₁(x,y)＝h₀+h₁g₂(r,c)+n₂(x,y)

r＝a^r ₀+a^r ₁r′+a^r ₂c′+δ_r

c＝b^c ₀+b^c ₁r′+b^c ₂c′+δ_c

式中，g₁与g₂为基准影像与待配准遥感影像连接点对应坐标位置的辐射值函数，n₁与n₂为基准影像与待配准遥感影像噪声的函数，(x,y)为连接点在基准影像上的坐标，h₀，h₁为基准影像与遥感影像间的辐射变形参数，a^r ₀……a^r ₂,b^c ₀……b^c ₂为几何变形参数(精化参数)，(r′,c′)为地面点利用RFM反投影到像方的投影坐标，(r,c)为反投影到像方的投影坐标精化后得到的像点坐标，δ_r,δ_c为模型误差；

2.5.2)根据全局约束最小二乘匹配模型建立误差方程：

v_g＝h₀+h₁g₂(r,c)-g₁(x,y)

v_r＝a^r ₀+a^r ₁r′+a^r ₂c′+δ_r-r

v_c＝b^c ₀+b^c ₁r′+b^c ₂c′+δ_c-c

式中，v_g为方程在配准影像间辐射值的误差，v_r为方程在列方向的误差，v_c为方程在行方向的误差；

对误差方程进行一阶泰勒展开，得到线性化的误差方程：

v_r＝v_r0+Δa^r ₀+r′·Δa^r ₁+c′·Δa^r ₂+Δδ_r

v_c＝v_c0+Δb^c ₀+r′·Δb^c ₁+c′·Δb^c ₂+Δδ_c

式中，

为泰勒级数0次项，

分别表示列方向与行方向坐标的一阶偏导数；

2.5.3)对每个输入的匹配初值点按照误差方程，逐点计算法方程；

2.5.4)求解法方程，对未知数进行改正，并判断平差后计算得出的相关系数，如果所有点相关系数的和满足要求，则输出平差结果，否则进入步骤2.5.2迭代求解；

2.5.5)输出遥感影像的辐射变形参数以及几何变形参数，完成遥感影像与基准影像间的全局最小二乘匹配。

6.根据权利要求1所述的卫星遥感影像大区域无缝正射影像制作方法，其特征在于，上述步骤3)遥感影像正射校正采用正解法和反解法。

7.根据权利要求6所述的卫星遥感影像大区域无缝正射影像制作方法，其特征在于，采用反解法对遥感影像正射校正的具体过程为：

3.1)设置要生成正射影像的X方向和Y方向分辨率(R_xdom，R_ydom)，利用原始影像大小及影像的几何模型，计算影像的覆盖范围，根据影像的范围设置影像的左上角点坐标(X₀,Y₀)，设预生成DOM上任意一点像素中心P的像点坐标为(x_dom,y_dom),计算P点对应的地面坐标(X,Y)；

X＝X₀+R_xdom*x_dom

Y＝Y₀+R_ydom*y_dom

3.2)利用(X,Y)及DEM，内插出该点的高程Z，利用建立并精化后的影像成像几何模型，计算原始图像上相应像点坐标p(x,y)；

x＝f1(X,Y,Z)

y＝f2(X,Y,Z)

式中，f1和f2代表精化后的影像成像几何模型，Z表示P点的高程；

3.3)灰度内插，求得像点p的灰度值g(x,y)；

3.4)将像点p的灰度值赋给纠正后的像元，也即DOM上的P点；

3.5)依次对DOM上每个像素完成上述运算，即能获得经过地形校正的DOM产品。

8.根据权利要求1所述的卫星遥感影像大区域无缝正射影像制作方法，其特征在于，上述步骤4.4)对每张影像进行校正和重采样采用正解法或反解法；其中，对每张影像进行纠正和重采样采用正解法的具体过程为：

4.4.1)确定纠正后正射影像的左上角点地面坐标以及分辨率：

式中，(X₀,Y₀)为原始正射影像左上角点地面坐标，(X₀',Y₀')为纠正后正射影像的左上角点地面坐标；

4.4.2)对于原始正射影像上每一点(x_dom,y_dom)，将其转换成地面坐标(X,Y)。

4.4.3)对X,Y进行坐标校正，得到其校正后的坐标X',Y'；

4.4.4)将校正后的地面坐标转换至纠正后DOM的影像坐标(x_dom',y_dom')：

x_dom'＝(X'-X₀')/m_x'

y_dom'＝(Y'-Y₀')/m_y'

式中，(m_x',m_y')为纠正后DOM的分辨率,(X₀’,Y₀’)为纠正后DOM的左上角点地面坐标；

4.4.5)将(x_dom,y_dom)的灰度值赋给(x_dom',y_dom')，逐像素完成上述校正过程，最后对DOM进行灰度重采样得到纠正后的正射影像。

Images