CN104732577A - 一种基于uav低空航测系统的建筑物纹理提取方法 - Google Patents

Abstract

本发明一种基于UAV低空航测系统的建筑物纹理提取方法，具体的步骤如下：第一步：建筑物轮廓线与影像自动配准；第二步：基于多视影像直线匹配的建筑物轮廓线与影像半自动配准；第三步：建筑物立面纹理提取；第四步：建筑物立面纹理几何纠正。本发明的有益效果是：(1)针对建筑物纹理映射存在的矢量线与影像不配准问题，提出了建筑物轮廓线与影像自动配准方法。(2)考虑到在多幅影像中建立建筑物空间轮廓线与其影像直线特征映射关系的复杂性，提出了基于多视几何约束的影像直线特征半自动匹配策略。

Description

一种基于UAV低空航测系统的建筑物纹理提取方法

技术领域

本发明涉及城市三维数据纹理提取方法，特别涉及一种基于UAV低空航测系统的建筑物纹理提取方法。

背景技术

建筑物纹理提取与处理手段、三维建模精度决定了城市三维模型生产技术与应用水平。目前城市三维建模技术存在的主要问题有:

①三维建模生产表现为低效率高成本。从精细三维建模技术分析，重建城市建筑物三维几何信息和纹理信息是一个困难的过程，甚至半自动的重建方式离实际运用也相差较远，基本上处于“计算机上绣花”的人工或者半自动生产水平，真三维模型的制作量巨大，是目前城市三维地理信息系统建设投资最多的一项内容，花费巨资的真三维模型其重复利用程度不高，投入产出不成比例，直接影响了城市三维地理信息系统建设进展与应用。

②三维模型应用水平较低，没有产生应有的工程价值。主要原因在于目前三维建模主要满足于地理空间信息的三维可视化浏览，侧重于用户真实体验感。从三维地理信息技术发展、数字城市/智慧城市三维(多维)地理空间框架建设需求等方面分析，应该发展具有高测图精度、真实影像纹理的三维建模技术。

近年来，无人飞行器(Unmanned Aerial Vehicle，UAV)低空航测遥感凭借其独特优势越来越受到摄影测量行业的重视，实际应用范围逐渐扩大。利用多相机组合形成等效大面阵相机，是近些年世界先进技术成就。UAV载组合相机系统在每个曝光位置，几个相机同时曝光，可从不同角度获取目标对象影像数据，并通过扩展航向视场角增大基高比，提高高程量测精度，为精细测绘提供可能性，旁向的宽视场可以增加航带影像的地面覆盖宽度，提高飞行作业效率以及减少野外控制点的布设数量。越来越多学者关注无人机低空飞行控制和组合相机影像拼接、匹配等难题，其应用也多局限于传统航测(如DEM、DOM、DLG)生产，而无人机载组合宽角相机影像在三维建模中的应用研究才刚刚开始。因此，以无人飞行器为平台，利用具有大倾角组合数码相机进行航空摄影，研究建筑墙面纹理的重建理论，对寻求经济、高效地重建建筑物立面真实纹理将具有实际的意义。

发明内容

本发明提出一种基于UAV低空航测系统的建筑物纹理提取方法，即建筑物轮廓线与影像配准方法。

为实现上述目的，本发明采用的技术方案为：

一种基于UAV低空航测系统的建筑物纹理提取方法，具体的步骤如下：

第一步：建筑物轮廓线与影像自动配准

1)首先进行“粗“定位；

利用共线条件方程(公式1)将建筑物轮廓矢量线反投于影像，获取粗略投影定位，并在周边开设一缓冲区，在此区域中进行直线提取；

2)然后进行直线提取；

利用经典Canny算子对进行边缘提取，利用建筑物轮廓线角度约束Hough直线提取算法进行直线提取；

3)建筑物轮廓线调整；

基于Canny边缘提取和Hough直线提取获取的直线特征，在建筑物轮廓线上任一空间直线上对应的是数目不等的备选直线段，首先进行候选直线段判定，其约束条件有两个，分别为：

条件一：角度限制

以任一条投影矢量线为基准，选择 $\left[\begin{array}{cc}c\tan \frac{y_1-y_2}{x_1-x_2}-A& c\tan \frac{y_1-y_2}{x_1-x_2}+A\end{array}\right]$ 作为备选直线段角度限制条件；其中，c为常数，取1；(x₁,y₁)(x₂,y₂)为矢量线上任意两点坐标；A为常数，初次循环A值设为15像素；

条件二：距离限制

以建筑物矢量线周围一定缓冲区内为基准，计算建筑物矢量线中心点坐标到备选直线段的距离作为距离限制条件；

得到候选直线段后，最后进行最佳直线段判定；

最佳直线段特征的判定，是逐一对每条建筑物轮廓线进行备选直线段特征评价，通过“带权置信度”来完成；

4)循环优化

利用获取的最佳候选直线段进行相交获取端点，以其中任一端点为标准，对应投影矢量线上端点进行坐标平移；以此矢量投影线为标准，依次进行候选直线段判断和最佳直线段判断；设立循环次数为2-4次，循环优化最佳直线段位置，保存矢量线结果；

5)修正影像内外方位元素

根据单片空间后方交会原理，以建筑物矢量轮廓线端点作为已知控制点，对应于影像上自动提取的特征线角点作为像坐标，利用共线条件方程(公式1)，修正该影像原先的3个内方位元素值和6个外方位元素值；

第二步：基于多视影像直线匹配的建筑物轮廓线与影像半自动配准

首先在第一步直线提取的基础上，在一张影像上以人工引导方式获取精确的直线特征信息；然后从其它影像上的备选直线段特征集中根据同名直线匹配方法自动获取目标直线特征，其实现方法如下：

1)基于核线束约束的候选直线组判断

利用目标直线端点所决定的核线束来缩小备选同名直线特征的搜索空间；

2)基于三视张量约束的候选直线判定

三视张量建立了三幅影像中同名直线间的约束关系。从待判定的两张影像候选直线组中依次选出两条直线与目标影像直线进行三视张量约束，计算候选直线在目标影像上的所有直线，作为对应同名直线的预测，定义预测直线与目标直线之间的误差公式，判定最佳候选直线；设目标直线为l_p，同名预测直线为l_i(i＝0,1,2...n)，表示l_i与l_p的夹角，取表示Q(x,y)(为l_p线段的中点坐标)到l_i的距离，则定义直线误差公式为：

$F\left(l_i\right)=0.5*(1-\frac{{\mathrm{\θ}}_{l_i}}{\mathrm{MAX}\left({\mathrm{\θ}}_{l_i}\right)})+0.5*(1-\frac{{\mathrm{LEN}}_{l_i{l}_p}}{\mathrm{MAX}\left({\mathrm{LMN}}_{l_i{l}_p}\right)})$

其中：MAX(*)：表示取预测直线特征中关于*的最大值。

3)基于核线约束的最小二乘影像匹配

通过三视张量约束得到每张影像对应直线，相邻直线依次相交得到建筑物轮廓线段；为了进一步提高建筑物轮廓线精度，对于每一个轮廓线顶点采用核线约束的最小二乘影像匹配；进行多重叠影像的核线几何约束，消除不位于核线几何上线段端点；

4)基于建筑物轮廓线的空间坐标调整

完成建筑物轮廓线在影像上的配准后，以影像内外方位元素和房屋边界线段端点像坐标为已知数，代入共线方程中，进行空间前方交会，计算各个角点精确的空间坐标；对于多视影像，空间前方交会的多余观测量增大，将共线方程线性化，代入最小二乘方程矩阵中，并根据建筑物本身具有的的几何约束条件进行空间前方交会，提高物方点计算稳定性和精度；

第三步：建筑物立面纹理提取

首先根据矢量数据里获取建筑物高，若原有矢量数据没有存储建筑物高度信息或高度信息不精确情况下，通过多视配准的影像建筑物边缘地面的同名点进行空间前方交会，计算建筑物边缘地面点高度值h，即为顶面高程与地面点高程差值；根据建筑物高度值h和建筑物顶面各角点坐标(X_i,Y_i,Z)，获取建筑物底座各点空间坐标(X_i,Y_i,Z-h)；通过共线方程(公式1)计算对应影像上建筑物底座点坐标值；以顶面轮廓线任一顶点为起点，依次以相邻两顶点和对应底座点组合为单位面元，作为立面纹理单位面，循环一周后可提取建筑物所有立面纹理区域；

对于高重叠度的序列影像，同一建筑物在6-10幅影像上可见，每一个建筑物立面对应着两幅以上影像，并且随着摄像角度的不同，立面纹理的清晰程度和受遮挡程度都各不相同，因此需从定向的纹理图像自动选择最佳的纹理；

对建筑物几何模型中的每个侧面T_K进行如下操作：

(1)计算该侧面的单位法向量n/|n|；

建筑物立面T_K对应顶点V₁、V₂、V₃的空间坐标为V₁(X₁，Y₁，Z₁)，V₂(X₂，Y₂，Z₂)，V₃(X₃，Y₃，Z₃)，构建的向量V₁V₂和V₁V₃垂直于法向量n，用公式表示为：

V₁V₂＝(X₂-X₁，Y₂-Y₁，Z₂-Z₁)

V₁V₃＝(X₃-X₁，Y₃-Y₁，Z₃-Z₁)

n＝V₁V₂*V₁V₃

(2)对于每幅图像i(0≤i≤M，共有M幅图像)，计算视点方向N_i，即与相机摄影中心连接矢量；

(3)计算单位法向量a和N_i的空间夹角θ_i；

通过计算获取所有目标影像的夹角θ_i，选取MIN(θ_i)对应的影像作为成像角度最好的纹理影像；目标影像的夹角：0°＜＝θ_i＜90°；

第四步：建筑物立面纹理几何纠正

1)纠正后纹理图像范围的确定

Length＝INT(L/GSD+0.5)

Width＝INT(H/GSD+0.5)

其中：Length和Width为影像长宽值，单位为像素；

GSD为地面分辨率；

L和H为建筑物空间长度和高度；

INT(*)表示取整运算；

纹理影像起点坐标(0,0)对应建筑物顶点坐标(X,Y,Z)，则纹理影像任意坐标(0+l,0+w)(其中l为x方向像素个数，w为y方向像素个数)对应建筑物空间坐标为(X+GSD*lcosθ,Y+GSD*lsinθ,Z+GSD*w)。

2)建立纠正函数

由共线方程可知，公式1其建立的是建筑物空间坐标与原始影像像平面坐标之间的数学关系；从而根据1)所建立的纠正后纹理影像与建筑物空间坐标关系，间接确立原始影像与纠正后影像间纠正函数。

$\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]=\left[\begin{array}{ccc}{a}_1& b_1& c_1\\ a_2& b_2& c_2\\ a_3& b_3& c_3\end{array}\right]\left[\begin{array}{c}X-X_s\\ Y-Y_s\\ Z-Z_s\end{array}\right]$

3)采用双线性内插，求得原始影像对应点处的灰度值并赋给纠正后像元素。

在上述方案的基础上，第一步的步骤2)中利用经典Canny算子对进行边缘提取时，设定低阈值和高阈值的比例为0.2，高阈值占图像像素总数的比例为0.9。

本发明方法所针对的建筑物对平顶矩形屋顶、人字形屋顶和半坡形屋顶的建筑物纹理提取效果更佳。

本发明的有益效果是：

(1)针对建筑物纹理映射存在的矢量线与影像不配准问题，提出了建筑物轮廓线与影像自动配准方法。对建筑物矢量线反投于影像上进行“粗”定位，利用Canny算子进行边缘提取，通过角度和距离限制，利用带权置信度自动快速获取最佳影像边缘线。对于简单建筑物自动化程度较高，可完全取代传统手工作业。

(2)考虑到在多幅影像中建立建筑物空间轮廓线与其影像直线特征映射关系的复杂性，提出了基于多视几何约束的影像直线特征自动匹配策略。在自动匹配方法的基础上，对第一张影像人工进行建筑物轮廓线配准；然后利用多视几何约束，包括核线约束、三视张量约束、相似性约束等进行其他影像与第一张影像匹配，获取其他影像上正确的建筑轮廓线。与传统手工进行建筑物纹理裁切相比，其自动化程度和效率都得到很大程度提高。

附图说明

图1是本发明建筑物纹理提取法总流程图；

图2是基于直线自动提取的建筑物轮廓线与影像配准流程图；

图3是最佳直线段的判定；

图4是基于多视影像直线匹配的建筑物轮廓线与影像配准流程图；

图5是基于核线约束的候选直线判断；

图6是三视张量约束图；

图7是四重叠度影像核线示意图；

图8地物空间坐标图；

图9建筑物视角法模型图；

图10数字图像纠正流程图；

图11实际立面空间矢量坐标图；

图12纠正后纹理影像图。

具体实施方式

第一步建筑物轮廓线与影像自动配准

如图1和2所示，

1)首先进行“粗“定位。利用共线条件方程(公式1)将建筑物轮廓矢量线反投于影像，获取粗略投影定位，并在周边开设一缓冲区，在此区域中进行直线提取。

$\begin{array}{c}x-x_0=-f\frac{a_1(X-X_s)+b_1(Y-Y_s)+c_1(Z-Z_s)}{a_3(X-X_s)+b_3(Y-Y_s)+c_3(Z-Z_s)}\\ y-y_0=-f\frac{a_2(X-X_s)+b_2(Y-Y_s)+c_2(Z-Z_s)}{a_3(X-X_s)+b_3(Y-Y_s)+c_3(Z-Z_s)}\end{array}$         (公式1)

式中：

x，y为像平面坐标系o-xy下的像方点坐标；

x₀，y₀，f为影像的内方位元素；

X_S，Y_S，Z_S为物方空间坐标系A-XYZ下的摄站点坐标；

X，Y，Z为空间坐标系A-XYZ下的物方点坐标；

a_i，b_i，c_i(i＝1，2，3)为影像的3个外方位角元素ω、κ组成的9个方向余旋。

2)然后进行直线提取。

利用经典Canny算子对进行边缘提取时，为保证有足够的边缘点参与直线特征的生成，又尽量减少各种伪边缘点和抑制无关的细节和噪声，设定低阈值和高阈值的比例为0.2，高阈值占图像像素总数的比例为0.9。为了提高直线提取的正确性和完整性，同时获取线段的端点坐标，利用建筑物轮廓线角度约束的Hough直线提取算法。

离散Hough变换公式为：

$H({\mathrm{\ρ}}_m,{\mathrm{\θ}}_n)=\underset{i=0,j=0}{\overset{i=\mathrm{IX},j=\mathrm{IY}}{\mathrm{\Σ}}}I(x_i,y_i)\mathrm{\δ}({\mathrm{\ρ}}_m,{{\mathrm{\ρ}}^{\′}}_m)$

ρ'_m＝x_i cosθ_n+y_i sinθ_n

$-\mathrm{\π}/2\≤\mathrm{\θ}\≤\mathrm{\π},0\≤\mathrm{\ρ}\≤\sqrt{{\mathrm{IX}}^2+{\mathrm{IY}}^2}$

m＝0,1,......,HX-1,n＝0,1......,HY-1

其等价权值矩阵表示为

$H({\mathrm{\ρ}}_m,{\mathrm{\θ}}_n)=\underset{i=0,j=0}{\overset{i=\mathrm{IX},j=\mathrm{IY}}{\mathrm{\Σ}}}I(x_i,y_i)h_{m,n}(i,j)$

$h_{m,n}(i,j)=\left\{\begin{array}{cc}1& |{\mathrm{\&rho;}}_m-x_i\cos {\mathrm{\&theta;}}_n-y_i\sin {\mathrm{\&theta;}}_n|<\mathrm{\&Delta;\&rho;}/2\\ 0& |{\mathrm{\&rho;}}_m-x_i\cos {\mathrm{\&theta;}}_n-y_i\sin {\mathrm{\&theta;}}_n|\&GreaterEqual;\mathrm{\&Delta;\&rho;}/2\end{array}\right.$

从上式可知，h_m,n(i,j)与输入图像的内容无关，是由输入图像的大小IX、IY唯一决定，可定义为霍夫变换的权值矩阵。这样，经过霍夫变换后的参数空间中的一个点等于该点对应的权值矩阵与输入图像的乘积。

权值矩阵中x_i、y_i和θ_n的计算方法为：

x_i＝i,(0≤i≤IX),y_i＝j,(0≤j≤IY)

θ_n＝-π/2+n×Δθ,Δθ＝arctan(1/max(IX,IY))

(0≤n＜1.5π/arctan(1/max(IX,IY))

ρ_m＝m×Δρ,(0≤m＜M)

上式中M和Δρ的取值为：

$M=\left\{\begin{array}{cc}(\mathrm{IX}\&times;\cos {\mathrm{\&theta;}}_n)/\left(\sqrt{2}\cos ({\mathrm{\&theta;}}_n+3\mathrm{\&pi;}/4)\right)& -\mathrm{\&pi;}/2\&le;{\mathrm{\&theta;}}_n<-\mathrm{\&pi;}/4\\ (\mathrm{IX}\&times;\cos {\mathrm{\&theta;}}_n)/\left(\sqrt{2}\cos ({\mathrm{\&theta;}}_n-\mathrm{\&pi;}/4)\right)& -\mathrm{\&pi;}/4\&le;\mathrm{\&theta;}<0\\ {\mathrm{\&rho;}}_{\mathrm{Max}}\cos (\mathrm{\&alpha;}-{\mathrm{\&theta;}}_n)/\left(\sqrt{2}\cos ({\mathrm{\&theta;}}_n+\mathrm{\&pi;}/4)\right)& 0\&le;{\mathrm{\&theta;}}_n<\mathrm{\&pi;}/4\\ {\mathrm{\&rho;}}_{\mathrm{Max}}\cos (\mathrm{\&alpha;}-{\mathrm{\&theta;}}_n)/\left(\sqrt{2}\cos ({\mathrm{\&theta;}}_n-3\mathrm{\&pi;}/4)\right)& \mathrm{\&pi;}/4\&le;{\mathrm{\&theta;}}_n<\mathrm{\&pi;}/2\\ (\mathrm{IY}\&times;\cos ({\mathrm{\&theta;}}_n-\mathrm{\&pi;}/2))/\left(\sqrt{2}\cos ({\mathrm{\&theta;}}_n-\mathrm{\&pi;}/4)\right)& \mathrm{\&pi;}/2\&le;{\mathrm{\&theta;}}_n<3\mathrm{\&pi;}/4\\ (\mathrm{IY}\&times;\cos ({\mathrm{\&theta;}}_n-\mathrm{\&pi;}/2))/\left(\sqrt{2}\cos ({\mathrm{\&theta;}}_n+3\mathrm{\&pi;}/4)\right)& 3\mathrm{\&pi;}/4\&le;{\mathrm{\&theta;}}_n<\mathrm{\&pi;}\end{array}\right.$

$\mathrm{\&Delta;\&rho;}=\left\{\begin{array}{cc}\sqrt{2}\cos ({\mathrm{\&theta;}}_n+3\mathrm{\&pi;}/4)& -\mathrm{\&pi;}/2\&le;{\mathrm{\&theta;}}_n<-\mathrm{\&pi;}/4\\ \sqrt{2}\cos ({\mathrm{\&theta;}}_n-\mathrm{\&pi;}/4)& -\mathrm{\&pi;}/4\&le;{\mathrm{\&theta;}}_n<0\\ \sqrt{2}\cos ({\mathrm{\&theta;}}_n+\mathrm{\&pi;}/4)& 0\&le;{\mathrm{\&theta;}}_n<\mathrm{\&pi;}/4\\ \sqrt{2}\cos ({\mathrm{\&theta;}}_n-3\mathrm{\&pi;}/4)& \mathrm{\&pi;}/4\&le;{\mathrm{\&theta;}}_n<\mathrm{\&pi;}/2\\ \sqrt{2}\cos ({\mathrm{\&theta;}}_n-\mathrm{\&pi;}/4)& \mathrm{\&pi;}/2\&le;{\mathrm{\&theta;}}_n<3\mathrm{\&pi;}/4\\ \sqrt{2}\cos ({\mathrm{\&theta;}}_n+3\mathrm{\&pi;}/4)& 3\mathrm{\&pi;}/4\&le;{\mathrm{\&theta;}}_n<\mathrm{\&pi;}\end{array}\right.$

上式中 ${\mathrm{\&rho;}}_{\mathrm{Max}}=\sqrt{{\mathrm{IX}}^2+{\mathrm{IY}}^2},$ α＝tan^-1(IY/IX)。

利用建筑物轮廓线角度约束的Hough直线提取算法主要步骤为：

(a)计算以每条“粗”定位投影线的角度值一定缓冲区域作为阈值TAi(i＝1，2，3.)-TBi(i＝1，2，3…)，同时计算对最短投影线所在的二值图像进行Hough变换，获取其峰值作为阈值TC；

(b)对累加器数组按照从大到小的顺序排序；

(c)若累加器数组中的最大值大于阈值TC，并且其参数空间中的角度位于阈值TAi(i＝1，2，3.)-TBi(i＝1，2，3…)，则记录参数空间中的角度所对应的m和距离所对应的n，最后从累加器数组中删除该最大值。若该值小于阈值TC，则不存在直线段，寻找次最大值进行步骤(c)处理。

(d)根据m和n值从权值矩阵中获得直线段的起点坐标和终点坐标。

3)建筑物轮廓线调整。

基于Canny边缘提取和Hough直线提取获取的直线特征，在建筑物轮廓线上任一空间直线上对应是数目不等的备选直线段，首先进行候选直线段判定，其约束条件有两个，分别为：

条件一：角度限制

以某条投影矢量线为基准，选择 $\left[\begin{array}{cc}c\tan \frac{y_1-y_2}{x_1-x_2}-A& c\tan \frac{y_1-y_2}{x_1-x_2}+A\end{array}\right]$ 作为备选直线段角度限制条件。其中，c为常数，取1；(x₁,y₁)(x₂,y₂)为矢量线上任意两点坐标；A为常数，初次循环A值设为15像素；

条件二：距离限制

以建筑物矢量线周围缓冲区内为基准(设定垂直于矢量线两个方向各20个像素距离作为缓冲区)，计算建筑物矢量线中心点坐标到备选直线段的距离作为距离限制条件。

得到候选直线段后，最后进行最佳直线段判定。

最佳直线特征的判定，是逐一对每条建筑物轮廓线进行备选直线特征评价，可通过“带权置信度”来完成，如图3所示，即：建筑物屋顶影像上任一直线段l_p对应于实际空间直线L，l_i(i＝0,1,2...n)为检测到的备选直线段，为l_i的长度，表示l_i与l_p的夹角，等于建筑物轮廓线的投影中心点P(x,y)到l_i距离。同时，为避免屋檐等双边缘的影响，取表示Q(x,y)(该点为l_p线段的中点坐标)到l_i的距离，则选取下列函数计算备选直线特征的“带权置信度”：

$\begin{array}{c}F\left(l_i\right)=p_1*\frac{{\mathrm{LEN}}_{l_i}}{\mathrm{MAX}\left({\mathrm{LEN}}_{l_i}\right)}+p_2*(1-\frac{{\mathrm{DIR}}_{l_i}}{\mathrm{MAX}\left({\mathrm{DIR}}_{l_i}\right)})\\ +p_3*(1-\frac{{\mathrm{DIST}}_{l_i}}{\mathrm{MAX}\left({\mathrm{DIST}}_{l_i}\right)})+p_4*(1-\frac{{\mathrm{DIST}}_{l_i{l}_p}}{\mathrm{MAX}\left({\mathrm{DIST}}_{l_i{l}_p}\right)})\end{array}$

其中： $\underset{j=1}{\overset{4}{\mathrm{\&Sigma;}}}{p}_j=1.0,(0<p_j<1.0);$

p_j(j＝1,2,3,4)为各备选线段的置信度，其取值分别为：p₁＝0.3,p₂＝0.3,p₃＝0.2,p₄＝0.2。

4)循环优化

利用获取的最佳候选直线段进行相交获取端点，以其中任一端点为标准，对应投影矢量线上端点进行坐标平移。以此矢量投影线为标准，依次进行候选直线段判断和最佳直线段判断。设立循环次数，循环优化最佳直线段位置，保存矢量线结果。

5)修正影像内外方位元素

根据单片空间后方交会原理，以建筑物矢量轮廓线端点作为已知控制点，对应于影像上自动提取的特征线角点作为像坐标，利用共线条件方程(公式1)，修正该影像原先的3个内方位元素值和6个外方位元素值。

第二步基于多视影像直线匹配的建筑物轮廓线与影像半自动配准

由于建筑物丰富的构建造型以及各部分细节表现差异，例如屋檐的存在导致的双边缘、高大建筑物在倾斜摄影时倾斜严重等，在自动建立建筑物空间轮廓线与其影像特征直线的对应关系时会产生歧异；同时，建筑物矢量轮廓线在量测过程中存在的误差以及数据的不完整性，例如相邻边缘不垂直、无房高，或高程值从DSM数据中获取时误差较大等，也会导致影像上自动边缘提取出现错误。因此对于以上情况，第一步第一次循环后建筑物矢量线与实际建筑物轮廓偏离较大，自动方法失败。基于纹理相机影像上目标物多视特点，利用多视几何约束进行同名直线匹配的原理可解决上述问题。在获取影像建筑物特征边缘基础上，可通过空间前方交会修正错误的建筑物矢量轮廓线，同时补充建筑物高程信息。其基本流程图如图4所示：

首先在第一步直线提取的基础上，在一张影像上以人工引导方式获取精确的直线特征信息。然后从其它影像上的备选直线特征集中根据同名直线匹配方法自动获取目标直线特征，其实现方法如下：

1)基于核线束约束的候选直线组判断

利用目标直线端点所决定的核线束来缩小备选同名直线特征的搜索空间，位于影像1上的线段I(X₁,X₂)，其端点在影像2上的核线分别为根据核点为所有核线的交点，因此以核点e和核线形成的核线束为基准，I对应的同名直线特征I′或者与核线相交，或者包含于核线束中，如图5所示。

2)基于三视张量约束的候选直线判定

三视张量建立了三幅影像中同名直线间的约束关系。从待判定的两张影像候选直线组中依次选出两条直线与目标影像直线进行三视张量约束，计算候选直线在目标影像上的所有直线，作为对应同名直线的预测，如图6所示。定义预测直线与目标直线之间的误差公式，判定最佳候选直线。误差公式的定义决定着预测直线与目标直线之间最大的相似程度，其中直线间夹角和距离是两个主要影响因素。设目标直线为l_p，同名预测直线为l_i(i＝0,1,2...n)，表示l_i与l_p的夹角，取表示Q(x,y)(为l_p线段的中点坐标)到l_i的距离，则定义直线误差公式为：

$F\left(l_i\right)=0.5*(1-\frac{{\mathrm{\&theta;}}_{l_i}}{\mathrm{MAX}\left({\mathrm{\&theta;}}_{l_i}\right)})+0.5*(1-\frac{{\mathrm{LEN}}_{l_i{l}_p}}{\mathrm{MAX}\left({\mathrm{LMN}}_{l_i{l}_p}\right)})$

其中：表示取预测直线特征中关于*的最大值。

3)基于核线约束的最小二乘影像匹配

通过三视张量约束得到每张影像对应直线，相邻直线依次相交得到建筑物轮廓线段。为了进一步提高建筑物轮廓线精度，对于每一个轮廓线顶点采用核线约束的最小二乘影像匹配。进行多重叠影像的核线几何约束，目的是消除不位于核线几何上线段端点。图7所示为4度重叠形成的6条核线，空间点对应的像点至少位于2条核线上才满足核线约束。

4)基于建筑物轮廓线的空间坐标调整

完成建筑物轮廓线在影像上的配准后，以影像内外方位元素和房屋边界线段端点像坐标为已知数，代入共线方程中，进行空间前方交会，计算各个角点精确的空间坐标。对于多视影像，空间前方交会的多余观测量增大，将共线方程线性化，代入最小二乘方程矩阵中，并根据建筑物本身具有的的几何约束条件，如垂直、平行(X方向或Y方向上两条边相等)、等距等进行空间前方交会，提高物方点计算稳定性和精度。

共线方程决定了摄影中心点、像点和物点严格的关系，由共线方程可以得：

$\left.\begin{array}{c}(x-x_0)[a_3(X-X_s)+b_3(Y-Y_s)+c_3(Z-Z_s)]\\ =-f[a_1(X-X_s)+b_1(Y-Y_s)+c_1(Z-Z_s)]\\ (y-y_0)[a_3(X-X_s)+b_3(Y-Y_s)+c_3(Z-Z_s)]\\ =-f[a_2(X-X_s)+b_2(Y-Y_s)+c_2(Z-Z_s)]\end{array}\right\}$

整理可得：

$\left.\begin{array}{c}{l}_{1}X+l_{2}Y+l_{3}Z-l_x=0\\ l_{4}X+l_{5}Y+l_{6}Z-l_y=0\end{array}\right\}$

其中：

l₁＝fa₁+(x-x₀)a₃；

l₂＝fb₁+(x-x₀)b₃；

l₃＝fc₁+(x-x₀)c₃；

l₄＝fa₂+(y-y₀)a₃；

l₅＝fb₂+(y-y₀)b₃；

l₆＝fc₂+(y-y₀)c₃；

l_x＝fa₁X_s+fb₁Y_s+fc₁Z_s+(x-x₀)X_s+(x-x₀)b₃Y_s+(x-x₀)c₃Z_s；

l_y＝fa₂X_s+fb₂Y_s+fc₂Z_s+(y-y₀)X_s+(y-y₀)b₃Y_s+(y-y₀)c₃Z_s；

将上述共线方程线性化后的带有几何约束条件平差模型为：

V＝BX-L

Cx-W_x＝0

对于约束条件方程，以建筑物外轮廓垂直约束条件为例：大多是建筑物相邻边存在直角约束，如图8所示，据矢量正交条件有：

(X_i+1-X)(X_i-X_i-1)+(Y_i+1-Y_i)(Y_i-Y_i-1)＝-l_xi＝0

线性化得

(X_i-X_i+1)dX_i-1+(X_i+1+X_i-1-2X_i)dX_i+(X_i-X_i-1)dX_i+1+

(Y_i-Y_i+1)dY_i-1+(Y_i+1+Y_i-1-2Y_i)dY_i+(Y_i-Y_i-1)dY_i+1-l_xi＝0

矩阵形式为:C₁x-W₁＝0

C₁＝[X_i-X_i+1 X_i+1+X_i-1-2X_i X_i-X_i-1 Y_i-Y_i+1Y_i+1+Y_i-1-2Y_i Y_i-Y_i-1]；

$x=\left[\begin{array}{c}{X}_{i-1}\\ X_i\\ X_{i+1}\\ Y_{i-1}\\ Y_i\\ Y_{i+1}\end{array}\right];$ W₁＝l_xi

根据带约束条件的最小二乘平差原理，迭代求解。

可得：

其中： $N_{\mathrm{bb}}=B^T\mathrm{PBW}=B^T\mathrm{PL}{N}_{\mathrm{cc}}={\mathrm{CN}}_{\mathrm{bb}}^{-1}{C}^T$

第三步建筑物立面纹理提取

利用建筑物顶面轮廓线与影像配准结果，可精确确定建筑物顶面纹理在影像上对应的纹理区域。建筑物立面纹理提取是基于建筑物顶面轮廓线配准结果。首先根据矢量数据里获取建筑物高，若原有矢量数据没有存储建筑物高度信息或高度信息不精确情况下，通过多视配准的影像建筑物边缘地面的同名点进行空间前方交会，计算建筑物边缘地面点高度值h，即为顶面高程与地面点高程差值。根据建筑物高度值h和建筑物顶面各角点坐标(X_i,Y_i,Z)，获取建筑物底座各点空间坐标(X_i,Y_i,Z-h)。通过共线方程(公式1)计算对应影像上建筑物底座点坐标值。以顶面轮廓线一顶点为起点，依次以相邻两顶点和对应底座点组合为单位面元，作为立面纹理单位面，循环一周后可提取建筑物所有立面纹理区域。

对于高重叠度的序列影像，同一建筑物往往在多幅影像上(6-10张)可见，每一个建筑物立面对应着两幅或者更多幅影像，并且随着摄像角度的不同，立面纹理的清晰程度和受遮挡程度都各不相同，因此需从定向的纹理图像自动选择最佳的纹理。

对建筑物几何模型中的每个侧面T_K进行如下操作：

(1)计算该侧面的单位法向量n/|n|；

建筑物立面T_K对应顶点V₁、V₂、V₃的空间坐标为V₁(X₁，Y₁，Z₁)，V₂(X₂，Y₂，Z₂)，V₃(X₃，Y₃，Z₃)，构建的向量V₁V₂和V₁V₃垂直于法向量n，用公式表示为：

V₁V₂＝(X₂-X₁，Y₂-Y₁，Z₂-Z₁)

V₁V₃＝(X₃-X₁，Y₃-Y₁，Z₃-Z₁)

n＝V₁V₂*V₁V₃

(2)对于每幅图像i(0≤i≤M，共有M幅图像),计算视点方向N_i，即与相机摄影中心连接矢量；

(3)计算单位法向量和N_i的空间夹角θ_i；

通过计算获取所有目标影像的夹角θ_i，选取MIN(θ_i)对应的影像作为成像角度最好的纹理影像。但实际摄影中由于遮挡情况的存在，不同视角下的遮挡程度各不相同，单纯的视角法获取的立面纹理未必是最合适的，会出现被其他纹理完全遮挡的恶劣情况，因此还需要对符合要求的纹理影像，即0°＜＝θ_i＜90°，进行遮挡检测。建筑物立面纹理存在这样的假设：不同影像上存在的相同地物越多，其相似性越大，即不同纹理影像相关性越强，表明建筑物纹理区域受遮挡的影响越小。根据这一特性，以纠正后的立面纹理作为影像匹配窗口，计算任意两张纹理影像间的相关系数ρ，对相关系数ρ最大的两幅纹理影像I₁,I₂，以视角最好的作为最优纹理影像。

$\mathrm{\&rho;}={\mathrm{\&sigma;}}_{\mathrm{xy}}/\sqrt{{\mathrm{\&sigma;}}_{\mathrm{xx}}\&CenterDot;{\mathrm{\&sigma;}}_{\mathrm{yy}}}$

$\mathrm{\&rho;}=(\mathrm{n\&Sigma;xy}-\mathrm{\&Sigma;x\&Sigma;y})/\sqrt{(\mathrm{n\&Sigma;xx}-\mathrm{\&Sigma;x\&Sigma;x})(\mathrm{n\&Sigma;yy}-\mathrm{\&Sigma;y\&Sigma;y})}$

第四步建筑物立面纹理几何纠正

建筑物立面纹理存在倾斜、不同比例投影变形等问题，导致直接作为纹理映射到三维模型时产生纹理失真，因此必须进行立面纹理几何纠正。

具体的纠正流程如图10所示；

1)纠正后纹理图像范围的确定

Length＝INT(L/GSD+0.5)

Width＝INT(H/GSD+0.5)

根据图11实际立面空间矢量坐标和图12纠正后纹理影像坐标，其中：Length和Width为影像长宽值，单位为像素。GSD为地面分辨率，根据航线规划可得，作为纠正影像采样间隔。L和H为建筑物空间长度和高度。INT(*)表示取整运算。

纹理影像起点坐标(0,0)对应建筑物顶点坐标(X,Y,Z)，则纹理影像任意坐标(0+l,0+w)(其中l为x方向像素个数，w为y方向像素个数)对应建筑物空间坐标为(X+GSD*lcosθ,Y+GSD*lsinθ,Z+GSD*w)。

2)建立纠正函数

由共线方程可知，公式1其建立的是建筑物空间坐标与原始影像像平面坐标之间的数学关系。从而根据1)所建立的纠正后纹理影像与建筑物空间坐标关系，间接确立原始影像与纠正后影像间纠正函数。

$\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]=\left[\begin{array}{ccc}{a}_1& b_1& c_1\\ a_2& b_2& c_2\\ a_3& b_3& c_3\end{array}\right]\left[\begin{array}{c}X-X_s\\ Y-Y_s\\ Z-Z_s\end{array}\right]$

3)采用双线性内插，求得原始影像对应点处的灰度值并赋给纠正后像元素。

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (2)

1.一种基于UAV低空航测系统的建筑物纹理提取方法，其特征在于：具体的步骤如下：

第一步：建筑物轮廓线与影像自动配准

1)首先进行“粗“定位；

利用共线条件方程(公式1)将建筑物轮廓矢量线反投于影像，获取粗略投影定位，并在周边开设一缓冲区，在此区域中进行直线提取；

2)然后进行直线提取；

利用经典Canny算子对进行边缘提取，利用建筑物轮廓线角度约束Hough直线提取算法进行直线提取；

3)建筑物轮廓线调整；

基于Canny边缘提取和Hough直线提取获取的直线特征，在建筑物轮廓线上任一空间直线上对应的是数目不等的备选直线段，首先进行候选直线段判定，其约束条件有两个，分别为：

条件一：角度限制

以任一条投影矢量线为基准，选择

$\left[\begin{array}{cc}c\tan \frac{y_1-y_2}{x_1-x_2}-A& c\tan \frac{y_1-y_2}{x_1-x_2}+A\end{array}\right]$ 作为备选直线段角度限制条件；其中，c为常数，取1；(x₁,y₁)(x₂,y₂)为矢量线上任意两点坐标；A为常数，初次循环A值设为15像素；

条件二：距离限制

以建筑物矢量线周围缓冲区内为基准，计算建筑物矢量线中心点坐标到备选直线段的距离作为距离限制条件；

得到候选直线段后，最后进行最佳直线段判定；

最佳直线段特征的判定，是逐一对每条建筑物轮廓线进行备选直线段特征评价，通过“带权置信度”来完成；

4)循环优化

利用获取的最佳候选直线段进行相交获取端点，以其中任一端点为标准，对应投影矢量线上端点进行坐标平移；以此矢量投影线为标准，依次进行候选直线段判断和最佳直线段判断；设立循环次数为2-4次，循环优化最佳直线段位置，保存矢量线结果；

5)修正影像内外方位元素

根据单片空间后方交会原理，以建筑物矢量轮廓线端点作为已知控制点，对应于影像上自动提取的特征线角点作为像坐标，利用共线条件方程(公式1)，修正该影像原先的3个内方位元素值和6个外方位元素值；

第二步：基于多视影像直线匹配的建筑物轮廓线与影像半自动配准

首先在第一步直线提取的基础上，在一张影像上以人工引导方式获取精确的直线特征信息；然后从其它影像上的备选直线段特征集中根据同名直线匹配方法自动获取目标直线特征，其实现方法如下：

1)基于核线束约束的候选直线组判断

利用目标直线端点所决定的核线束来缩小备选同名直线特征的搜索空间；

2)基于三视张量约束的候选直线判定

三视张量建立了三幅影像中同名直线间的约束关系。从待判定的两张影像候选直线组中依次选出两条直线与目标影像直线进行三视张量约束，计算候选直线在目标影像上的所有直线，作为对应同名直线的预测，定义预测直线与目标直线之间的误差公式，判定最佳候选直线；设目标直线为l_p，同名预测直线为l_i(i＝0,1,2...n)_， 表示l_i与l_p的夹角，取表示Q(x,y)(为l_p线段的中点坐标)到l_i的距离，则定义直线误差公式为：

$F\left(l_i\right)=0.5*(1-\frac{{\mathrm{\&theta;}}_{l_i}}{\mathrm{MAX}\left({\mathrm{\&theta;}}_{l_i}\right)})+0.5*(1-\frac{{\mathrm{LEN}}_{l_i{l}_p}}{\mathrm{MAX}\left({\mathrm{LEN}}_{l_i{l}_p}\right)})$

其中：MAX(*)：表示取预测直线特征中关于*的最大值。

3)基于核线约束的最小二乘影像匹配

通过三视张量约束得到每张影像对应直线，相邻直线依次相交得到建筑物轮廓线段；为了进一步提高建筑物轮廓线精度，对于每一个轮廓线顶点采用核线约束的最小二乘影像匹配；进行多重叠影像的核线几何约束，消除不位于核线几何上线段端点；

4)基于建筑物轮廓线的空间坐标调整

完成建筑物轮廓线在影像上的配准后，以影像内外方位元素和房屋边界线段端点像坐标为已知数，代入共线方程中，进行空间前方交会，计算各个角点精确的空间坐标；对于多视影像，空间前方交会的多余观测量增大，将共线方程线性化，代入最小二乘方程矩阵中，并根据建筑物本身具有的的几何约束条件进行空间前方交会，提高物方点计算稳定性和精度；

第三步：建筑物立面纹理提取

首先根据矢量数据里获取建筑物高，若原有矢量数据没有存储建筑物高度信息或高度信息不精确情况下，通过多视配准的影像建筑物边缘地面的同名点进行空间前方交会，计算建筑物边缘地面点高度值h，即为顶面高程与地面点高程差值；根据建筑物高度值h和建筑物顶面各角点坐标(X_i,Y_i,Z)，获取建筑物底座各点空间坐标(X_i,Y_i,Z-h)；通过共线方程(公式1)计算对应影像上建筑物底座点坐标值；以顶面轮廓线任一顶点为起点，依次以相邻两顶点和对应底座点组合为单位面元，作为立面纹理单位面，循环一周后可提取建筑物所有立面纹理区域；

对于高重叠度的序列影像，同一建筑物在6-10幅影像上可见，每一个建筑物立面对应着两幅以上影像，并且随着摄像角度的不同，立面纹理的清晰程度和受遮挡程度都各不相同，因此需从定向的纹理图像自动选择最佳的纹理；

对建筑物几何模型中的每个侧面T_K进行如下操作：

(1)计算该侧面的单位法向量n/|n|；

建筑物立面T_K对应顶点V₁、V₂、V₃的空间坐标为V₁(X₁，Y₁，Z₁)，V₂(X₂，Y₂，Z₂)，V₃(X₃，Y₃，Z₃)，构建的向量V₁V₂和V₁V₃垂直于法向量n，用公式表示为：

V₁V₂＝(X₂-X₁，Y₂-Y₁，Z₂-Z₁)

V₁V₃＝(X₃-X₁，Y₃-Y₁，Z₃-Z₁)

n＝V₁V₂*V₁V₃

(2)对于每幅图像i(0≤i≤M，共有M幅图像)，计算视点方向N_i，即与相机摄影中心连接矢量；

(3)计算单位法向量a和N_i的空间夹角θ_i；

通过计算获取所有目标影像的夹角θ_i,选取MIN(θ_i)对应的影像作为成像角度最好的纹理影像；目标影像的夹角：0°＜＝θ_i＜90°；

第四步：建筑物立面纹理几何纠正

1)纠正后纹理图像范围的确定

Length＝INT(L/GSD+0.5)

Width＝INT(H/GSD+0.5)

其中：Length和Width为影像长宽值，单位为像素；

GSD为地面分辨率；

L和H为建筑物空间长度和高度；

INT(*)表示取整运算；

纹理影像起点坐标(0,0)对应建筑物顶点坐标(X,Y,Z)，则纹理影像任意坐标(0+l,0+w)(其中l为x方向像素个数，w为y方向像素个数)对应建筑物空间坐标为(X+GSD*lcosθ,Y+GSD*lsinθ,Z+GSD*w)。

2)建立纠正函数

由共线方程可知，公式1其建立的是建筑物空间坐标与原始影像像平面坐标之间的数学关系；从而根据1)所建立的纠正后纹理影像与建筑物空间坐标关系，间接确立原始影像与纠正后影像间纠正函数。

$\left[\begin{array}{c}x-x_0\\ y-y_0\\ -f\end{array}\right]=\left[\begin{array}{ccc}{a}_1& b_1& c_1\\ a_2& b_2& c_2\\ a_3& b_3& c_3\end{array}\right]\left[\begin{array}{c}X-X_s\\ Y-Y_s\\ Z-Z_s\end{array}\right]$

3)采用双线性内插，求得原始影像对应点处的灰度值并赋给纠正后像元素。

2.根据权利要求1所述的一种基于UAV低空航测系统的建筑物纹理提取方法，其特征在于：第一步的步骤2)中利用经典Canny算子对进行边缘提取时，设定低阈值和高阈值的比例为0.2，高阈值占图像像素总数的比例为0.9。