CN103473782B - 基于物方垂直双面元的最小二乘匹配方法 - Google Patents

Abstract

本发明是一种基于物方垂直双面元的最小二乘匹配方法，主要用于城市地区遥感影像中建筑棱角线上特征点的匹配包括：首先，在密集匹配的点云中，建立每个点的初始法向量，包括两条垂直的初始法向量；然后基于各三维点的初始法向量，建立误差方程，进行双面元参数的迭代优化，搜索最佳匹配点位。本发明相比原有的单面元算法，能够匹配出建筑物棱角线上的特征点，大幅提高了遥感影像中建筑物的模型精度。

Description

基于物方垂直双面元的最小二乘匹配方法

技术领域

本发明属于遥感图像处理影像匹配技术领域，涉及一种基于物方垂直双面元的最小二乘匹配方法。

背景技术

影像匹配技术目前已经被广泛地应用于三维重建、拼接、检索、目标跟踪等方面，在民用、军事医疗领域具有重要的应价值。由于物体表面在大多数情况下并不是理想的。比如外界各种遮挡关系、光照变化的影响，即使是空间场景物体上的同一点在不同的二维图像上的对应投影点的亮度值或颜色值并不是完全相同；同样地，具有相同亮度值或颜色值的二维图像投影点也未必对应物体表面的同一点。影像匹配技术是计算机视觉及摄影测量领域中最基础、也是最困难的工作之一。

随着城市化速度不断加快，建立精确的城市三维模型，实现由数字城市到智慧城市的过度，已成为3S技术发展的主要趋势。其中，高精度的密集匹配，是建立精确城市三维模型的基础。城市建筑最大的特点是存在大量的棱角线、点，而对于这些建筑物棱角线上的点，由于其在大部分匹配算法中会被模糊处理，如何匹配出建筑物棱角上的点，一直是密集匹配与三维重建中的难点。

最小二乘匹配是Ackerman教授于1983年提出的一种高精度的影像匹配方法，由于其具有灵活、可靠、高精度的特点，很快推广到了多点和多片的情况，包括基于SIFT的宽基线立体影像最小二乘匹配以及共线条件约束的多片最小二乘匹配算法MPGC等。为了改进MPGC算法中多幅待匹配影像间存在独立变形参数这一不足，江万寿教授提出了基于物方面元的最小二乘匹配，该算法利用共线条件方程作为几何约束关系，利用物方面元参数以及单中心投影的投影方程消除参数与匹配的各片之间独立的变形系数，减少未知数的个数，同时得到空间面元的法向量。但这些方法在优化的过程中会模糊包括建筑物顶部边缘在内的垂直边缘点，使得优化后的点云无法得到清晰的棱角线。具体来说，传统的单面元最小二乘匹配方法中，斜面元的几何意义是三维点的切平面。将整个点云的三维模型看作一个连续曲面，设其曲面方程为F(x,y,z)=0(即z=f(x,y))，面元的法向量（同样也是曲面的法向量）可表示为 （也就是。但在城市地区遥感影像中存在大量的建筑物，使得地面高层起伏大，并不是一个连续的曲面，此时城市建筑物中存在大量的棱角线，面元的法向量不止一个，传统的单面元最小二乘匹配方法会模糊这样的棱角线。尤其当匹配的窗口变大，这种模糊现象就越严重。

发明内容

针对现有技术的缺点，本发明的目的是提供一种能够突出建筑物棱角线上的特征点，提高遥感影像中建筑物的模型精度的匹配方法。

本发明的技术方案为一种基于物方垂直双面元的最小二乘匹配方法，包括以下步骤：

步骤1，在密集匹配的三维点云数据中，建立每个三维点的初始法向量，

包括以每个三维点为初始法向量的起点，以(α,β)、(π/2-α,π/2-β)为初始方向角，建立两条垂直的初始法向量，

a₁＝cosαcosβ

b₁＝sinαcosβ

c₁＝sinβ(1)

a₂＝cos(π/2-α)cos(π/2-β)＝sinαsinβ

b₂＝sin(π/2-α)cos(π/2-β)＝cosαsinβ

c₂＝sin(π/2-β)＝cosβ

其中，(a₁,b₁,c₁)、(a₂,b₂,c₂)表示两条垂直的初始法向量的坐标；

步骤2，基于各三维点的初始法向量，建立误差方程，进行双面元参数的迭代优化，搜索最佳匹配点位，实现包括以下步骤，

步骤2.1，将三维点(X_c,Y_c,Z_c)投影到参考影像上，得到参考影像上相应像点坐标(x_0c,y_0c)；

步骤2.2，以(x_0c,y_0c)为中心，一个像素为采样间隔，从参考影像上取出一个μ×μ个像素大小的影像窗口，计算窗口内每个像素的像点坐标(x₀,y₀)，其中μ为预设值；

步骤2.3，为三维点在物方建立一对互相垂直的面元P1、P2，面元中心为三维点的物方坐标值(X_c,Y_c,Z_c)，P1、P2的法向量为(a₁,b₁,c₁)、(a₂,b₂,c₂)；

步骤2.4，分别计算步骤2.2所取影像窗口中的μ×μ个点投影到P1、P2上每个投影点的物方坐标(X₁,Y₁,Z₁)、(X₂,Y₂,Z₂)；

步骤2.5，分别计算P1、P2上每个投影点投影到各搜索影像上的像点坐标(x_i,y_i)，得到搜索影像上的相应影像窗口，其中投影到第i张搜索影像上的像点坐标记为(x_i,y_i)；

步骤2.6，分别对参考影像和搜索影像上影像窗口中的点采用双线性内插法进行重采样，得到参考影像上影像窗口中的点重采样后的像素亮度值g₀(x₀,y₀)和P1上的投影点投影到第i张搜索影像上的像点重采样后的像素亮度值g_i(x_i1,y_i1)、P2上的投影点投影到第i张搜索影像上的像点重采样后的像素亮度值g_i(x_i2,y_i2)；

步骤2.7，根据共线方程约束和最小二乘匹配原理建立误差方程如下，

$v_{P1}={\mathrm{dh}}_{0i}+g_i(x_{i1},y_{i1})\·{\mathrm{dh}}_{1i}+h_{1i}(\frac{{\∂g}_{i1}}{{\∂x}_{i1}}{\mathrm{dx}}_{i1}+\frac{{\∂g}_{i1}}{{\∂y}_{i1}}{\mathrm{dy}}_{i1})-(g_0(x_0,y_0)-h_{0i}-h_{1i}\·g_i(x_{i1},y_{i1}))---\left(15\right)$

$v_{P2}={\mathrm{dh}}_{0i}+g_i(x_{i2},y_{i2})\·{\mathrm{dh}}_{1i}+h_{1i}(\frac{{\∂g}_{i2}}{{\∂x}_{i2}}{\mathrm{dx}}_{i2}+\frac{{\∂g}_{i2}}{{\∂y}_{i2}}{\mathrm{dy}}_{i2})-(g_0(x_0,y_0)-h_{0i}-h_{1i}\·g_i(x_{i2},y_{i2}))$

其中，v_P1、v_P2为投影误差；h_0i、h_1i为第i张搜索影像到参考影像的辐射畸变系数；

步骤2.8，进行辐射畸变改正，解求得到未知数的改正值，包括辐射畸变改正值dh_0i、dh_1i，面元角度改正值2个dα、dβ和物点坐标改正值dZ_c；

步骤2.9，当所有改正值小于预设改正值阈值时，停止迭代，取搜索影像窗口中心点为最佳匹配点位；否则，根据改正值得到改正后的参数h_0i、h_1i、Z_c、α、β，基于改正后的参数α、β根据公式(1)计算改正后的面元法向量作为新的(a₁,b₁,c₁)、(a₂,b₂,c₂)，重复步骤2.3至2.8，并将迭代次数计算值+1，当迭代次数大于预设次数阈值仍未出现所有改正值小于预设改正值阈值时，判断迭代失败，没有找到该三维点的匹配点。

而且，步骤2.8中解求未知数的改正值方式为，用差分代替微分来计算亮度值在x、y方向的导数值，将解求的投影系数带入光束方程，联立共线方程解求搜索影像中像点坐标的改正值；简化未知数，求解改正值。

而且，步骤2.9中，停止迭代后根据下式计算影像间的相关系数r_i1、r_i2，

$r_{i1}=\frac{\mathrm{\Σ}(g_0(x_{0r},y_{0s})-\stackrel{\‾}{g_0(x_0,y_0)})\left((g_{i1}(x_{i1r},y_{i1s})-\stackrel{\‾}{g_{i1}(x_{i1},y_{i1})})\right)}{\sqrt{\mathrm{\Σ}{(g_0(x_{0r},y_{0r})-\stackrel{\‾}{g_0(x_0,y_0)})}^2\·\mathrm{\Σ}{(g_{i1}(x_{i1r},y_{i1s})-\stackrel{\‾}{g_{i1}(x_{i1},y_{i1})})}^2}}$

$r_{i2}=\frac{\mathrm{\Σ}(g_0(x_{0r},y_{0s})-\stackrel{\‾}{g_0(x_0,y_0)})\left((g_{i2}(x_{i2r},y_{i2s})-\stackrel{\‾}{g_{i2}(x_{i2},y_{i2})})\right)}{\sqrt{\mathrm{\Σ}{(g_0(x_{0r},y_{0s})-\stackrel{\‾}{g_0(x_0,y_0)})}^2\·\mathrm{\Σ}{(g_{i2}(x_{i2r},y_{i2s})-\stackrel{\‾}{g_{i2}(x_{i2},y_{i2})})}^2}}---\left(21\right)$

其中，r_i1、r_i2分别为参考影像上影像窗口与从面元P1、P2投影到第i张搜索影像所得影像窗口间的相关系数，r、s为像素在影像窗口中的行列号；g₀(x_0r,y_0s)为参考影像上影像窗口中的点重采样后的像素亮度值，其中(x_0r,y_0s)表示参考影像中影像窗口内的第r行第s列像点坐标；g_i1(x_i1r,y_i1s)、g_i2(x_i2r,y_i2s)分别为面元P1、P2上的投影点投影到第i张搜索影像上的像点重采样后的像素亮度值，其中(x_i1r,y_i1s)、(x_i2r,y_i2s)分别为从面元P1、P2投影到第i张搜索影像所得影像窗口内第r行第s列的像点坐标； 分别为g₀(x_0r,y_0s)、g_i1(x_i1r,y_i1s)和g_i2(x_i2r,y_i2s)的均值。

而且，第k次迭代执行步骤2.9时，根据改正值得到改正后的参数h_0i、h_1i、Z_c、α、β如下式，

h_0i ^k＝h_0i ^k-1+dh_0i ^k+h_0i ^k-1dh_1i ^k(22)

h_1i ^k＝h_1i ^k-1+h_1i ^k-1dh_1i ^k

Z_c ^k＝Z_c ^k-1+dZ_c ^k

α^k＝α^k-1+dα^k(23)

β^k＝β^k-1+dβ^k

其中，h_0i ^k、h_1i ^k表示第k次迭代求得的参数h_0i、h_1i的值，dh_0i ^k、dh_1i ^k表示第k次求得的参数h_0i、h_1i的改正值；Z_c ^k、α^k、β^k分别表示第k次迭代求得的参数Z_c、α、β的值，dZ_c ^k、dα^k、dβ^k分别表示第k次求得的参数Z_c、α^k、β^k的改正值。

本发明的技术方案在物方逐个三维点建立两个互相垂直的面元，根据参考影像、面元、搜索影像的投影关系建立几何约束条件，以最小二乘思想建立误差方程，通过迭代求得匹配结果。本发明中以两个互相垂直的面元作为三维点的切平面克服了棱角线上点切平面的法向量不唯一的问题，对其具有很好的优化效果，实现突出建筑物棱角线，提高建筑物密集匹配精度。

附图说明

图1为本发明实施例的流程图。

图2为本发明实施例中建立初始法向量的示意图。

图3为本发明实施例中相关参数的几何关系图。

具体实施方式

为了更好地理解本发明的技术方案，下面结合附图对本发明做进一步的详细说明。

本发明技术方案可采用计算机软件技术实现自动运行流程。本发明的实施例是对两幅城市区域无人机影像进行匹配，参照图1，本发明实施例流程的步骤如下：

步骤1，建立初始法向量。

实施例首先输入影像相关参数，采用SGM或其他的基于图像的密集点云生成算法生成密集的三维点云数据，输出的数据包括三维点云数据中任一三维点的物方坐标值P(X_c,Y_c,Z_c)，然后参照图2为三维点云数据中每个三维点建立两条互相垂直的初始法向量（图中S为摄影中心；p(u_0c,v_0c,w_0c)为像点在像空间辅助坐标系中的坐标）。在密集匹配的点云中，逐个点建立初始法向量后，每个三维点位于两条垂直法向量的交点。通过每个三维点两条垂直法向量的交点，可以搜索出位于建筑棱角线附近的特征点。

实施例中，对每个三维点建立初始法向量的具体规则为：

以三维点为初始法向量的起点，设两条垂直法向量的初始方向角为(α,β)、(π/2-α,π/2-β)，由法向量与方向角的关系得：

a₁＝cosαcosβ

b₁＝sinαcosβ

c₁＝sinβ(1)

a₂＝cos(π/2-α)cos(π/2-β)＝sinαsinβ

b₂＝sin(π/2-α)cos(π/2-β)＝cosαsinβ

c₂＝sin(π/2-β)＝cosβ

(a₁,b₁,c₁)、(a₂,b₂,c₂)表示为三维点建立的两条垂直的法向量坐标。

步骤2，参照图3的几何关系，采用共线方程约束建立误差方程，逐个三维点进行双面元参数的迭代优化。实施例设阈值为10^-5，即当求解的每个参数的改正值均小于10^-5时停止迭代，否则，采用求得的改正值改正相关参数继续迭代。记录迭代次数，当迭代次数超过1000次时，迭代失败，该点无效。

实施例的步骤b对每个三维点分别进行迭代优化的具体实现步骤如下：

步骤2.1，将三维点(X_c,Y_c,Z_c)投影到参考影像上，按式(2)-(5)计算参考影像上相应像点坐标(x_0c,y_0c)：

$\left(\begin{array}{c}{\mathrm{ku}}_{0c}\\ {\mathrm{kv}}_{0c}\\ {\mathrm{kw}}_{0c}\end{array}\right)=R^{-1}\×\left(\begin{array}{c}{X}_c-X_{s0}\\ Y_c-Y_{s0}\\ Z_c-Z_{s0}\end{array}\right)---\left(2\right)$

$\mathrm{\λ}=R^{-1}\left(3\right)\×{\left(\begin{array}{ccc}{\mathrm{ku}}_{0c}& {\mathrm{kv}}_{0c}& {\mathrm{kw}}_{0c}\end{array}\right)}^T/(-f_0)---\left(3\right)$

$\left(\begin{array}{c}{u}_{0c}\\ v_{0c}\\ w_{0c}\end{array}\right)=\left(\begin{array}{c}{\mathrm{ku}}_{0c}\\ {\mathrm{kv}}_{0c}\\ {\mathrm{kw}}_{0c}\end{array}\right)/\mathrm{\λ}---\left(4\right)$

$\left(\begin{array}{c}{x}_{0c}\\ y_{0c}\end{array}\right)=\left(\begin{array}{c}{R}^{-1}\left(1\right)\\ R^{-1}\left(2\right)\end{array}\right)\×\left(\begin{array}{c}{u}_{0c}\\ v_{0c}\\ w_{0c}\end{array}\right)---\left(5\right)$

式(2)中，R为参考影像的旋转矩阵，(X_c,Y_c,Z_c)为三维点的物方坐标值，(X_s0,Y_s0,Z_s0)为参考影像外方位元素中的线元素，(ku_0c,kv_0c,kw_0c)为中间变量；

式(3)中，R^-1(3)为参考影像旋转矩阵逆矩阵的第三行，f₀为参考影像主距，λ为三维点投影到参考影像的投影系数；

式(4)中，(u_0c,v_0c,w_0c)为像点(x_0c,y_0c)的像空间辅助坐标；

式(5)中，(x_0c,y_0c)为三维点(X_c,Y_c,Z_c)投影到参考影像上的像点坐标，R^-1(1)、R^-1(2)为参考影像旋转矩阵逆矩阵的第一、二行。

步骤2.2，以(x_0c,y_0c)为中心，一个像素为采样间隔，从参考影像上取出一个μ×μ个像素（μ可由用户预设，一般μ为7、9或11）大小的影像窗口，计算窗口内每个像素的像点坐标(x₀,y₀)，实施例以毫米为单位进行计算。

步骤2.3，为三维点在物方建立一对互相垂直的面元P1、P2，面元中心为三维点的物方坐标值(X_c,Y_c,Z_c)，第一次执行步骤2.3时P1、P2法向量分别为步骤1.2中求得的(a₁,b₁,c₁)、(a₂,b₂,c₂)，之后执行步骤2.3时，采用上一次迭代执行步骤2.9时计算改正后的面元法向量所得新的(a₁,b₁,c₁)、(a₂,b₂,c₂)。

按公式(6)计算垂直面元的方程：

P1：a₁(X-X_c)+b₁(Y-Y_c)+c₁(Z-Z_c)＝0(6)

P2：a₂(X-X_c)+b₂(Y-Y_c)+c₂(Z-Z_c)＝0

式(6)中，(X,Y,Z)表示面元P1、P2上的点。

步骤2.4，按公式(7)-(9)，分别计算步骤2.2所取影像窗口中的μ×μ个点投影到P1、P2上每个投影点的物方坐标(X₁,Y₁,Z₁)、(X₂,Y₂,Z₂)，为简化公式表达起见，实施例采用1(2)作为下标表示P1或P2的情况：

$\left(\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right)=R\×\left(\begin{array}{c}{x}_0\\ y_0\\ -f_0\end{array}\right)---\left(7\right)$

${\mathrm{\λ}}_{1\left(2\right)}=\frac{a_{1\left(2\right)}(X_c-X_{s0})+b_{1\left(2\right)}(Y_c-Y_{s0})+c_{1\left(2\right)}(Z_c-Z_{s0})}{a_{1\left(2\right)}{u}_0+b_{1\left(2\right)}{v}_0+c_{1\left(2\right)}{w}_0}---\left(8\right)$

$\left(\begin{array}{c}{X}_{1\left(2\right)}\\ Y_{1\left(2\right)}\\ Z_{1\left(2\right)}\end{array}\right)=\left(\begin{array}{c}{X}_{s0}\\ Y_{s0}\\ Z_{s0}\end{array}\right)+{\mathrm{\λ}}_{1\left(2\right)}\left(\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right)---\left(9\right)$

式(7)中，(u₀,v₀,w₀)为影像窗口中各点的像空间辅助坐标；

式(8)中，λ₁₍₂₎表示λ₁或λ₂，分别为影像窗口投影到P1、P2上的投影系数；a₁₍₂₎、b₁₍₂₎、c₁₍₂₎，即a₁或a₂,b₁或b₂,c₁或c₂；

式(9)中，(X_s0,Y_s0,Z_s0)为参考影像外方位元素中的线元素，(X₁₍₂₎,Y₁₍₂₎,Z₁₍₂₎)表示影像窗口投影到P1、P2上的物方坐标(X₁,Y₁,Z₁)或(X₂,Y₂,Z₂)。

步骤2.5，按公式(10)-(13)，分别计算P1、P2上每个投影点投影到各搜索影像上的像点坐标(x_i,y_i)，得到搜索影像上的相应影像窗口。其中从P1、P2分别投影到搜索影像上各形成一个影像窗口。相应计算公式如下：

$\left(\begin{array}{c}{\mathrm{ku}}_{i1\left(2\right)}\\ {\mathrm{kv}}_{i1\left(2\right)}\\ {\mathrm{kw}}_{i1\left(2\right)}\end{array}\right)={R_i}^{-1}\×\left(\begin{array}{c}{X}_{1\left(2\right)}-X_{\mathrm{si}}\\ Y_{1\left(2\right)}-Y_{\mathrm{si}}\\ Z_{1\left(2\right)}-Z_{\mathrm{si}}\end{array}\right)---\left(10\right)$

${\mathrm{\λ}}_{i1\left(2\right)}={R_i}^{-1}\left(3\right)\×{\left({\mathrm{ku}}_{i1\left(2\right)}{\mathrm{kv}}_{i1\left(2\right)}{\mathrm{kw}}_{i1\left(2\right)}\right)}^T/(-f_i)---\left(11\right)$

$\left(\begin{array}{c}{u}_{i1\left(2\right)}\\ v_{i1\left(2\right)}\\ w_{i1\left(2\right)}\end{array}\right)=\left(\begin{array}{c}{\mathrm{ku}}_{i1\left(2\right)}\\ {\mathrm{kv}}_{i1\left(2\right)}\\ {\mathrm{kw}}_{i1\left(2\right)}\end{array}\right)/{\mathrm{\λ}}_{i1\left(2\right)}---\left(12\right)$

$\left(\begin{array}{c}{x}_{i1\left(2\right)}\\ y_{i1\left(2\right)}\end{array}\right)=\left(\begin{array}{c}{R_i}^{-1}\left(1\right)\\ {R_i}^{-1}\left(2\right)\end{array}\right)\×\left(\begin{array}{c}{u}_{i1\left(2\right)}\\ v_{i1\left(2\right)}\\ w_{i1\left(2\right)}\end{array}\right)---\left(13\right)$

式(10)中，R_i为第i张搜索影像的旋转矩阵，R_i ^-1为第i张搜索影像旋转矩阵的逆矩阵，(X_si,Y_si,Z_si)为第i张搜索影像外方位元素中的线元素，(ku_i1(2),kv_i1(2),kw_i1(2))为中间变量；

式(11)中，R_i ^-1(3)为第i张搜索影像旋转矩阵逆矩阵的第三行，f_i为第i张搜索影像主距，λ_i1(2)为三维点投影到第i张搜索影像的投影系数；

式(12)中，(u_i1(2),v_i1(2),w_i1(2))分别为P1、P2上的投影点投影到第i张搜索影像上像点(x_i1(2),y_i1(2))的像空间辅助坐标；

式(13)中，(x_i1(2),y_i1(2))分别为第i张搜索影像上的像点坐标，R_i ^-1(1)、R_i ^-1(2)为第i张搜索影像旋转矩阵逆矩阵的第一、二行。

i的取值取决于搜索影像的数量。

步骤2.6，分别对参考影像和搜索影像上影像窗口中的点采用双线性内插法进行重采样，得到g₀(x₀,y₀)和g_i(x_i1(2),y_i1(2))；

步骤2.7，按共线方程约束（公式(14)）和最小二乘匹配原理建立误差方程（公式(15)）；

$x_{i1\left(2\right)}=-f_i\frac{a_{1i}(X_{1\left(2\right)}-X_{\mathrm{si}})+b_{1i}(Y_{1\left(2\right)}-Y_{\mathrm{si}})+c_{1i}(Z_{1\left(2\right)}-Z_{\mathrm{si}})}{a_{3i}(X_{1\left(2\right)}-X_{\mathrm{si}})+b_{3i}(Y_{1\left(2\right)}-Y_{\mathrm{si}})+c_{3i}(Z_{1\left(2\right)}-Z_{\mathrm{si}})}$

$y_{i1\left(2\right)}=-f_i\frac{a_{2i}(X_{1\left(2\right)}-x_{\mathrm{si}})+b_{2i}(Y_{1\left(2\right)}-Y_{\mathrm{si}})+c_{2i}(Z_{1\left(2\right)}-Z_{\mathrm{si}})}{a_{3i}(X_{1\left(2\right)}-X_{\mathrm{si}})+b_{3i}(Y_{1\left(2\right)}-Y_{\mathrm{si}})+c_{3i}(Z_{1\left(2\right)}-Z_{\mathrm{si}})}---\left(14\right)$

其中： $\left(\begin{array}{ccc}{a}_{1i}& a_{2i}& a_{3i}\\ b_{1i}& b_{2i}& b_{3i}\\ c_{1i}& c_{2i}& c_{3i}\end{array}\right)=R_i$

$v_{P1}={\mathrm{dh}}_{0i}+g_i(x_{i1},y_{i1})\·{\mathrm{dh}}_{1i}+h_{1i}(\frac{{\∂g}_{i1}}{{\∂x}_{i1}}{\mathrm{dx}}_{i1}+\frac{{\∂g}_{i1}}{{\∂y}_{i1}}{\mathrm{dy}}_{i1})-(g_0(x_0,y_0)-h_{0i}-h_{1i}\·g_i(x_{i1},y_{i1}))---\left(15\right)$

$v_{P2}={\mathrm{dh}}_{0i}+g_i(x_{i2},y_{i2})\·{\mathrm{dh}}_{1i}+h_{1i}(\frac{{\∂g}_{i2}}{{\∂x}_{i2}}{\mathrm{dx}}_{i2}+\frac{{\∂g}_{i2}}{{\∂y}_{i2}}{\mathrm{dy}}_{i2})-(g_0(x_0,y_0)-h_{0i}-h_{1i}\·g_i(x_{i2},y_{i2}))$

式(14)中， $R_i=\left(\begin{array}{ccc}{a}_{1i}& a_{2i}& a_{3i}\\ b_{1i}& b_{2i}& b_{3i}\\ c_{1i}& c_{2i}& c_{3i}\end{array}\right)$ 为搜索影像的旋转矩阵，a_1i、a_2i、a_3i、b_1i、b_2i、b_3i、c_1i、c_2i、c_3i为其中的参数；

式(15)中，v_P1、v_P2为投影误差；h_0i、h_1i为第i张搜索影像到参考影像的辐射畸变系数；g_i(x_i1(2),y_i1(2))分别为P1、P2上的投影点投影到第i张搜索影像上的像点重采样后的像素亮度值；g₀(x₀,y₀)为参考影像上影像窗口中的点重采样后的像素亮度值。

步骤2.8，辐射畸变改正。解求未知数的改正值方法如下：

用差分代替微分来计算亮度值在x、y方向的导数值；

将解求的投影系数带入光束方程，联立共线方程解求搜索影像中像点坐标的改正值；

简化未知数，求解改正值。

当迭代所求得的改正值均小于阈值时，所求得的搜索影像的像点坐标即为参考影像的匹配点，将相关数据记录到txt文档，当所有三维点的匹配点求解完毕，得到一个完整的三维模型。

实施例的具体实现如下：

由于影像是不连续的且窗口中每个像素的采样间隔是相同的，可用差分代替微分来计算亮度值在x、y方向的导数值 即公式(16)：

$\frac{{\∂g}_{i1\left(2\right)}}{{\∂x}_{i1}}=\frac{1}{2}(g_i(x_{i1\left(2\right)}+1,y_{i1})-g_{i1}(x_{i1\left(2\right)}-1,y_{i1}))---\left(16\right)$

$\frac{{\∂g}_{i1\left(2\right)}}{{\∂y}_{i1\left(2\right)}}=\frac{1}{2}(g_i(x_{i1\left(2\right)},y_{i1}+1)-g_i(x_{i1\left(2\right)},y_{i1}-1))$

注：式(16)中的坐标为像素坐标。

由于参考影像面元中心(u_0c,v_0c,w_0c)（即像点(x_0c,y_0c)的像空间辅助坐标）也在光束方程(17)上，将参考影像面元中心带入光束方程中(18)，由求得的中间变量λ_c，求解投影系数λ₁₍₂₎：

X＝X_s0+λu₀

Y＝Y_s0+λv₀(17)

Z＝Z_s0+λw₀

X_c＝X_s0+λ_cu_0c

Y_c＝Y_s0+λ_cv_0c(18)

Z_c＝Z_s0+λ_cw_0c

λ_c＝(Z_c-Z_s0)/w_0c

${\mathrm{\λ}}_{1\left(2\right)}=\frac{a_{1\left(2\right)}{u}_{0c}+b_{1\left(2\right)}{v}_{0c}+c_{1\left(2\right)}{w}_{0c}}{a_{1\left(2\right)}{u}_0+b_{1\left(2\right)}{v}_0+c_{1\left(2\right)}{w}_0}\·{\mathrm{\λ}}_c=\frac{a_{1\left(2\right)}{u}_{0c}+b_{1\left(2\right)}{v}_{0c}+c_{1\left(2\right)}{w}_{0c}}{a_{1\left(2\right)}{u}_0+b_{1\left(2\right)}{v}_0+c_{1\left(2\right)}{w}_0}\·\frac{(Z_c-Z_{s0})}{w_{0c}}---\left(19\right)$

将公式(17)、(19)带入公式(14)计算(x_i1(2),y_i1(2))，由公式(1)，(a₁,b₁,c₁)、(a₂,b₂,c₂)可以用参数α、β来表示，因此，误差方程(15)中的(dx_i1,dy_i1)、(dx_i2,dy_i2)可以按公式(20)计算：

${\mathrm{dx}}_{i1}=\frac{{\∂x}_{i1}}{{\∂Z}_c}{\mathrm{dZ}}_c+\frac{{\∂x}_{i1}}{\∂\mathrm{\α}}\mathrm{d\α}+\frac{{\∂x}_{i1}}{\∂\mathrm{\β}}\mathrm{d\β}$

${\mathrm{dy}}_{i1}=\frac{{\∂y}_{i1}}{{\∂Z}_c}d{Z}_c+\frac{{\∂y}_{i1}}{\∂\mathrm{\α}}\mathrm{d\α}+\frac{{\∂y}_{i1}}{\∂\mathrm{\β}}\mathrm{d\β}---\left(20\right)$

${\mathrm{dx}}_{i2}=\frac{{\∂x}_{i2}}{{\∂Z}_c}{\mathrm{dZ}}_c+\frac{{\∂x}_{i2}}{\∂\mathrm{\α}}\mathrm{d\α}+\frac{{\∂x}_{i2}}{\∂\mathrm{\β}}\mathrm{d\β}$

${\mathrm{dy}}_{i2}=\frac{{\∂y}_{i2}}{{\∂Z}_c}{\mathrm{dZ}}_c+\frac{{\∂y}_{i2}}{\∂\mathrm{\α}}\mathrm{d\α}+\frac{{\∂y}_{i2}}{\∂\mathrm{\β}}\mathrm{d\β}$

其中，i用于标识第i张搜索影像。

此时，误差方程(15)中的未知数已经化简为2n+3个：辐射畸变改正值2n个，即dh_0i、dh_1i；面元角度改正值2个，即dα、dβ；物点坐标改正值1个，即dZ_c，n为搜索影像的个数。误差方程(15)中的未知数化简后与单面元的未知数个数相同。

步骤2.9，具体实施时可根据精度要求预设改正值阈值和次数阈值，实施例中设改正值阈值为10^-5，即当所有改正值小于10^-5时，停止迭代，取第i张搜索影像窗口中心点为该搜索影像的最佳匹配点位。

否则，以所求改正值按公式(22)、(23)改正参数h_0i、h_1i、Z_c、α、β，基于改正后的参数α、β根据公式(1)计算改正后的面元法向量作为新的(a₁,b₁,c₁)、(a₂,b₂,c₂)，重复步骤2.3-2.8，并将迭代次数计算值+1（可设迭代次数的初始值为0），当迭代次数大于1000仍未出现任一搜索影像所有改正值小于10^-5时，认为迭代失败，该三维点无效，即没有找到该点的匹配点。

为便于评价匹配精度，本发明进一步提出可在迭代结束后按公式(21)计算影像间的相关系数r_i1、r_i2：

$r_{i1}=\frac{\mathrm{\Σ}(g_0(x_{0r},y_{0s})-\stackrel{\‾}{g_0(x_0,y_0)})\left((g_{i1}(x_{i1r},y_{i1s})-\stackrel{\‾}{g_{i1}(x_{i1},y_{i1})})\right)}{\sqrt{\mathrm{\Σ}{(g_0(x_{0r},y_{0r})-\stackrel{\‾}{g_0(x_0,y_0)})}^2\·\mathrm{\Σ}{(g_{i1}(x_{i1r},y_{i1s})-\stackrel{\‾}{g_{i1}(x_{i1},y_{i1})})}^2}}$

$r_{i2}=\frac{\mathrm{\Σ}(g_0(x_{0r},y_{0s})-\stackrel{\‾}{g_0(x_0,y_0)})\left((g_{i2}(x_{i2r},y_{i2s})-\stackrel{\‾}{g_{i2}(x_{i2},y_{i2})})\right)}{\sqrt{\mathrm{\Σ}{(g_0(x_{0r},y_{0s})-\stackrel{\‾}{g_0(x_0,y_0)})}^2\·\mathrm{\Σ}{(g_{i2}(x_{i2r},y_{i2s})-\stackrel{\‾}{g_{i2}(x_{i2},y_{i2})})}^2}}---\left(21\right)$

其中，r_i1、r_i2分别为参考影像上影像窗口与从面元P1、P2投影到第i张搜索影像所得影像窗口间的相关系数，r、s为像素在影像窗口中的行列号；g₀(x_0r,y_0s)为参考影像上影像窗口中的点重采样后的像素亮度值，其中(x_0r,y_0s)表示参考影像中影像窗口内的第r行第s列像点坐标；g_i1(x_i1r,y_i1s)、g_i2(x_i2r,y_i2s)分别为面元P1、P2上的投影点投影到第i张搜索影像上的像点重采样后的像素亮度值，其中(x_i1r,y_i1s)、(x_i2r,y_i2s)分别为从面元P1、P2投影到第i张搜索影像所得影像窗口内第r行第s列的像点坐标； 分别为g₀(x_0r,y_0s)、g_i1(x_i1r,y_i1s)和g_i2(x_i2r,y_i2s)的均值。

参数改正可依据以下公式：

h_0i ^k＝h_0i ^k-1+dh_0i ^k+h_0i ^k-1dh_1i ^k(22)

h_1i ^k＝h_1i ^k-1+h_1i ^k-1dh_1i ^k

Z_c ^k＝Z_c ^k-¹+dZ_c ^k

α^k＝α^k-1+dα^k(23)

β^k＝β^k-1+dβ^k

式(22)中，h_0i ^k、h_1i ^k表示第k次迭代求得的参数h_0i、h_1i的值，dh_0i ^k、dh_1i ^k表示第k次求得的参数h_0i、h_1i的改正值；

式(23)中，Z_c ^k、α^k、β^k分别表示第k次迭代求得的参数Z_c、α、β的值，dZ_c ^k、dα^k、dβ^k分别表示第k次求得的参数Z_c、α^k、β^k的改正值。

第一次执行步骤2.9时，改正h_0i ¹、h_1i ¹、Z_c ¹、α¹、β¹所用h_0i ⁰、h_1i ⁰、Z_c ⁰、α⁰、β⁰即初始的h_0i、h_1i、Z_c、α、β。

将迭代完成后的结果记录为txt文档，文档内容可提供全面的结果，包括：初始三维点坐标、改正后三维点坐标、参考影像像点坐标、搜索影像像点坐标、初始面元法向量、改正后面元法向量、初始匹配点相关系数、改正后匹配点相关系数、迭代次数、是否迭代成功等。具体实施时，可根据投影公式计算改正后的三维点坐标(X_c,Y_c,Z_c)。

以下通过仿真实验来验证本发明的有效性：

仿真实验采用两组航向/旁向重叠的城市区域无人机影像，并得到准确的内外方位元素，相邻影像间的重叠度均大于80%，标准影像大小为8206×6078。采用SGM密集匹配方法生成三维点云数据。

评价指标：对三维模型中建筑物棱角附近点的匹配成功率和匹配后影像窗口间相关系数两个方面进行评价。

（1）三维点的匹配成功率：统计完整的三维模型中的点云匹配成功的概率与单面元最小二乘匹配中点云匹配进行比较，比较迭代成功的概率。

比较了完整的三维模型中点云匹配成功的概率后，三维模型中的点云数据分为建筑物棱角点和非棱角点两类，分别统计匹配成功的概率，并与单面元最小二乘匹配进行比较。

（2）窗口间相关系数：影像窗口间的相关系数是评价两个影像窗口是否相关的重要评价指标，影像窗口相关系数越高，说明影像中心点是同名点的概率越高。与评价（1）相类似，通过整体和分类两方面来比较本发明与单面元最小二乘匹配法方法的匹配效果。

仿真结果：根据评价指标，仿真实验的实验数据如下表所示：

表1：实验结果统计表

从表1的实验结果可以看出，由于实验模型采取的垂直的双面元，对于垂直边缘（实验中的房顶边缘）具有很好的突出效果，同时由于本发明采取的几何约束和误差方程的建立方式与单面元最小二乘匹配相同，因此并没有损失非边缘区域的匹配精度。

综上所述，本发明具有如下的优点：

（1）采用垂直双面元代替传统的单面元最小二乘匹配构造几何模型，突出了垂直边缘。

（2）两张面元有固定的几何关系，可由一张面元的参数简单解求另一张面元的参数，没有增加计算的复杂度。

（3）采用与单面元最小二乘匹配相同的几何约束和误差方程的建立方式，在突出垂直边缘的情况下没有损失非边缘区域的匹配精度。

与传统的最小二乘匹配方法相比，本发明的方法都具有很明显的优势，既有较少的未知数个数，又有较高的匹配精度，更突出了三维模型中的垂直边缘，是一种可行的影响匹配方法。

以上内容是结合遥感影像和最佳实施方案对本发明说做的进一步详细说明，不能认定本发明的具体实施只限于这些说明。本领域的技术人员应该理解，在不脱离由所附权利要求书限定的情况下，可以在细节上进行各种修改，都应当视为属于本发明的保护范围。

Claims (4)

1.一种基于物方垂直双面元的最小二乘匹配方法，其特征在于包括以下步骤：

步骤1，在密集匹配的三维点云数据中，建立每个三维点的初始法向量，

包括以每个三维点为初始法向量的起点，以(α,β)、(π/2-α,π/2-β)为初始方向角，建立两条垂直的初始法向量，

a₁＝cosαcosβ

b₁＝sinαcosβ

c₁＝sinβ(1)

a₂＝cos(π/2-α)cos(π/2-β)＝sinαsinβ

b₂＝sin(π/2-α)cos(π/2-β)＝cosαsinβ

c₂＝sin(π/2-β)＝cosβ

其中，(a₁,b₁,c₁)、(a₂,b₂,c₂)表示两条垂直的初始法向量的坐标；

步骤2，基于各三维点的初始法向量，建立误差方程，进行双面元参数的迭代优化，搜索最佳匹配点位，实现包括以下步骤，

步骤2.1，将三维点(X_c,Y_c,Z_c)投影到参考影像上，得到参考影像上相应像点坐标(x_0c,y_0c)；

步骤2.2，以(x_0c,y_0c)为中心，一个像素为采样间隔，从参考影像上取出一个μ×μ个像素大小的影像窗口，计算窗口内每个像素的像点坐标(x₀,y₀)，其中μ为预设值；

步骤2.3，为三维点在物方建立一对互相垂直的面元P1、P2，面元中心为三维点的物方坐标值(X_c,Y_c,Z_c)，P1、P2的法向量为(a₁,b₁,c₁)、(a₂,b₂,c₂)；建立时采用计算垂直面元的方程如下，

P1：a₁(X-X_c)+b₁(Y-Y_c)+c₁(Z-Z_c)＝0

P2：a₂(X-X_c)+b₂(Y-Y_c)+c₂(Z-Z_c)＝0

其中，(X,Y,Z)表示面元P1、P2上的点；

步骤2.4，分别计算步骤2.2所取影像窗口中的μ×μ个点投影到P1、P2上每个投影点的物方坐标(X₁,Y₁,Z₁)、(X₂,Y₂,Z₂)；计算公式如下，

$\left(\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right)=R\×\left(\begin{array}{c}{x}_0\\ y_0\\ -f_0\end{array}\right)$

${\mathrm{\λ}}_{1\left(2\right)}=\frac{a_{1\left(2\right)}(X_c-X_{s0})+b_{1\left(2\right)}(Y_c-Y_{s0})+c_{1\left(2\right)}(Z_c-Z_{s0})}{a_{1\left(2\right)}{u}_0+b_{1\left(2\right)}{v}_0+c_{1\left(2\right)}{w}_0}$

$\left(\begin{array}{c}{X}_{1\left(2\right)}\\ Y_{1\left(2\right)}\\ Z_{1\left(2\right)}\end{array}\right)=\left(\begin{array}{c}{X}_{s0}\\ Y_{s0}\\ Z_{s0}\end{array}\right)+{\mathrm{\λ}}_{1\left(2\right)}\×\left(\begin{array}{c}{u}_0\\ v_0\\ w_0\end{array}\right)$

其中，

(u₀,v₀,w₀)为影像窗口中各点的像空间辅助坐标；

λ₁₍₂₎表示λ₁或λ₂，分别为影像窗口投影到P1、P2上的投影系数；a₁₍₂₎、b₁₍₂₎、c₁₍₂₎，表示a₁或a₂,b₁或b₂,c₁或c₂；

(X_s0,Y_s0,Z_s0)为参考影像外方位元素中的线元素，(X₁₍₂₎,Y₁₍₂₎,Z₁₍₂₎)表示影像窗口投影到P1、P2上的物方坐标(X₁,Y₁,Z₁)或(X₂,Y₂,Z₂)；

步骤2.5，分别计算P1、P2上每个投影点投影到各搜索影像上的像点坐标(x_i,y_i)，得到搜索影像上的相应影像窗口，其中投影到第i张搜索影像上的像点坐标记为(x_i,y_i)；

步骤2.6，分别对参考影像和搜索影像上影像窗口中的点采用双线性内插法进行重采样，得到参考影像上影像窗口中的点重采样后的像素亮度值g₀(x₀,y₀)和P1上的投影点投影到第i张搜索影像上的像点重采样后的像素亮度值g_i(x_i1,y_i1)、P2上的投影点投影到第i张搜索影像上的像点重采样后的像素亮度值g_i(x_i2,y_i2)；

步骤2.7，根据共线方程约束和最小二乘匹配原理建立误差方程如下，

$v_{P1}={\mathrm{dh}}_{0i}+g_i(x_{i1},y_{i1})\·{\mathrm{dh}}_{1i}+h_{1i}(\frac{\∂g_{i1}}{\∂x_{i1}}{\mathrm{dx}}_{i1}+\frac{\∂g_{i1}}{\∂y_{i1}}{\mathrm{dy}}_{i1})-\left(g_0\right(x_0,y_0)-h_{0i}-h_{1i}\·g_i(x_{i1},y_{i1}\left)\right)---\left(15\right)$

$v_{P2}={\mathrm{dh}}_{0i}+g_i(x_{i2},y_{i2})\·{\mathrm{dh}}_{1i}+h_{1i}(\frac{\∂g_{i2}}{\∂x_{i2}}{\mathrm{dx}}_{i2}+\frac{\∂g_{i2}}{\∂y_{i2}}{\mathrm{dy}}_{i2})-\left(g_0\right(x_0,y_0)-h_{0i}-h_{1i}\·g_i(x_{i2},y_{i2}\left)\right)$

其中，v_P1、v_P2为投影误差；h_0i、h_1i为第i张搜索影像到参考影像的辐射畸变系数；

步骤2.8，进行辐射畸变改正，解求得到未知数的改正值，包括辐射畸变改正值dh_0i、dh_1i，面元角度改正值2个dα、dβ和物点坐标改正值dZ_c；

步骤2.9，当所有改正值小于预设改正值阈值时，停止迭代，取搜索影像窗口中心点为最佳匹配点位；否则，根据改正值得到改正后的参数h_0i、h_1i、Z_c、α、β，基于改正后的参数α、β根据公式(1)计算改正后的面元法向量作为新的(a₁,b₁,c₁)、(a₂,b₂,c₂)，重复步骤2.3至2.8，并将迭代次数计算值+1，当迭代次数大于预设次数阈值仍未出现所有改正值小于预设改正值阈值时，判断迭代失败，没有找到该三维点的匹配点。

2.根据权利要求1所述基于物方垂直双面元的最小二乘匹配方法，其特征在于：步骤2.8中解求未知数的改正值方式为，用差分代替微分来计算亮度值在x、y方向的导数值，将解求的投影系数带入光束方程，联立共线方程解求搜索影像中像点坐标的改正值；简化未知数，求解改正值。

3.根据权利要求1或2所述基于物方垂直双面元的最小二乘匹配方法，其特征在于：步骤2.9中，停止迭代后根据下式计算影像间的相关系数r_i1、r_i2，

$r_{i1}=\frac{\mathrm{\Σ}\left(g_0\right(x_{0r},y_{0s})-\stackrel{\‾}{g_0\left(x_0,y_0\right)})\left(\right(g_{i1}\left(x_{i1r},y_{i1s}\right)-\stackrel{\‾}{g_{i1}\left(x_{i1},y_{i1}\right)}\left)\right)}{\sqrt{\mathrm{\Σ}{\left(g_0\right(x_{0r},y_{0s})-\stackrel{\‾}{g_0\left(x_0,y_0\right)})}^2\·\mathrm{\Σ}{\left(g_{i1}\right(x_{i1r},y_{i1s})-\stackrel{\‾}{g_{i1}\left(x_{i1},y_{i1}\right)})}^2}}$

$r_{i2}=\frac{\mathrm{\Σ}\left(g_0\right(x_{0r},y_{0s})-\stackrel{\‾}{g_0\left(x_0,y_0\right)})\left(\right(g_{i2}\left(x_{i2r},y_{i2s}\right)-\stackrel{\‾}{g_{i2}\left(x_{i2},y_{i2}\right)}\left)\right)}{\sqrt{\mathrm{\Σ}{\left(g_0\right(x_{0r},y_{0s})-\stackrel{\‾}{g_0\left(x_0,y_0\right)})}^2\·\mathrm{\Σ}{\left(g_{i2}\right(x_{i2r},y_{i2s})-\stackrel{\‾}{g_{i2}\left(x_{i2},y_{i2}\right)})}^2}}---\left(21\right)$

其中，r_i1、r_i2分别为参考影像上影像窗口与从面元P1、P2投影到第i张搜索影像所得影像窗口间的相关系数，r、s为像素在影像窗口中的行列号；g₀(x_0r,y_0s)为参考影像上影像窗口中的点重采样后的像素亮度值，其中(x_0r,y_0s)表示参考影像中影像窗口内的第r行第s列像点坐标；g_i1(x_i1r,y_i1s)、g_i2(x_i2r,y_i2s)分别为面元P1、P2上的投影点投影到第i张搜索影像上的像点重采样后的像素亮度值，其中(x_i1r,y_i1s)、(x_i2r,y_i2s)分别为从面元P1、P2投影到第i张搜索影像所得影像窗口内第r行第s列的像点坐标；分别为g₀(x_0r,y_0s)、g_i1(x_i1r,y_i1s)和g_i2(x_i2r,y_i2s)的均值。

4.根据权利要求1或2所述基于物方垂直双面元的最小二乘匹配方法，其特征在于：第k次迭代执行步骤2.9时，根据改正值得到改正后的参数h_0i、h_1i、Z_c、α、β如下式，

h_0i ^k＝h_0i ^k-1+dh_0i ^k+h_0i ^k-1dh_1i ^k(22)

h_1i ^k＝h_1i ^k-1+h_1i ^k-1dh_1i ^k

Z_c ^k＝Z_c ^k-1+dZ_c ^k

α^k＝α^k _- ¹+dα^k(23)

β^k＝β^k-1+dβ^k

其中，h_0i ^k、h_1i ^k表示第k次迭代求得的参数h_0i、h_1i的值，dh_0i ^k、dh_1i ^k表示第k次求得的参数h_0i、h_1i的改正值；Z_c ^k、α^k、β^k分别表示第k次迭代求得的参数Z_c、α、β的值，dZ_c ^k、dα^k、dβ^k分别表示第k次求得的参数Z_c、α^k、β^k的改正值。