提取光流的最佳算法（论文推荐田峰联合星载光学与SAR图像的城市大面积建筑物高度快速提取）

《测绘学报》

构建与学术的桥梁 拉近与权威的距离

田峰^1,2, 李虎^2,3

1. 新疆大学资源与环境科学学院, 新疆 乌鲁木齐 830000;

2. 新疆维吾尔自治区卫星应用中心, 新疆 乌鲁木齐 830000;

3. 滁州学院, 安徽 滁州 239000

收稿日期：2016-11-24；修回日期：2017-05-31

基金项目：高分辨率对地观测系统重大专项（95-Y40B02-9001-13/15）；新疆青年科技创新人才培养工程（2013701035）；安徽省高等学校自然科学研究项目（KJ2016A526）；滁州学院科研启动基金项目（2014QD029）

第一作者简介：田峰(1992—), 男, 硕士生, 研究方向为城市多光谱与SAR遥感应用, 国产遥感卫星应用。E-mail: 815240572@qq.com

通信作者：李虎. E-mail: lihu2881@aliyun.com

摘要：星载高分辨率光学图像与SAR图像广泛应用于城市建筑物高度提取，但光学图像存在缺少相关卫星参数的情况，而SAR图像则存在散射特征不完整以及提取效率低等缺陷。针对以上问题，本文提出一种联合高分辨率星载光学与SAR图像的城市大面积建筑物高度快速提取方法。首先，结合支持向量机（SVM）和形态学阴影指数（MSI）快速提取光学图像中的阴影并自动测量阴影长度；之后选择多个合适样本，基于模型匹配法从SAR图像中提取高度；最后将高度与阴影长度作线性回归分析，建立数学模型来提取其他建筑物的高度。该方法将不同卫星系统的数据和特征相结合，互相弥补各自缺陷，不仅提高了效率、降低了成本，同时满足精度要求。

Fusion of Spaceborne Optical and SAR Images for Building Height Quick Extraction in Big Urban Areas

TIAN Feng^1,2, LI Hu^2,3

Abstract: The spaceborne high-resolution optical images and synthetic aperture radar (SAR) images are applied to extract building height in urban areas widely. But the lack of optical satellite parameters along with the SAR images' incomplete scattering characteristics and inefficient extraction make the application flawed. To cure the above problems, we investigated the joint use of the spaceborne high-resolution optical images and SAR images to extract building height information quickly in big urban areas. The chain is decomposed into the main following steps: First, the building shadows are extracted by integrating support vector machines (SVM) with morphological shadow index (MSI) and their lengths are measured automatically. Then, the height extraction from SAR images based on a model matching technique for some appropriate samples. Finally, obtain the other heights based on the simple linear regression analysis. This approach which combines the data and feature from different satellite systems to make up the flaws for each other is not only efficient and low-cost, but also satisfy the basic accuracy requirement.

Key words: spaceborne high-resolution image building height big areas quick extraction MSI model match technique simple linear regression

随着卫星遥感技术的发展，呈现出多平台、多数据源、多分辨率以及全天时全天候观测等特点，并且成本相对较低，成为城市动态监测与信息提取的重要手段。建筑物高度作为城市遥感监测的重要内容，从侧面反映了城市的纵向信息，为研究城市空间发展过程与城市扩张驱动力因子提供了依据^[1]。

目前，用于城市建筑物高度提取的卫星数据源主要包括光学图像与合成孔径雷达(synthetic aperture radar, SAR)图像。由于成像原理不同，导致提取方式与特征选取各不相同。光学图像提取建筑物高度主要有立体像对与阴影长度两种方式，基于阴影的方式因其只需单幅数据而被广泛研究和使用。早期，文献[2]基于SPOT全色图像直接根据阴影估计建筑物高度；之后文献[3]则利用阴影成像几何模型提取建筑物高度；文献[4-5]提出了形态学阴影指数(morphological shadow index，MSI)，相较于监督分类有较好的提取效果；文献[6]则将MSI与面向对象分类相结合来提取建筑物高度。SAR图像获取建筑物高度主要分为单幅图像提取、立体像对反演^[7]、干涉测量^[8]以及后向散射强度^[9]定量估计4种方式。其中立体像对与干涉测量需要两幅以上的图像，而基于后向散射强度则需要知道建筑物表面的粗糙度以及介电常数，无疑增加了提取难度；单幅SAR图像，不仅拥有较高的分辨率与较低的成本，而且无需事先获取相关参数，并能反演出精度较高的结果。基于单幅图像提取包括直接测量^[10-11]与模型匹配^[12-14]两种方式，其中模型匹配法因其进行了多次迭代计算，可以自动获得最佳的结果，精度和抗噪性能都优于手动测量。传统的模型匹配法常涉及SAR成像模拟，增加了计算复杂度，文献[15-16]提出了通过正交投影模型实现匹配模型的快速生成。文献[17]则提出了新的思路，将光学图像中获取的参数结合假设的高度生成仿真图像，与原图匹配来获取高度。上述的提取方式各有优劣，基于阴影的光学图像提取的缺陷在于某些卫星的卫星高度角与方位角无从获取，难以利用阴影几何模型获取高度，需实地测量高度并与阴影长度建立线性关系来获得结果，该方法费时费力。而基于模型匹配的SAR图像高度提取的缺陷有以下几点：第一，城市区域结构复杂，容易造成SAR图像中的散射特征不完整；第二，某些被遮挡的建筑物无法形成强散射特征；第三，模型匹配法目前只能对仅存在一栋建筑物的图像切片进行操作，增加了大面积提取难度。

针对以上问题，本文提出一种联合星载高分辨率光学与SAR图像的大面积建筑物高度快速提取方法：首先，结合支持向量机(support vector machines, SVM)与MSI快速提取光学图像建筑物阴影并自动测量阴影长度；然后，选定样本并生成平行四边形高亮特征几何模型，基于遗传模拟退火算法迭代匹配模型与图像以获得高度；最后，将样本高度与阴影长度进行线性回归分析，依据线性相关参数提取其他建筑物高度。采用GF-2和TerraSAR-X卫星数据对本方法进行试验验证，结果显示本文方法可以有效且高效地提取大面积建筑物高度。本文方法不仅改进了算法和流程，并且采用了国产高分辨率卫星数据，为今后的国产高分辨率卫星应用进行探索实践。

1 光学图像阴影提取及长度测量1.1 阴影提取

1.1.1 阴影提取方法及流程

本文将SVM与MSI结合提取光学图像建筑物阴影。SVM是建立在统计学习理论的VC维理论和结构风险最小化准则基础上的一种学习方法, 可以较好地解决小样本、非线性、高维数据和局部极小等实际问题^[18]。然而SVM提取的结果中会出现大量错分现象，需要采用约束条件进行完善。面积约束可剔除错分的小斑块；MSI的约束则使建筑物阴影更加完整，方法是将剩下的阴影区域膨胀5个像素，再结合MSI获取最终结果。如此将监督分类和非监督分类相结合，省去大量特征阈值确定过程，使阴影提取高效而准确，提取流程如图 1所示。

图 1 阴影提取流程Fig. 1 Flow chart of shadow extraction

图选项

1.1.2 形态学阴影指数特征构建

数学形态学可以有效描述区域特征，能提取图像中感兴趣的成分。形态学开、闭运算常用来提取图像的明暗特征，并保持图像整体特征相对稳定，形态学开、闭重构则能更好地减少形状损失，适合应用于高分辨率图像^[19]。MSI建立在光谱特征和一系列形态学操作的基础之上，充分考虑了形状、尺寸、亮度、对比度和方向等特征。主要操作为形态学闭运算的黑帽重构，能够平滑图像中小于结构体的暗细节信息，提取暗结构而剔除亮结构。MSI的构建过程是：① 计算亮度，即每个像素的最大灰度值；② 对亮度图像进行闭运算，即实现黑帽重构；③ 设定线性结构体元素，有利于区分道路对阴影的混淆；④ 构建微分形态学剖面(differential morphological profiles，DMP)；⑤ 取DMP的均值，构建MSI。

1.1.3 阴影长度测量

阴影长度测量采用了一种简单可行且高效的方式^[20]：① 将提取的阴影进行标记；② 沿着太阳方位角等距生成平行线；③ 提取与阴影部分重叠的平行线并注以相同的标记；④ 自动寻找标记相同的平行线，分别计算平行线的像元个数；⑤ 并运用算法获取其众数，去除无效值，求取平均值；⑥ 利用像元个数乘以图像分辨率的方式计算出最终的阴影长度^[21]。该方式充分考虑了太阳的影响，且能快速提取大面积建筑物阴影长度。

2 星载SAR图像建筑物高度提取2.1 建筑物散射特性分析

高分辨率星载SAR图像中，建筑物主要表现为强散射的高亮特征和阴影特征，但由于城市区域结构复杂，一般阴影特征不明显。高亮特征包括亮线和亮条带，亮线由墙地二次散射引起，可以用来确定墙角的位置；而亮条带则由建筑物墙面上多个复杂结构形成的三面角反射器引起亮点组成。图 2为某建筑物的GF-2图像、TerraSAR-X原始SAR图像以及经过预处理的SAR图像切片，其中TerraSAR-X卫星沿方位向降轨右视成像，斜距图像与实际地面呈现左右镜像关系，经过预处理后可与实际地面相对应。由于主墙面上存在大量窗户结构和金属材料，形成明显的平行四边形的亮条带特征，而墙地二面角形成的二次亮线则淹没在亮条带中无法识别；侧墙面上无过多复杂结构，形成的亮点特征不明显，但二次亮线较明显^[22]。图 2(a)中的a、b、c3点分别对应建筑物墙角的3个拐点，在SAR图像中所对应的位置为图 2(b)、(c)中的a′、b′、c′ 3点和a″、b″、c″ 3点。

图 2 同一建筑物光学与SAR图像切片Fig. 2 Optical image and SAR image of the same building

图选项

2.2 模型匹配反演建筑物高度

通过建立高亮特征几何模型，并与真实SAR图像进行匹配，获得最佳匹配度的模型参数即可以反映实际建筑物尺寸。如此就将建筑物尺寸提取问题转换成模型匹配问题。模型匹配的过程分3个步骤：① 模型的建立；② 匹配函数设计；③ 迭代匹配寻优。匹配函数是判断模型与图像匹配程度的标准，是整个过程的关键，本文从区域和轮廓两个方面设计匹配函数；迭代匹配寻优算法则选用遗传算法，遗传算法对一些非线性、多模型、多目标的函数优化问题均有较好的效果，为防止遗传算法出现早熟现象，引入退火机制，兼顾了全局搜索和局部搜索且避免了早熟现象的发生^[23]。

2.2.1 高亮散射特征几何模型建立

建筑物墙面与SAR图像中的高亮特征之间存在仿射变换关系，通过建立建筑物几何参数与SAR图像特征之间的定量关系，即可反演建筑物的几何参数。图 3为建筑物与预处理后SAR图像成像几何关系，平行四边形B′F′G′C′表示亮条带，多段线A′B′C′表示二次亮线。线段B′C′的长度代表建筑物长度l，线段A′B′的长度代表建筑物宽度w，角ω反映建筑物的走向，建筑物高度h=|B′F′|×tanθ，θ为电磁波的入射角。

图 3 建筑物与SAR图像成像几何关系Fig. 3 The geometrical relationship of building and SAR image

图选项

2.2.2 匹配函数设计

匹配函数从区域和轮廓两个方面入手。区域包括灰度均值和灰度分布，轮廓则采用边界信息。

2.2.2.1 灰度均值

高亮特征区域与周围的灰度差异较大，因此可用灰度均值的差值作为匹配函数，公式^[16]如下

(1)

式中，g₁和g₂分别表示模型区域内部与外部的灰度均值。

2.2.2.2 灰度分布

Hellinger距离是用于衡量概率分布间的相似度的一种距离度量方式，因此可以用来描述灰度分布差异，若模型内部的灰度分布一致性最好，与周围区域的分布差异最大，Hellinger距离也达到最大，公式^[16]如下

(2)

式中，p_i和q_i分别表示模型内部与外部灰度值为i的像元出现的概率。

2.2.2.3 边界信息

高亮区域与周边区域间存在明显的边界，可以利用梯度信息来反映匹配的程度，公式^[16]如下

(3)

式中，L表示模型轮廓；grad(j)表示像元j的梯度值。

以上匹配函数单独使用无法全面的判断匹配程度，将3个匹配函数加权组合形成新的匹配函数可以满足要求，公式如下

(4)

式中，w₁、w₂、w₃分别表示权重因子。对匹配函数进行归一化处理，使3个匹配函数保持相同的数量级，权重因子可视数量级自适应调节。

2.2.3 遗传模拟退火算法迭代寻优

完成模型建立和匹配函数设计后，将运用遗传模拟退火算法进行迭代寻优来获得最佳结果，由于本文只关注高度信息，建筑物长度、宽度和走向信息直接从高分辨率光学图像中获取，因此只需假设高度与模型的位置信息即可。遗传模拟退火算法步骤如下：

(1) 染色体编码。采用实数编码的方式，每条染色体包括3个基因，分别表示高度h和位置(x，y)，根据公式h=|B′F′|×tanθ可得到图 2模型中假设的线段B′F′的长度，(x，y)表示点F′的位置，如此可生成确定的几何模型。

(2) 初始化种群。种群内随机初始化N_p个长度为3的染色体。

(3) 设置初温T及退温参数a，温度衰减函数为Tn 1=aT_n，其中n表示进化代数。

(4) 适应度判断。适应度用来评价个体的优劣程度，用匹配函数来判断几何模型的适应度。

(5) 遗传模拟退火操作。对种群内的个体进行选择、交叉、变异操作形成新种群，并基于Metropolis准则^[24]判断是否保留当代最优个体：Mn 1和M_n分别表示新种群和上一代的最佳适应度，当Mn 1＞M_n时，则保留新种群；当Mn 1＜M_n时，则按概率exp[Mn 1-M_n/kT]接收新种群，其中k为Boltzmann常数。

(6) 终止条件。算法是否终止将依据每代种群中最佳适应度的变化情况来判断或者是否达到最大迭代次数MaxGen，当连续r代不发生变化则认为算法收敛，可以终止算法，此时最佳适应度的染色体基因值h即为反演的建筑物高度。

3 建筑物高度联合提取方法及流程3.1 提取方法

本文提出的城市大面积建筑物高度快速提取方法，在高效快速的基础上兼顾精度要求。快速高效的特点在阴影提取、阴影测量、高度获取中都有所体现。阴影提取过程中，SVM与MSI的结合省去大量特征阈值确定的过程；阴影长度测量采取自动计算方式，相比于逐个手动测量，大大提高了大面积计算的效率；模型匹配的SAR图像高度提取，同样比多次手动测量并取均值的方式更高效；联合高度提取则无需先验的实际高度测量和大量的SAR图像切片匹配计算，也显示了快速高效的特点。对于精度方面，由于所有的信息皆从遥感图像中获取，需尽可能保证信息最大限度接近实际情况，才能获取有效的线性相关参数。为满足该要求，样本的选取则显得十分重要，样本的阴影呈现要完整，在SAR图像中其亮条带与亮线特征均要明显且完整，这是该方法的关键。获取样本阴影长度与高度后，依据成像时太阳和卫星参数的确定性^[25]，阴影长度与建筑物高度近似线性关系，对长度和高度进行线性回归分析并获取系数和常数，依此线性关系进行大面积建筑物高度提取。

3.2 提取流程

主要流程包括图像预处理、阴影提取、阴影测量、样本高度提取、线性回归分析、大面积高度提取。技术路线如图 4所示。

图 4 技术路线Fig. 4 Flowchart of technology

图选项

4 试验结果及分析4.1 数据及研究区

高分辨率星载光学图像来自国产的GF-2卫星，将多光谱图像与全色图像融合后分辨率可达0.8m，成像时间为2016-04-17 13:25:49。高分辨率SAR图像来自德国的TerraSAR-X卫星，分辨率为1m的HH极化聚束式斜距图像，入射角为43.45°，成像时间为2016-05-21 00:17:48。两类数据成像时间相差大约一个月，由于各自图像的特征及参数提取是独立完成的，并不参考其他图像的成像参数，因此该因素不会影响本文方法的使用。实际建筑物高度数据来自全站仪实地测量。研究区位于乌鲁木齐市头屯河区王家沟油库东侧建筑区，具体范围如图 5所示，标记的建筑物为本次的试验对象。

图 5 研究区范围Fig. 5 Range of study area

图选项

4.2 光学图像阴影提取及结果分析

将研究区均分为3个子区域，分别提取建筑物阴影，如图 6所示为第1个子区域的提取结果。整个提取流程大约用时3.5min，而结合多尺度分割与面向对象分类的非监督方式的结果与本文类似，但用时远远大于本文方法。从图 6中可以看出SVM初提取的结果中出现大量错分现象，主要受树木与围墙的影响，而且建筑物的阴影并不完整；约束后剔除了错分现象，并使阴影更加完整。引入检测率、虚警率、漏警率来定量的分析提取结果^[26]，如表 1所示。约束后虚警率大大降低，不过稍微降低了检测率，提高了漏警率，这是由于某些建筑物阴影被植被遮挡后形成极小的阴影斑块，在面积约束时被剔除，但对整体结果影响不大。

图 6 阴影提取结果Fig. 6 Results of shadow extraction

图选项

表 1 方法改进前后结果对比Tab. 1 Comparison of the results before and after improvement

表选项

4.3 样本高度提取

依据上文提到的样本选取原则，从SAR图像中选取合适的样本，并基于模型匹配法提取建筑物高度。由于建筑物长、宽以及走向信息可以由光学图像获取，而且某些建筑物的侧墙面亮线也不明显，可能影响匹配结果，因此本文简化了模型的建立，仅建立平行四边形的几何模型匹配亮条带特征即可。本文选取5个样本，下面为其中一个样本的高度提取试验。从光学图像中获取样本建筑物长度为56m，走向方位角为21°，建立模型并采用遗传模拟退火算法迭代匹配样本原图，种群个体数N_p=100，初温T=300，退温参数a=0.2，选择操作选择轮盘赌法，交叉操作选择单点交叉法，变异操作选择单点变异法，最大迭代次数MaxGen=200，r=50。提取的高度为20.8408m，实地测量的高度为21.4080m。图 7为样本的原图和匹配结果叠加图，从中可以看出亮线特征比模型稍长，这是由于墙面比地面光滑，电磁波第一次入射墙面发生镜面反射到地面，之后在地面发生漫反射后被传感器接收，此时没有正面朝向入射方向的墙面，其第一次反射的电磁波将落在稍微靠后的地面上再发生二次散射，因而造成如此现象。图 8为进化过程中平均适应度与最佳适应度曲线。表 2为所有样本的提取高度与实际高度以及阴影长度。将提取高度与阴影长度做线性回归分析，得到系数和常数为1.016、3.213。

图 7 匹配结果示意图Fig. 7 Results of matching

图选项

图 8 匹配函数演化曲线Fig. 8 Matching function curve evolution

图选项

表 2 样本提取结果Tab. 2 Extraction results of samplesm

表选项

4.4 大面积提取结果及分析

从研究区内随机挑选24栋建筑物，这些建筑物均匀分布在研究区内，运用本文方法提取建筑物高度，结果如表 3所示。结果显示，最大绝对误差不超过4m，绝大多数误差不超过3m，表明本方法能保证精度要求。测量速率方面，如果手动测量一栋建筑物的高度和阴影长度，且获得相对准确的结果，大约需要4min的时间，则24栋建筑物则需要不断测量96min，而且是在保证SAR图像特征完整的前提下；而本方法的操作大约花费了25min，用时远远低于手动测量，而且对图像特征不完整的情况有很好的针对效果，不仅短时高效，更能最大限度地提高区域内可测建筑物的数量。针对几个误差较大的结果进行分析，主要有以下4个原因：① 阴影区域存在高大植被，影响了阴影的呈现，进而影响了阴影长度测量结果；② 某些建筑物其底部所处平面高于阴影呈现平面，造成阴影长度过长；③ 部分区域建筑物密度较大，造成阴影的重叠或遮挡，影响阴影长度的测量；④ 结构复杂的建筑物，其阴影也同样复杂，也会影响阴影长度的测量。

表 3 大面积提取结果与误差Tab. 3 Extraction results and errors of big area

表选项

5 结论

针对星载光学图像或SAR图像单独使用提取大面积建筑物高度时存在过程复杂、低效、抗噪性能低等问题，本文提出一种联合高分辨率星载光学和SAR图像的城市大面积建筑物高度快速提取方法。将支持向量机与形态学阴影指数相结合快速提取光学图像建筑物阴影，同时基于模型匹配法提取SAR图像建筑物高度，将高度与阴影长度作线性回归分析，并据此线性关系提取大面积建筑物高度。该方法将两类图像的优势相结合，弥补了各自的缺陷，改进算法以简化了提取流程，模型的建立更加适合星载图像，半自动化的方式提高了效率，同时还兼顾了精度要求，可以为大面积城市研究和决策提供良好的数据支撑，并为国产卫星数据应用进行探索实践。但有两点需要指出：① 本研究区位于市区边缘，区域内不存在高层建筑物，因此缺少对该类建筑物的试验和分析，将在今后的工作中予以完善；② 针对各种原因引起的阴影长度测量不准确现象，应该再次利用SAR图像的对该建筑物进行高度提取，进一步对结果进行约束和完善，使该方法能更好地应用于城市信息提取。

【引文格式】田峰，李虎。联合星载光学与SAR图像的城市大面积建筑物高度快速提取[J]. 测绘学报，2017，46(7)：891-899. DOI: 10.11947/j.AGCS.2017.20160573

更多精彩内容：

,