方差分析引导的高阶马尔可夫网络及其在点云建筑物提取中的应用

郝娇娇， 倪欢， 管海燕

南京信息工程大学 遥感与测绘工程学院, 南京 210044

1 引 言

应三维城市建模、GIS 数据库更新等应用需求，建筑物提取已成为近些年的研究热点（Khoshelham 等，2010；Tomljenovic 等，2015）。ALS（Airborne Laser Scanning）技术是获取空间数据最直接的手段之一（Du等，2017），快速高效地从海量ALS 点云数据所包含的复杂空间信息中提取建筑物，成为了一系列地学应用的首要任务。

目前，ALS点云建筑物提取可分为融合提取和直接提取两类（邓飞 等，2018）。第一类方法融合了图像光谱和点云三维几何信息；优点是通过融合多源数据，实现信息互补，进而提高建筑物提取精度（Gerke 和Xiao，2014；Zarea 和Mohammadzadeh，2016）。例如，点云与图像纹理的特征匹配，不仅能有效区分相似高程地物的种类（Lai 等，2019），还可提高区域分析和监测的效率（Ullo等，2020）。在融合两类数据的基础上，常使用由粗到精的提取策略。首先结合多重特征获取初始建筑物，然后根据一定范围内的正射影像颜色信息优化建筑物边缘（杜守基 等，2018），或利用形态学运算优化建筑物候选区域（Chen 等，2020），以提高建筑物提取精度。此外，深度学习技术由于其在视觉数据上的诸多优势，被引入到该类方法。郭峰等（2020）提出一种基于Encoder-decoder 语义模型的建筑物提取方法（SegNet），将两类数据（点云与图像）嵌入到深度学习模型。有别于SegNet 的监督学习过程，Nguyen 等（2020）提出了一种融合点云与图像的大范围建筑物无监督提取方法，具有较高的自动化水平。然而，大范围图像和点云的自动配准过程会增加建筑物提取的误差传播风险（Parmehr 等，2014）；同时，两种数据分辨率和获取时间的不同，也会产生几何和光谱特征不协调问题。

第二类方法仅使用点云数据，包括监督和非监督方法。传统的监督分类方法一般会提取建筑物和其他地物的高程、纹理、法线向量、回波次数和强度特征，根据其中某一特征、特征组合或者衍生特征（邓飞 等，2018），利用SVM（Support Vector Machine；Chen 等，2019）、CRF （Conditional Random Field；Niemeyer 等，2012）、RF（Random Forest；Niemeyer 等，2013）、GMM（Gaussian Mixed Model；赵 泉 华 等，2015）、Adaboost （Wei 等，2012）等单个或多种经典机器学习算法进行建筑物提取。这些方法充分利用了建筑物的先验信息，对大量标注样本的依赖程度较低，但分类时未学习到更高层次语义，分类精度难以进一步提高（Niemeyer 等，2014）。随着深度学习技术的不断推广，基于卷积神经网络（Maltezos 等，2019）、深度残差网络（赵传 等，2020）、图几何矩卷积（Li 等，2020）的建筑物提取方法在逐步更新和完善。这类方法提高了建筑物识别精度，但前提是需要大量标注数据进行模型训练。在数据量有限的应用场景下，深度学习算法无法对数据的规律进行无偏差估计。此外，训练过程耗时较长且依赖于高性能计算设备。在这种情况下，学者们应该重视非监督的、不依赖于人工标注的识别方法。非监督方法一般根据地物的不同属性设计出特定规则，并引入一系列模式识别算法，实现建筑物与其他类别的自适应分割。其中，逆迭代数学形态学（RIMM）算法（Cheng 等，2013）、差分形态剖面（DMPs）映射（Mongus 等，2014）、投票模型（赵传 等，2017）、非流形理论（Gilani等，2018）、自顶向下（Top-down）策略（Huang 等，2018）、MRF（Markov Random Field；Wang 等，2020）、最小割算法（Min-cut；Liu 等，2020）、层次角感知（Widyaningrum 等，2020）等被引入到点云建筑物非监督提取过程中，并均在相应的实验数据上取得了较好结果。

非监督方法目前虽精度不及监督学习方法，但其不需要借助大量人工标记数据进行训练，对硬件设备要求不高，使用更为灵活。本研究聚焦非监督识别方法，在统计学框架下，实现点云建筑物提取。相比于深度学习方法，本文方法的优势和特色主要体现在3 个方面，即（1）无人工标注依赖，不存在过拟合问题，对符合预设规则的数据，泛化能力更强；（2）算法运行时间短，不依赖于高性能计算设备，能够更为快速地获取建筑物识别结果；（3）从数理统计角度发掘点云空间上下文关系，理论可解释性更强，所发现的相关关系具备可回溯性。方差分析ANOVA（Analysis of Variance）是对两组及两组以上样本进行显著性检验的统计方法（Pulido-González 等，2020），能够准确检验样本的组间差异。因此，ANOVA 方法为分析多维度特征之间的相关性提供了一条便捷途径。

基于以上分析，本文提出了一种方差分析引导的高阶马尔可夫网络，构建了一种基于方差分析P值检验的高阶无向图势函数计算方法，实现点云超体素上下文信息建模，并精准提取建筑物。

2 方差分析引导的高阶超体素马尔可夫网络

点云数据量通常较大，难以用每一个点作为节点进行相关性建模。本文将点云数据进行滤波处理，以剔除地面点；构建超体素结构，建立对象级网络与分析框架。其中地面点剔除和超体素构建采用CSF（Zhang 等，2016）和SVLA（Ni 和Niu，2020）算法。方差分析引导的高阶超体素马尔可夫网络建筑物提取过程如图1所示。该方法首先以超体素为基元计算特征，并将超体素作为网络节点。接着，使用贝叶斯高斯混合模型实现初始聚类，以获取马尔可夫网络各节点的初始状态。然后，检索超体素邻域，采用方差分析原理，结合预定义的特征，构建高阶因子，并计算节点、边和高阶势函数，形成超体素高阶马尔可夫网络。最后，采用信念传播算法，实现高阶马尔可夫网络的近似推理，以提取建筑物。

图1 方差分析引导的高阶超体素马尔可夫网络建筑物提取流程图Fig. 1 The flowchart of building extraction using ANOVA guided high-order supervoxel Markov network

2.1 基于贝叶斯高斯混合模型的节点初始状态捕捉

高斯混合模型（Chen 等，2019）采用多个高斯分布线性叠加的方式来表示样本分布情况，其收敛性受到初始假设影响，并且其常用的期望最大化EM（Expectation Maximum）解法经常出现奇异值。贝叶斯高斯混合模型可以解决以上问题，其采用变分推理思想，使用形式简单的分布来逼近高斯混合模型，并得到一种局部最优且具有确定解的近似后验分布。

假设超体素集合为S={S1，S2，…，SM}，其相应的特征集合为X={x1，x2，…，xM}，并且马尔可夫网络的节点状态空间是二值的。贝叶斯高斯混合模型的概率密度函数定义如下：

式中，

是一个多元高斯分布，C是类别数，αc为混合系数，μc为第c个高斯分布的均值向量，Σc代表第c个高斯分布的协方差矩阵。该模型的解对应参数值αc，μc，Σc。其求解通过优化证据下界ELBO（Evidence Lower BOund）过程完成。更新后的参数带入式（2），求出节点取得各个状态的概率，即每个超体素从属于各状态的概率分布表，其对应马尔可夫网络中节点的初始状态。在本文的研究中，该初始状态给出了每一个节点从属于建筑物和非建筑物的初始概率。

2.2 超体素邻域检索

超体素邻域Ni指第i个超体素Si周围的最邻近区域，其由多个与Si相邻的超体素，k=1，…，K构成。本文采用由点邻域到超体素邻域拓展的方式，确定超体素Si的邻域信息。首先采用KD 树（Buitinck 等，2013）检索Si中每一个点的邻接点。若邻接点的集合中含有其他超体素中的点，则认为该邻接点从属的超体素位于Si邻域，此时，记录该邻接点对之间的距离，并作为与Si之间的距离dk。以此类推，确定与Si相邻的所有超体素，并将其集合作为初始邻域。如图2所示，从检索点的邻接关系出发，扩展到超体素的邻接关系，最终得到超体素A的邻域包含超体素B-H。此外，本文采用与点邻域检索类似的邻域大小限制，即根据每一个邻接体素到Si之间的距离dk，取距离Si最近的Kmax个超体素构成Si的最终邻域Ni。引入该阈值Kmax，可以增强本文方法在不同点密度和超体素尺寸条件下的调节能力。例如，若超体素尺寸较大，将设置较小的Kmax，反之，则设置较大。

图2 超体素邻域示意图Fig. 2 The schematic diagram of the neighborhood of a supervoxel

Kmax个邻接超体素和当前超体素Si共同构成了高阶马尔可夫网络中每一个高阶因子，该结构填补了点云上下文信息的缺失。

2.3 基于方差分析原理的高阶因子构建

假设预先计算了T个维度的超体素特征X={x1，x2，…，xM}，其中xi={xi1，xi2，…，xiT}T。依方差分析原理，对X进行线性最小二乘方差分析检验，其观测值为超体素各维度特征向量。首先计算各维度（组）观测值均值xˉi以及全部观测值的总体均值：

式中，Kmax是当前邻域内的超体素个数，Kmax+ 1表示同时考虑了当前超体素及其邻域，则Kmax相当于该过程的自由度。n表示全部观测值个数。接着，根据均值计算各误差平方和，包括全部观测值与总均值的误差平方和SST（Sum of Squares for Total），组间误差平方和SSA（Sum of Squares for Factor A）以及组内误差平方和SSE（Sum of Squares for Error）：

为了得到无偏估计结果，需要用各平方和除以自由度（即Kmax），该结果在方差分析中被称作方差。相应的，组间方差MSA（Mean Squares for Factor A）和组内方差MSE（Mean Squares for Error）定义为

进而，检验统计量Fv表示为

最后，根据检验统计量Fv，利用多元F分布检验，找出Fv的伴随P值。该统计量Fv和P值可以有效反映超体素之间的统计相关性，将用以实现高阶因子构建，并作为势函数计算的依据。

P值的大小反应了该邻域内所有超体素特征服从同一概率分布的可能性大小，若其较大，本文将该邻域的超体素结构均纳入到同一个因子中。因该因子所包含的节点（即超体素）多于2个，所以称作高阶因子。被划分到同一高阶因子的节点之间，使用无向边相连。

本文在经典的成对马尔可夫网络基础上，增加高阶因子，构成高阶马尔可夫网络，以提高网络上下文信息表达能力。此时，高阶马尔可夫网络的联合概率分布形式如下：

式中，X={x1，x2，…，xM}是网络中所有节点构成的集合，Z*是配分函数，ψF(·)、φF(·)、φF(·)分别是节点、边、高阶势函数，E是边集合，FSH是高阶因子集合，Hl是第l个高阶因子。接着，使用信念传播算法对式（11）的联合概率分布进行最大后验MAP（Maximum A Posterior）状态推理，得到最佳状态分布：

2.4 P值引导的势函数计算

2.4.1 节点势函数

经过贝叶斯高斯混合模型的独立预测过程，各节点（超体素）从属于各类别的概率p(xm∈Lc)被计算出来，其中Lc是第c个类别对应的节点集合。则节点势函数定义如下：

式中，ψF(xm∈Lc)对应于未归一化的概率分布表，其事件空间包含C个事件（即类别）。马尔可夫网络的归一化概率，最后由配分函数Z*归一化求得，本文采用信念传播算法求取Z*。

2.4.2 边势函数

边势函数反映了相邻节点所对应超体素之间的特征相关性，相关性越高，从属于同一类别的概率越大。本文采用前述方差分析所得到的P值来确定该先验相关性。具体地，P值越大，说明该邻域内超体素之间的同质性越高。因此，若P值大于某一阈值，则该邻域内的所有超体素从属于同一类别。若P值较小，则需要对邻域内的每一对超体素进行多重对比检验（Test for multiple comparisons）。本文采用Tukey HSD（Tukey’s Honestly Significant Difference）检验方法，计算每一对超体素之间的HSD 检验Ptukey值，以确定每一对超体素之间的相关性大小。综上所述，每一条边上的先验相关信息定义如下：

式中，TRP是P值的显著性阈值，研究一般认为，TRP值小于0.05为有统计学差异。基于此，本文取TRP= 0.05。Ptukey是Tukey HSD 检验方法计算得到的P值。

P值虽然刻画了不同超体素服从同一分布的联合概率，但其仍然无法在局部区域上衡量边上相关性的大小。本文利用一种局部相关性先验概率加权该联合分布，以增强其表达能力。假设当前待确定的边为E(i，j)，则其邻域为其连接的两个节点邻域集合的并集，即S(i&j) =S(i) ∪S(j)。因各节点存在多个状态（对应于分类问题的类别数），本文为每一类别计算了一个先验概率以衡量各对邻接超体素的局部相关性。

式中，I(xm∈Lc)是指示函数：

进而，边势函数定义如下：

式（18）表达xi和xj具有相同状态时的势函数，式（19）表达xi和xj具有不同状态时的势函数。ratee是调节系数。

2.4.3 高阶势函数

高阶势函数中包含的变量对应多个节点（即超体素），参考节点和边势函数的定义，其具体形式定义如下：

式中，L(x1) =L(x2) = … =L(xKmax+1)对应于各节点取得相同状态的情况。rateh是调节系数。其中P是方差分析计算得到的P值。

2.5 基于信念传播算法的高阶马尔可夫网络推理

本文使用信念传播算法（Hazan 和Shashua，2010）对上述高阶马尔可夫网络进行最大后验状态推理。信念传播算法首先将高阶马尔可夫网络转化成因子图GF，其由因子节点和变量节点构成，因子节点和变量节点之间由无向边相连。信念传播算法通过在因子和变量节点之间迭代传递消息的方式，不断更新因子和势函数直到收敛。

信念传播算法是一种求解图模型推理问题的高效算法，与基于能量最小化原理的图割法相比，灵活度更高。信念传播算法可以对任意势函数进行优化，包括一些图割算法所不能处理的非规则（不满足次模性的）势函数（Koller 和Friedman，2009；Le^-Huu，2019；徐胜军 等，2013）。由于信念传播算法在机器学习任务中被广泛使用，本文不对其原理进行冗余叙述。

2.6 面向建筑物提取的特征计算

本文旨在利用高阶马尔可夫网络提取建筑物，因此面向建筑物提取相应特征。本文所提取的特征包括：（1）离散度，（2）法线向量倾角，（3）法线方向方差，（4）曲率和（5）维度特征。

（1）离散度。离散度衡量超体素所包含点的空间分布离散程度，其计算超体素内各点pi=(xi，yi，zi)到中心点(xm，ym，zm)距离的平均值作为离散度。即：

式中，N是超体素中点的总数。

（2）法线向量倾角。本文以超体素为单元，利用各超体素所包含的点拟合空间平面，进而计算超体素法线向量。本文以角度的余弦值来表达法线向量倾角。设第m个超体素内的点集是Sm={pi|i= 1，2，…，N}，pi=(xi，yi，zi)，那么超体素的协方差矩阵为

（3）法线方向方差。本文引用文献 （Du 等，2017）的定义，检索超体素内部每一个点的邻域点集，利用式（22）和（23）计算协方差矩阵和法线向量倾角。再计算该超体素内部所有点的法线向量倾角方差σm，用以描述超体素的平整度。

（4）曲率。对平整度的衡量，通常还可以使用曲率，其基于式（22）得到的协方差矩阵来计算。即计算Cm的3 个特征值λ0≤λ1≤λ2，则曲率为

（5）维度特征。维度特征（Demantké等，2011）使用广泛，且效果较好，本文引用其中两种，具体如下：

式中，λ0、λ1、λ2与式（24）中的定义相同。

3 实验分析与对比

3.1 实验数据

本文使用国际摄影测量与遥感学会（ISPRS）提供的参考数据集，包括Vaihingen 和Toronto 两组机载激光雷达点云（Liu 等，2020）。为了便于与其他方法对比，本文对两个数据集中的各标注区域采用类似的命名方式，即区域1—3（来自Vaihingen 数据集）和区域4—5（来自Toronto 数据集）。区域1—3的平均点密度约为4—7个点/m2，位于住宅区。区域4—5的平均点密度约为6个点/m2，位于商业区。图3展示了各区域按高程着色的点云场景。

图3 高程着色的测试数据Fig. 3 Testing data colored by elevation

区域1 包括37 栋建筑，主要由密集的古老建筑组成，形状相对复杂。区域2 包括14 栋建筑，主要由屋顶水平的高层住宅组成。区域3 包括56栋建筑，主要由独立建筑组成，屋顶结构简单，沿路有植被。区域4 包括58 栋建筑，由不同高度建筑组成，建筑物结构复杂，顶部有附着物。区域5 包括38 栋建筑，建筑物结构多样，形状复杂，顶部带有附着物；其中包括一个显著高于其他地物的高层建筑。

3.2 精度评价指标

本文采用ISPRS 指标体系（Rottensteiner 等，2014）衡量结果质量，包括Completeness、Correctness、Quality和F1指数4个评价指标，计算方法如下：

以上4 种评价指标，在采用不同量化因子时，含义存在细微差别。其中，基于面积和基于对象的指标体系最为常用。基于面积的评价指标记为Compar、Corrar、Qar、F1ar，以投影到二维平面的面积作为量化因子，该面积通常使用像素个数衡量，精度与正确分类的区域相关。此时，TP表示正确检测的建筑物面积，FP表示错误检测的建筑物面积，FN表示遗漏的建筑物面积。基于对象的评价指标记为Compobj、Corrobj、Qobj、F1obj，以建筑物个数作为量化因子，若提取对象与参考对象重叠度大于50%，则认为识别到该对象。此时，TP表示正确检测的建筑物个数，FP表示错误检测的建筑物个数，FN表示遗漏的建筑物个数。此外，Rottensteiner 等（2014）还提供了用于大型建筑物检测的评价指标，但是面积大于50 m 的建筑物提取精度普遍较高，可比性较弱，因此，不单独对大面积建筑物进行精度评价。为更细致地评估建筑物提取性能，本文增加了基于点的评价指标，即逐点根据式（27）—式（30）计算精度，记为Comppoi、Corrpoi、Qpoi、F1poi。此时，TP表示正确检测的建筑物点数，FP表示错误检测的建筑物点数，FN表示遗漏的建筑物点数。

3.3 参数设置

本文使用Ni 和Niu（2020）提出的SVLA 算法构建超体素，该算法可以获得紧致的超体素簇，并通过调整种子分辨率（Rseed）来控制超体素大小。体素的大小是影响建筑物提取结果的重要因素，若构建的点云超体素较小，则意味着点云数据中所含超体素过多，这会增加算法的计算复杂度和运行时间，降低建筑物提取效率。反之，若构建的超体素较大，虽充分顾及了上下文信息，却忽略了局部细节特征，进而影响建筑物的提取精度。表1 记录了区域1 中Rseed设置与建筑物提取精度的相互影响关系。为了简化表达，本文以基于点的评价指标为例进行分析，其他基于面积和对象的精度变化情况与基于点的精度变化情况类似。从表中可以看到，当0.9≤Rseed≤1.3 时，建筑物提取精度逐渐上升，Rseed取1.3 时，各项精度评价指标达到最大值。当Rseed＞1.3 时，随着Rseed的增大，建筑物提取的精度逐渐下降。因此，构建合适大小的超体素至关重要。本文经过实验对比，最终将区域1—3的种子分辨率设置为1.3，区域4—5的种子分辨率设置为2.0。此外，在输入数据中不可避免地会存在一部分较小的超体素，包括噪声点和小型地物，一定程度上影响了结果的精度。因此，在超体素分割完成后，需要剔除较小的超体素，即超体素内的点数小于t时，删除该超体素。如表2 所示，区域4 和区域5 高层建筑物居多，建筑物体积较大，因此，设置t的值也最大。

表1 Rseed对建筑物提取精度的影响（区域1）Table 1 The relationship between Rseed and accuracy （Area 1）

邻域检索参数Kmax（2.2 节）控制邻域内超体素的个数，其对应了高阶马尔可夫网络中高阶因子所包含的节点个数。由于对区域1—3 的超体素分割较为细碎，邻域内包含的超体素较多，因此，设置Kmax为9。相应的，区域4—5 的体素块较大，相邻区域包含的超体素较少，设置Kmax为7。在2.4 节计算边和团势函数的过程中，经过粗对比，调节系数ratee和rateh在0.1—1.5时，建筑物提取精度较高。图4表示不同ratee和rateh取值对应的平均F1poi值。由图4 可知，在区域1—3 上，ratee和rateh的值为0.3 和1.0 时，平均F1poi值最高；在区域4—5 上，ratee和rateh的值为0.5 和1.5 时，平均F1poi值最高。各区域参数设置见表2。

表2 实验参数设置Table 2 Parameter settings of this experiment

图4 ratee和rateh的参数分析Fig. 4 Parameter analysis of ratee and rateh

3.4 实验结果

图5展示了本文方法在不同区域的建筑物提取结果。从图中可以看到，场景中的建筑物提取完整，屋顶内部无缺失，建筑物边界清晰。大型建筑物（区域4、区域5）内部无明显噪声。植被覆盖密集区域（区域2）建筑物边界明显，建筑物完整。小型建筑物（区域3）提取完整，周围无明显噪声，建筑物与篱笆、道路的边界明显。图6为本文建筑物提取的投影面积可视化结果情况，该结果对应了精度评价指标Compar、Corrar、Qar和F1ar；其中黄色、红色、蓝色分别代表TP、FN和FP，可以看到，本文方法正确提取了大部分建筑物。区域1 和2 中FN的面积最小，说明此区域建筑物提取最为完整。然而，提取结果中还是存在少量误差，多余提取和提取不完整的误差分别为FP误差和FN误差。其中FP误差主要来源是场景边缘的树木点（区域1 蓝框）、离散点（区域3、4 蓝框）以及低矮灌木点（区域2蓝框），FN误差主要存在于结构复杂的建筑物区域，其中包括部分建筑物顶部（图6 的绿框）和部分建筑物立面（图6黑框）。这是因为建筑物顶部附件较复杂的区域，不符合预设建筑物表面由平面构成的前提；此外，部分建筑物立面点分布稀疏，与预设中建筑物点密度较高有所出入。另一方面，方差分析的引入使得不同超体素特征差异更为明显，显著提高了非监督识别的整体精度。但与此同时，方差分析的引入使本文方法的识别过程对超体素特征差异的依赖度更高，进一步增加了不符合预设物体被误判的概率。因此，如何更好地进行特征选择和特征权重设置，是后续研究的重点。表3列出了这两个数据集上建筑物提取结果的各项精度值。可以看到，所有精度值都保持在80%以上，区域1—3的精度优于区域4—5，说明本文方法对住宅区的提取效果更佳。综上所述，本文方法能有效提取不同场景下的建筑物且精度较高，对于低层建筑主导的住宅区提取效果最佳。

表3 本文方法建筑物提取结果精度情况Table 3 The accuracy values of building extraction results produced by our method

图5 基于高阶超体素马尔可夫网络的建筑物提取结果Fig. 5 Building extraction results based on the high-order supervoxel Markov network

图6 建筑物提取投影面积可视化结果Fig. 6 Visualization of building extraction results at the per-area level

3.5 消融分析

本文提出方差分析引导的高阶马尔可夫网络，并将其应用于点云建筑物提取，该网络可以分解为成对马尔可夫网络、高阶因子、贝叶斯高斯混合先验和方差分析4个模块，其中成对马尔可夫网络结构是本文方法框架的基础。本节选取区域2数据进行消融实验，通过对比建筑物提取的可视化结果和精度，验证各模块的有效性。

图7展示了删除特定模块后建筑物提取结果和投影面积可视化结果，其中为了便于展示结果细节，图中选取区域2的一个局部区域进行分析。方法1采用成对马尔可夫网络实现建筑物提取，删除了本文方法中的方差分析、高阶因子和贝叶斯高斯混合先验模块，其建筑物提取效果较差，投影面积缺失严重。方法2 在方法1 的基础上加入贝叶斯高斯混合先验模块，极大提高了检测到的建筑物面积，但仍然存在建筑物边界混乱、多余检测（FP值较高）问题，这是由于贝叶斯高斯混合模型仅在独立假设条件下进行预测，未考虑到超体素的邻接关系。方法3 在方法2 的基础上加入高阶因子模块，高阶因子聚合了相邻的超体素，为实验增加了邻域信息，其建筑物提取结果更接近于真值，但建筑物边界依然不够清晰，且建筑物与植被粘连区域的FN 和FP 值较大。方法4 即本文方法，在方法3基础上加入了方差分析模块，对比前3 部分实验结果，其建筑物边界更加清晰，投影面积误差最小，这说明方差分析可以有效区分不同质的高阶因子，提高建筑物提取精度，并优化建筑物边界。此外，表4列出了消融实验中各方法的建筑物提取精度。为简化表达，本文以基于点的评价指标为例，分析各模块对建筑物提取精度的影响，基于投影面积和对象的精度变化情况与基于点的精度变化情况类似。从表4中可以看到，随着各模块的加入，各项精度逐步提升。其中，每增加一个模块，表示整体建筑物提取质量的精度评价指标Q和F1 增幅保持在10%—29%。加入所有模块后，建筑物提取精度最高，效果也最佳。因此，方差分析引导的高阶马尔可夫网络中各模块对建筑物提取均至关重要。其中，贝叶斯高斯混合先验模块提供了较为准确的样本先验状态，为后续优化提供了良好数据基础。高阶因子模块补充了样本之间的邻接关系，提高了建筑物提取结果的平整度。方差分析模块则充分考虑了样本之间的统计差异，通过P值检验，有效区分不同质的高阶因子，得到更为精确的点云建筑物提取结果，并优化了建筑物边界。

表4 消融实验的建筑物提取精度（区域2）Table 4 The accuracy of ablation studies on building extraction （Area 2）

图7 消融实验结果可视化对比Fig. 7 Visual comparisons of the ablation studies

3.6 对比分析

为充分验证本文方法的优越性和不足，本文与现有方法进行了比较。由于现有方法的实验过程所使用的数据集不尽相同，本文首先与利用Vaihingen 数据集（区域1—3）验证的8 种方法进行对比，再与利用Toronto 数据集（区域4—5）验证的5种方法进行对比，从而充分论证本文方法的优势和特点。其中，CSU（Du 等，2017）是融合点特征和网格特征的建筑物提取方法，其采用图割算法检测建筑物，并通过形态学滤波优化建筑物区域。AISeg（Chen 等，2020）是基于数据融合的建筑物提取方法，其采用自适应迭代分割（Adaptive Iterative Segmentation） 和分层分割（Hierarchical Segmentation）原理。MDFC（Maltezos 等，2019）是基于多维特征向量（Multi-Dimensional Feature Vector）的建筑物提取方法。MPP（Liu 等，2020）是基于最小割（Min-cut）算法的建筑物提取方法，其后处理过程使用上下文信息和一致性约束消除建筑物的内部异质性。HANC2（Niemeyer 等，2013）是基于CRF 框架的建筑物提取方法，其纳入多尺度特征改善提取结果。ITCR（Gerke 和Xiao，2014）是基于RF 算法的建筑物提取方法，其实验数据融合了ALS 点云信息和光谱图像信息。TUM（Wei 等，2012）是基于Adaboost 分类器的建筑物提取方法。Mar1（Mongus 等，2013）是基于差分形态学的点云建筑物提取方法。MON 和MON2（Awrangjeb 和Fraser，2014）是基于点云共面特征和邻域信息的建筑物提取方法，两种方法均使用DEM（Digital Elevation Model）分离地面点和非地面点。

表5和表6列出了不同方法在区域1—3和区域4—5 的平均精度结果。从精度对比来看，本文提出的非监督方法在面积和对象级别上精度普遍较高。在区域1—3上，本文方法的Compar、Qar、F1ar、Compobj、Qobj和F1obj值均高于其他方法，这意味着在区域1—3，本文方法更易识别出建筑物。在区域4—5，本文方法结果在以面积为量化因子的情况下，精度低于一部分方法（如MPP）。通过仔细分析，其原因在于区域4—5 内高层建筑立面点分布稀疏，这与本文先验假设——建筑物离散度较低有所出入，因此，该区域检测的建筑物立面缺失较严重。其次，对于部分顶部附件复杂的建筑物，其不符合建筑物点云由平面组成这一预设，因此这部分建筑物检测效果不佳。但是在以对象为量化因子的情况下，本文方法在区域4—5 的整体精度处于领先水平。综上所述，方差分析引导的超体素高阶马尔可夫网络，在不需要标注数据训练的情况下，可以有效提取点云建筑物，方法的创新性在理论上和实验上，均得到了验证。

表5 不同建筑物提取方法的平均精度对比（区域1-3）Table 5 The comparison between average accuracy scores of different building extraction methods （Area 1 to 3）

表6 不同建筑物提取方法的平均精度对比（区域4-5）Table 6 The comparison between average accuracy scores of different building extraction methods （Area 4 to 5）

4 结 论

本文提出一种方差分析引导的高阶超体素马尔可夫网络，并将其应用于ALS 点云建筑物提取任务，提高了非监督点云建筑物提取精度。首先，本文构建了一种由粗到精的无监督建筑物提取框架。该框架在独立假设条件下，利用贝叶斯高斯混合模型捕捉超体素类别先验，实现建筑物粗提取。接着结合高阶马尔可夫网络和方差分析为超体素邻域的相关性建模，对点云建筑物提取结果予以优化。第二，本文构建了一种基于方差分析P值检验的高阶无向图势函数计算方法，使超体素邻域相关性建模从传统的经验、几何分析，转化成依概率统计属性的精确推理过程，增强了势函数对邻接超体素之间相关性的表达能力。第三，本文将高阶因子引入到点云超体素马尔可夫网络模型，构建了一种高阶无向图网络，实现了三维超体素局部邻域相关性到全局依赖关系的拓展。并采用信念传播算法，实现最大后验状态推理，以精准识别建筑物。

实验采用国际摄影测量与遥感学会（ISPRS）提供的Vaihingen 和Toronto 机载激光雷达点云数据集，结合4种通用精度评价指标，对本文方法进行性能分析，并与现有方法进行对比。结果表明，本文方法能有效提取不同场景下的建筑物且精度较高，特别在建筑物面积较大的住宅区，取得了优异的提取效果。此外，为分析各子模块的有效性，本文进行了消融实验，从定量和定性角度出发，肯定了各子模块的积极作用。

本文方法的局限性在于尚未知如何在高阶无向图网络框架下，选择特征并确定其在网络构建过程中的权重。在后续的研究中，将聚焦该问题并寻求突破。