CN108830889A - 基于全局几何约束的遥感影像与基准影像的匹配方法 - Google Patents

Abstract

本发明涉及一种基于全局几何约束的遥感影像与基准影像的匹配方法，其特征在于包括以下内容：构建遥感影像成像几何模型；读取参考影像数据和待匹配遥感影像，对参考影像数据预处理得到基准影像；对待匹配遥感影像进行预处理；将基准影像与待配准遥感影像进行仿射‑尺度不变特征变换匹配，提取特征点，并获取最小二乘匹配的初值；基于最小二乘匹配的初值通过全局约束最小二乘匹配模型进行最小二乘匹配，输出待配准遥感影像的辐射变形参数以及几何变形参数。

Description

基于全局几何约束的遥感影像与基准影像的匹配方法

技术领域

本发明是关于一种基于全局几何约束的遥感影像与基准影像的匹配方法，涉及摄影测量与遥感影像处理技术领域。

背景技术

在航空航天及深空摄影测量中，利用遥感影像进行几何定位，往往需要将影像纠正到统一的基准。遥感影像的纠正中需要建立几何模型，即为构建影像坐标与基准坐标系间的关系。几何模型可以由遥感影像获取时的轨道、姿态等测控参数获取，同时也可以由影像与基准数据间匹配建立的对应关系拟合获取。由于各种测量的系统误差及随机误差，例如卫星轨道和姿态测量误差、相机畸变等，导致利用测控星历数据获取的遥感影像几何模型计算存在一定偏差，使得卫星影像直接纠正后，影像之间存在着较大的空间不一致性及整体地形偏移等问题。例如不同传感器获取的多源遥感影像间、同一传感器不同时间获取的遥感影像间都存在着几何不一致性。因此需要将影像与空间参考数据进行配准，获取控制点进行影像模型的精化。遥感影像与基准数据的配准，可以为影像高精度几何处理提供精确的地面控制点，为几何模型的精化，仪器校正，自动校正等领域服务。同时也可以直接作为控制点进行遥感影像的纠正，因此为了充分利用多源遥感影像，实现不同数据源在定位精度、分辨率、覆盖范围等各方面的优势互补，将这些影像与参考影像进行配准，使其统一到一致的空间参考中则显得尤为重要。

传统的遥感影像与参考影像间的匹配，往往只顾及影像的灰度信息，例如利用SIFT、特征提取与匹配等常规的方案进行影像与影像间的匹配。在影像已经具有初始几何模型时，几何模型往往只用于获取初始匹配点或限定匹配搜索范围，并没有充分利用几何信息的约束来提高匹配的精度。在参考影像分辨率较高、影像质量较好情况下，遥感影像与参考影像间的匹配精度往往能达到较高精度，此时几何模型内部一致性约束下的精度低于影像间的匹配精度，几何模型全局约束不会提高匹配精度。然而在参考影像或参考数据信息量不足，分辨率较遥感影像更低的时候，遥感影像自身几何模型的一致性精度会优于影像匹配定位精度，此时利用遥感影像自身几何模型的全局约束进行匹配，将能提高匹配精度。

发明内容

针对上述问题，本发明的目的是提供一种能够提高遥感影像与参考数据的匹配精度的基于全局几何约束的遥感影像与基准影像的匹配方法。

为实现上述目的，本发明采取以下技术方案：一种基于全局几何约束的遥感影像与基准影像的匹配方法，其特征在于包括以下内容：构建遥感影像成像几何模型；读取参考影像数据和待匹配遥感影像，对参考影像数据预处理得到基准影像；对待匹配遥感影像进行预处理；将基准影像与待配准遥感影像进行仿射-尺度不变特征变换匹配，提取特征点，并获取最小二乘匹配的初值；基于最小二乘匹配的初值通过全局约束最小二乘匹配模型进行最小二乘匹配，输出待配准遥感影像的辐射变形参数以及几何变形参数。

进一步地，当参考影像数据只有DEM时，对DEM进行预处理构建基准影像，具体过程为：

2.1)在DEM上，逐像素计算在x与y方向的梯度；

2.2)根据梯度值，逐像素计算坡度S、坡向值A；

2.3)根据参考影像数据的头文件信息，创建与DEM相同大小的空影像，逐像素计算其对应的模拟辐射值信息并存入空影像，构建得到基准影像：

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

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

进一步地，当参考影像数据为DOM时，DOM即为基准影像。

进一步地，对待配准遥感影像进行预处理的具体过程为：

2.1)计算待配准遥感影像的直方图p_r(r_k)，具体过程为：

a)将待配准遥感影像的灰度级归一化：

式中，L为灰度级层次数，r_k为遥感影像第k级的灰度，k为遥感影像的灰度级别；

b)计算各灰度级出现的概率p_r(r_k)：

式中，n_k为像素值为k的像素的频数，p_r(r_k)为其出现的概率，N为遥感影像的像素总数。

2.2)将r_k根据变换函数T(r_k)进行变换得到的灰度级s_k，将s_k四舍五入为范围[0,L-1]内的整数：

式中，n_j为灰度值为r_j的像素数；

2.3)将步骤2.1)和2.2)计算出的s_k与G(z_q)进行匹配，具体过程为：

A)采用步骤2.1)计算基准影像的灰度级直方图p_z(z_q)，并对q＝0,1,2,…,L-1计算变换函数G的所有值，把G的值四舍五入为范围[0,L-1]内的整数，并将G存储在一个表中：

式中，z_q为基准影像直方图匹配后第q级的灰度，z_i为基准影像直方图匹配后第i级的灰度；

B)对每一个值s_k,k＝0,1,2,…,L-1,使用步骤A)存储的G值寻找相应的z_q值，以使G(z_q)最接近s_k，并存储这些从s到z的映射，当满足给定s_k的z_q值多于一个时，按惯例选择最小的值；

C)使用步骤B)找到的映射将待配准遥感影像的每个变换后的像素值s_k映射为基准影像的相应z_q值，形成直方图匹配后的待配准遥感影像。

进一步地，计算待配准遥感影像所有灰度级s_k的具体过程为：

a)将待配准遥感影像的灰度级归一化，若遥感图像的灰度级为0,1，…，L-1，则令：

式中，L为灰度级层次数，r_k为遥感影像第k级的灰度，k为遥感影像的灰度级别；

b)计算各灰度级的像素频数

假设n_k为灰度级为r_k的像素个数，N为总的像素个数，令：

式中，n_k为像素值为k的像素的频数，p_r(r_k)为其出现的概率；

c)计算第k级的灰度级

式中，n_j为灰度值为r_j的像素数。

进一步地，将基准影像与待配准遥感影像进行仿射-尺度不变特征变换匹配，提取特征点，并获取最小二乘匹配的初值的具体过程为：

3.1)采用仿射-尺度不变特征变换算法获得待配准遥感影像与基准影像之间的初始匹配点；

3.2)根据初始匹配点结合DEM计算有理函数模型精化参数

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

3.3)在基准影像上提取特征点；

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

进一步地，在基准影像上提取特征点的具体过程为：

1)确定一个n×n大小的影像窗口，对窗口内的每一个像素点进行一阶差分运算，求得在x，y方向的梯度g_x，g_y；

2)对梯度值进行高斯滤波：

式中，Ga为高斯卷积模板，σ取0.4，g_x为灰度在x方向的梯度，g_y为y方向的梯度；

3)计算矩阵M及兴趣值I：

I＝det(M)-k·tr²(M)

式中，det是矩阵的行列式，tr是矩阵的迹，k是默认常数，取0.04；

4)选取兴趣值的局部极值点，并在设定窗口内取最大值作为特征点。

进一步地，基于最小二乘匹配初值通过全局约束最小二乘匹配模型进行最小二乘匹配，输出待配准遥感影像的辐射变形参数以及几何变形参数的具体过程为：

4.1)引入全局约束最小二乘匹配模型：

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

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

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

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

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

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

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

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

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

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

式中，为泰勒级数0次项，分别表示列方向与行方向坐标的一阶偏导数；

即针对每个特征点，误差方程为：

v＝Bx-L

其中，

x＝(Δa₀ Δa₁ Δa₂ Δb₀ Δb₁ Δb₂ Δh_0i Δh_1i Δδ_ri Δδ_ci)^T

式中，m为连接点的点对数，n为窗口大小；

4.3)对每个输入的匹配初值点代入线性化的误差方程，逐点计算法方程；

4.4)根据最小二乘原理通过对x＝(B^TB)^-1B^TL进行求解可以计算出所有未知数，并根据相关系数计算公式计算相关系数ρ²：

判断平差后计算得出的相关系数，如所有点相关系数的和满足平差精度要求则输出平差结果，平差精度不满足要求，则重复步骤4.2)～4.4)，迭代求解；

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

本发明由于采取以上技术方案，其具有以下优点：1、本发明提出了基于全局约束的遥感影像与基准影像的最小二乘配准方法，由于全局几何约束的引入，可以提升遥感影像与基准影像匹配的可靠性和精度。2、在匹配的同时，本发明对遥感影像几何模型进行了精化，使得本发明具有更良好的扩展性和更广泛的使用范围。综上所述，本发明可以广泛应用于遥感影像与基准影像的配准中。

附图说明

图1是本发明的待配准遥感影像与参考影像匹配方法流程示意图；

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

具体实施方式

以下结合附图来对本发明进行详细的描绘。然而应当理解，附图的提供仅为了更好地理解本发明，它们不应该理解成对本发明的限制。

如图1所示，本发明提供的基于全局几何约束的遥感影像与基准影像的匹配方法，包括以下步骤：

1、构建遥感影像成像几何模型，目前常用的行星遥感影像成像几何模型主要分为严格成像几何模型和通用成像几何模型。

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

以月球侦察轨道器(lunar reconnaissance orbiter,LRO)窄角相机(narrowangle camera,NAC)为例，详细说明本发明的严格成像几何模型和通用几何模型的构建过程。严格成像几何模型的构建一般包括内定向和外定向两个过程，而通用几何模型的构建则需要以构建的严格成像几何模型为基础，具体过程为：

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

1.1)LRO NAC内定向

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

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.2)LRO NAC外定向

1.2.1)建立共线方程

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

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

1.2.2)初始外方位元素的读取

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

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

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

X_s(t)＝a₀+a₁t+a₂t²+a₃t³

Y_s(t)＝b₀+b₁t+b₂t²+b₃t³

Z_s(t)＝c₀+c₁t+c₂t²+c₃t³

ω(t)＝e₀+e₁t+e₂t²+e₃t³

κ(t)＝f₀+f₁t+f₂t²+f₃t³ (3)

式中，X_s(t),Y_s(t),Z_s(t)表示t时刻摄影中心在星固坐标系中的坐标，即外方位线元素；ω(t),κ(t)表示t时刻焦平面在星固坐标系中的姿态角，即外方位角元素；a₀...f₃表示对应参数的多项式系数，这些系数可使用最小二乘法根据轨道测量数据求解。

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

2)建立LRO NAC有理函数模型

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

2.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为高程。

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

2、读取参考影像数据和待匹配遥感影像，并对参考影像数据和待配准遥感影像数据分别进行预处理。

本发明中的参考影像数据是指已经具备了一定投影的影像产品，包含数字高程模型(Digital Elevation Model,DEM)和数字正射影像(Digital Orthophoto Map,DOM)。当参考影像数据分别只有DEM或DOM时，本发明采用不同的预处理方法进行处理，具体过程为：

2.1)参考影像数据只有DEM时，对DEM进行预处理，构建基准影像，具体过程为：

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

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

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

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

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

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

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

2.2)当参考影像数据为DOM时，DOM即为基准影像。

2.3)当参考影像数据为DOM时，将DOM作为参考将待配准遥感影像进行预处理，与DOM通过软件将进行直方图匹配处理，具体过程为：

(1)计算待配准遥感影像的直方图p_r(r_k),即遥感影像中灰度级r_k出现的概率：

a)将待配准遥感影像的灰度级归一化。若遥感图像的灰度级为0,1，…，L-1，则令：

式中，L为灰度级层次数，r_k为遥感影像第k级的灰度，k为遥感影像的灰度级别。

b)计算各灰度级的出现的概率p_r(r_k),与r_k相对的p_r(r_k)图形即为直方图。

假设n_k为灰度级为r_k的像素个数，N为总的像素个数，令：

式中，n_k为像素值为k的像素的频数，p_r(r_k)为其出现的概率，N为遥感影像的像素总数。

(2)将r_k根据变换(映射)函数T(r_k)进行变换，计算图像经过T(r_k)变换后得到的灰度级s_k，将s_k四舍五入为范围[0,L-1]内的整数，从而将遥感影像中灰度级为r_k的各像素映射到经过T(r_k)变换的输出影像中灰度级为s_k的对应像素得到。

式中，n_j为灰度值为r_j的像素数。

(3)将步骤(1)和(2)计算出的s_k与G(z_q)进行匹配，具体过程为：

A)采用步骤(1)的方法计算基准影像的灰度级直方图p_z(z_q)，并采用公式(13)对q＝0,1,2,…,L-1计算变换函数G的所有值，把G的值四舍五入为范围[0,L-1]内的整数，将G存储在一个表中。

式中，z_q为基准影像直方图匹配后第q级的灰度，z_i为基准影像直方图匹配后第i级的灰度。

B)对每一个值s_k,k＝0,1,2,…,L-1,使用步骤A)存储的G值寻找相应的z_q值，以使G(z_q)最接近s_k，并存储这些从s到z的映射，当满足给定s_k的z_q值多于一个时，按惯例选择最小的值。

C)使用步骤B)找到的映射将待配准遥感影像的每个变换后的像素值s_k映射为基准影像的相应z_q值，形成直方图匹配后的待配准遥感影像。

3、将基准影像与遥感影像进行仿射-尺度不变特征变换匹配，提取特征点，获取最小二乘匹配初值，具体过程为：

3.1)将基准影像与待处理的遥感影像进行初始匹配

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

3.2)计算有理函数模型精化参数

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

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

3.3)在基准影像上提取特征点，本发明以Harris特征点为例，但不仅限于此，具体过程为：

1)确定一个n×n大小的影像窗口，对窗口内的每一个像素点进行一阶差分运算，求得在x，y方向的梯度g_x，g_y。

2)对梯度值进行高斯滤波：

式中，Ga为高斯卷积模板，σ取0.4，g_x为灰度在x方向的梯度，g_y为y方向的梯度。

3)计算矩阵M及兴趣值I：

式中，det是矩阵的行列式，tr是矩阵的迹，k是默认常数，取0.04。

4)选取兴趣值的局部极值点，并在设定窗口内取最大值作为局部极值点。

局部极值点的数目可根据特征点数提取的数目要求，对所有的极值点进行降序排序，根据特征点提取的数目(num)要求(在此取特征点数目num＝500个)，选出降序排序的兴趣值中前num个点作为最后的结果，确定一个n×n大小的影像窗口。

3.4)计算全局最小二乘匹配初值

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

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

式中，Y与X分别为特征点的纬度与经度坐标,Y₀与X₀分别为基准影像左上角点的纬度与经度坐标，f_y，f_x为特征点在基准影像上的行列号，R_y，R_x为基准影像在y方向以及x方向的分辨率。根据计算得出的经纬度坐标，在数字高程模型上提取对应经纬度的高程值Z。

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

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

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

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

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

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

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

式中，为泰勒级数0次项，分别表示列方向与行方向坐标的一阶偏导数。

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

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

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

4.5)输出遥感影像的辐射变形参数以及几何变形参数，完成遥感影像与基准影像间的全局最小二乘匹配，通过解求得出的变形参数可以为大区域制图提供高精度控制点，与此同时，通过解求的参数可以将影像制作生成正射影像，为大区域制图提供基础。

上述各实施例仅用于说明本发明，其中方法的各实施步骤等都是可以有所变化的，凡是在本发明技术方案的基础上进行的等同变换和改进，均不应排除在本发明的保护范围之外。

Claims (8)

1.一种基于全局几何约束的遥感影像与基准影像的匹配方法，其特征在于包括以下内容：

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

读取参考影像数据和待匹配遥感影像，对参考影像数据预处理得到基准影像；

对待匹配遥感影像进行预处理；

将基准影像与待配准遥感影像进行仿射-尺度不变特征变换匹配，提取特征点，并获取最小二乘匹配的初值；

基于最小二乘匹配的初值通过全局约束最小二乘匹配模型进行最小二乘匹配，输出待配准遥感影像的辐射变形参数以及几何变形参数。

2.根据权利要求1所述的基于全局几何约束的遥感影像与基准影像的匹配方法，其特征在于，当参考影像数据只有DEM时，对DEM进行预处理构建基准影像，具体过程为：

2.1)在DEM上，逐像素计算在x与y方向的梯度；

2.2)根据梯度值，逐像素计算坡度S、坡向值A；

2.3)根据参考影像数据的头文件信息，创建与DEM相同大小的空影像，逐像素计算其对应的模拟辐射值信息并存入空影像，构建得到基准影像：

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

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

3.根据权利要求1所述的基于全局几何约束的遥感影像与基准影像的匹配方法，其特征在于，当参考影像数据为DOM时，DOM即为基准影像。

4.根据权利要求3所述的基于全局几何约束的遥感影像与基准影像的匹配方法，其特征在于，对待配准遥感影像进行预处理的具体过程为：

2.1)计算待配准遥感影像的直方图p_r(r_k)，具体过程为：

a)将待配准遥感影像的灰度级归一化：

式中，L为灰度级层次数，r_k为遥感影像第k级的灰度，k为遥感影像的灰度级别；

b)计算各灰度级出现的概率p_r(r_k)：

式中，n_k为像素值为k的像素的频数，p_r(r_k)为其出现的概率，N为遥感影像的像素总数；

2.2)将r_k根据变换函数T(r_k)进行变换得到的灰度级s_k，将s_k四舍五入为范围[0,L-1]内的整数：

式中，n_j为灰度值为r_j的像素数；

2.3)将步骤2.1)和2.2)计算出的s_k与G(z_q)进行匹配，具体过程为：

A)采用步骤2.1)计算基准影像的灰度级直方图p_z(z_q)，并对q＝0,1,2,…,L-1计算变换函数G的所有值，把G的值四舍五入为范围[0,L-1]内的整数，并将G存储在一个表中：

式中，z_q为基准影像直方图匹配后第q级的灰度，z_i为基准影像直方图匹配后第i级的灰度；

B)对每一个值s_k,k＝0,1,2,…,L-1,使用步骤A)存储的G值寻找相应的z_q值，以使G(z_q)最接近s_k，并存储这些从s到z的映射，当满足给定s_k的z_q值多于一个时，按惯例选择最小的值；

C)使用步骤B)找到的映射将待配准遥感影像的每个变换后的像素值s_k映射为基准影像的相应z_q值，形成直方图匹配后的待配准遥感影像。

5.根据权利要求4所述的基于全局几何约束的遥感影像与基准影像的匹配方法，其特征在于，计算待配准遥感影像所有灰度级s_k的具体过程为：

a)将待配准遥感影像的灰度级归一化，若遥感图像的灰度级为0,1，…，L-1，则令：

式中，L为灰度级层次数，r_k为遥感影像第k级的灰度，k为遥感影像的灰度级别；

b)计算各灰度级的像素频数

假设n_k为灰度级为r_k的像素个数，N为总的像素个数，令：

式中，n_k为像素值为k的像素的频数，p_r(r_k)为其出现的概率

c)计算第k级的灰度级

式中，n_j为灰度值为r_j的像素数。

6.根据权利要求1所述的基于全局几何约束的遥感影像与基准影像的匹配方法，其特征在于，将基准影像与待配准遥感影像进行仿射-尺度不变特征变换匹配，提取特征点，并获取最小二乘匹配的初值的具体过程为：

3.1)采用仿射-尺度不变特征变换算法获得待配准遥感影像与基准影像之间的初始匹配点；

3.2)根据初始匹配点结合DEM计算有理函数模型精化参数

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

3.3)在基准影像上提取特征点；

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

7.根据权利要求6所述的基于全局几何约束的遥感影像与基准影像的匹配方法，其特征在于，在基准影像上提取特征点的具体过程为：

1)确定一个n×n大小的影像窗口，对窗口内的每一个像素点进行一阶差分运算，求得在x，y方向的梯度g_x，g_y；

2)对梯度值进行高斯滤波：

式中，Ga为高斯卷积模板，σ取0.4，g_x为灰度在x方向的梯度，g_y为y方向的梯度；

3)计算矩阵M及兴趣值I：

I＝det(M)-k·tr²(M)

式中，det是矩阵的行列式，tr是矩阵的迹，k是默认常数，取0.04；

4)选取兴趣值的局部极值点，并在设定窗口内取最大值作为特征点。

8.根据权利要求1所述的基于全局几何约束的遥感影像与基准影像的匹配方法，其特征在于，基于最小二乘匹配初值通过全局约束最小二乘匹配模型进行最小二乘匹配，输出待配准遥感影像的辐射变形参数以及几何变形参数的具体过程为：

4.1)引入全局约束最小二乘匹配模型：

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

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

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

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

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

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

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

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

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

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

式中，为泰勒级数0次项，分别表示列方向与行方向坐标的一阶偏导数；

即针对每个特征点，误差方程为：

v＝Bx-L

其中，

x＝(Δa₀ Δa₁ Δa₂ Δb₀ Δb₁ Δb₂ Δh_0i Δh_1i Δδ_ri Δδ_ci)^T

式中，m为连接点的点对数，n为窗口大小；

4.3)对每个输入的匹配初值点代入线性化的误差方程，逐点计算法方程；

4.4)根据最小二乘原理通过对x＝(B^TB)^-1B^TL进行求解可以计算出所有未知数，并根据相关系数计算公式计算相关系数ρ²：

判断平差后计算得出的相关系数，如所有点相关系数的和满足平差精度要求则输出平差结果，平差精度不满足要求，则重复步骤4.2)～4.4)，迭代求解；

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