郝嬌嬌， 倪歡， 管海燕

南京信息工程大學(xué) 遙感與測繪工程學(xué)院, 南京 210044

1 引 言

應(yīng)三維城市建模、GIS 數(shù)據(jù)庫更新等應(yīng)用需求，建筑物提取已成為近些年的研究熱點(diǎn)（Khoshelham 等，2010；Tomljenovic 等，2015）。ALS（Airborne Laser Scanning）技術(shù)是獲取空間數(shù)據(jù)最直接的手段之一（Du等，2017），快速高效地從海量ALS 點(diǎn)云數(shù)據(jù)所包含的復(fù)雜空間信息中提取建筑物，成為了一系列地學(xué)應(yīng)用的首要任務(wù)。

目前，ALS點(diǎn)云建筑物提取可分為融合提取和直接提取兩類（鄧飛 等，2018）。第一類方法融合了圖像光譜和點(diǎn)云三維幾何信息；優(yōu)點(diǎn)是通過融合多源數(shù)據(jù)，實(shí)現(xiàn)信息互補(bǔ)，進(jìn)而提高建筑物提取精度（Gerke 和Xiao，2014；Zarea 和Mohammadzadeh，2016）。例如，點(diǎn)云與圖像紋理的特征匹配，不僅能有效區(qū)分相似高程地物的種類（Lai 等，2019），還可提高區(qū)域分析和監(jiān)測的效率（Ullo等，2020）。在融合兩類數(shù)據(jù)的基礎(chǔ)上，常使用由粗到精的提取策略。首先結(jié)合多重特征獲取初始建筑物，然后根據(jù)一定范圍內(nèi)的正射影像顏色信息優(yōu)化建筑物邊緣（杜守基 等，2018），或利用形態(tài)學(xué)運(yùn)算優(yōu)化建筑物候選區(qū)域（Chen 等，2020），以提高建筑物提取精度。此外，深度學(xué)習(xí)技術(shù)由于其在視覺數(shù)據(jù)上的諸多優(yōu)勢，被引入到該類方法。郭峰等（2020）提出一種基于Encoder-decoder 語義模型的建筑物提取方法（SegNet），將兩類數(shù)據(jù)（點(diǎn)云與圖像）嵌入到深度學(xué)習(xí)模型。有別于SegNet 的監(jiān)督學(xué)習(xí)過程，Nguyen 等（2020）提出了一種融合點(diǎn)云與圖像的大范圍建筑物無監(jiān)督提取方法，具有較高的自動化水平。然而，大范圍圖像和點(diǎn)云的自動配準(zhǔn)過程會增加建筑物提取的誤差傳播風(fēng)險(xiǎn)（Parmehr 等，2014）；同時(shí)，兩種數(shù)據(jù)分辨率和獲取時(shí)間的不同，也會產(chǎn)生幾何和光譜特征不協(xié)調(diào)問題。

第二類方法僅使用點(diǎn)云數(shù)據(jù)，包括監(jiān)督和非監(jiān)督方法。傳統(tǒng)的監(jiān)督分類方法一般會提取建筑物和其他地物的高程、紋理、法線向量、回波次數(shù)和強(qiáng)度特征，根據(jù)其中某一特征、特征組合或者衍生特征（鄧飛 等，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）等單個(gè)或多種經(jīng)典機(jī)器學(xué)習(xí)算法進(jìn)行建筑物提取。這些方法充分利用了建筑物的先驗(yàn)信息，對大量標(biāo)注樣本的依賴程度較低，但分類時(shí)未學(xué)習(xí)到更高層次語義，分類精度難以進(jìn)一步提高（Niemeyer 等，2014）。隨著深度學(xué)習(xí)技術(shù)的不斷推廣，基于卷積神經(jīng)網(wǎng)絡(luò)（Maltezos 等，2019）、深度殘差網(wǎng)絡(luò)（趙傳 等，2020）、圖幾何矩卷積（Li 等，2020）的建筑物提取方法在逐步更新和完善。這類方法提高了建筑物識別精度，但前提是需要大量標(biāo)注數(shù)據(jù)進(jìn)行模型訓(xùn)練。在數(shù)據(jù)量有限的應(yīng)用場景下，深度學(xué)習(xí)算法無法對數(shù)據(jù)的規(guī)律進(jìn)行無偏差估計(jì)。此外，訓(xùn)練過程耗時(shí)較長且依賴于高性能計(jì)算設(shè)備。在這種情況下，學(xué)者們應(yīng)該重視非監(jiān)督的、不依賴于人工標(biāo)注的識別方法。非監(jiān)督方法一般根據(jù)地物的不同屬性設(shè)計(jì)出特定規(guī)則，并引入一系列模式識別算法，實(shí)現(xiàn)建筑物與其他類別的自適應(yīng)分割。其中，逆迭代數(shù)學(xué)形態(tài)學(xué)（RIMM）算法（Cheng 等，2013）、差分形態(tài)剖面（DMPs）映射（Mongus 等，2014）、投票模型（趙傳 等，2017）、非流形理論（Gilani等，2018）、自頂向下（Top-down）策略（Huang 等，2018）、MRF（Markov Random Field；Wang 等，2020）、最小割算法（Min-cut；Liu 等，2020）、層次角感知（Widyaningrum 等，2020）等被引入到點(diǎn)云建筑物非監(jiān)督提取過程中，并均在相應(yīng)的實(shí)驗(yàn)數(shù)據(jù)上取得了較好結(jié)果。

非監(jiān)督方法目前雖精度不及監(jiān)督學(xué)習(xí)方法，但其不需要借助大量人工標(biāo)記數(shù)據(jù)進(jìn)行訓(xùn)練，對硬件設(shè)備要求不高，使用更為靈活。本研究聚焦非監(jiān)督識別方法，在統(tǒng)計(jì)學(xué)框架下，實(shí)現(xiàn)點(diǎn)云建筑物提取。相比于深度學(xué)習(xí)方法，本文方法的優(yōu)勢和特色主要體現(xiàn)在3 個(gè)方面，即（1）無人工標(biāo)注依賴，不存在過擬合問題，對符合預(yù)設(shè)規(guī)則的數(shù)據(jù)，泛化能力更強(qiáng)；（2）算法運(yùn)行時(shí)間短，不依賴于高性能計(jì)算設(shè)備，能夠更為快速地獲取建筑物識別結(jié)果；（3）從數(shù)理統(tǒng)計(jì)角度發(fā)掘點(diǎn)云空間上下文關(guān)系，理論可解釋性更強(qiáng)，所發(fā)現(xiàn)的相關(guān)關(guān)系具備可回溯性。方差分析ANOVA（Analysis of Variance）是對兩組及兩組以上樣本進(jìn)行顯著性檢驗(yàn)的統(tǒng)計(jì)方法（Pulido-González 等，2020），能夠準(zhǔn)確檢驗(yàn)樣本的組間差異。因此，ANOVA 方法為分析多維度特征之間的相關(guān)性提供了一條便捷途徑。

基于以上分析，本文提出了一種方差分析引導(dǎo)的高階馬爾可夫網(wǎng)絡(luò)，構(gòu)建了一種基于方差分析P值檢驗(yàn)的高階無向圖勢函數(shù)計(jì)算方法，實(shí)現(xiàn)點(diǎn)云超體素上下文信息建模，并精準(zhǔn)提取建筑物。

2 方差分析引導(dǎo)的高階超體素馬爾可夫網(wǎng)絡(luò)

點(diǎn)云數(shù)據(jù)量通常較大，難以用每一個(gè)點(diǎn)作為節(jié)點(diǎn)進(jìn)行相關(guān)性建模。本文將點(diǎn)云數(shù)據(jù)進(jìn)行濾波處理，以剔除地面點(diǎn)；構(gòu)建超體素結(jié)構(gòu)，建立對象級網(wǎng)絡(luò)與分析框架。其中地面點(diǎn)剔除和超體素構(gòu)建采用CSF（Zhang 等，2016）和SVLA（Ni 和Niu，2020）算法。方差分析引導(dǎo)的高階超體素馬爾可夫網(wǎng)絡(luò)建筑物提取過程如圖1所示。該方法首先以超體素為基元計(jì)算特征，并將超體素作為網(wǎng)絡(luò)節(jié)點(diǎn)。接著，使用貝葉斯高斯混合模型實(shí)現(xiàn)初始聚類，以獲取馬爾可夫網(wǎng)絡(luò)各節(jié)點(diǎn)的初始狀態(tài)。然后，檢索超體素鄰域，采用方差分析原理，結(jié)合預(yù)定義的特征，構(gòu)建高階因子，并計(jì)算節(jié)點(diǎn)、邊和高階勢函數(shù)，形成超體素高階馬爾可夫網(wǎng)絡(luò)。最后，采用信念傳播算法，實(shí)現(xiàn)高階馬爾可夫網(wǎng)絡(luò)的近似推理，以提取建筑物。

圖1 方差分析引導(dǎo)的高階超體素馬爾可夫網(wǎng)絡(luò)建筑物提取流程圖Fig. 1 The flowchart of building extraction using ANOVA guided high-order supervoxel Markov network

2.1 基于貝葉斯高斯混合模型的節(jié)點(diǎn)初始狀態(tài)捕捉

高斯混合模型（Chen 等，2019）采用多個(gè)高斯分布線性疊加的方式來表示樣本分布情況，其收斂性受到初始假設(shè)影響，并且其常用的期望最大化EM（Expectation Maximum）解法經(jīng)常出現(xiàn)奇異值。貝葉斯高斯混合模型可以解決以上問題，其采用變分推理思想，使用形式簡單的分布來逼近高斯混合模型，并得到一種局部最優(yōu)且具有確定解的近似后驗(yàn)分布。

假設(shè)超體素集合為S={S1，S2，…，SM}，其相應(yīng)的特征集合為X={x1，x2，…，xM}，并且馬爾可夫網(wǎng)絡(luò)的節(jié)點(diǎn)狀態(tài)空間是二值的。貝葉斯高斯混合模型的概率密度函數(shù)定義如下：

式中，

是一個(gè)多元高斯分布，C是類別數(shù)，αc為混合系數(shù)，μc為第c個(gè)高斯分布的均值向量，Σc代表第c個(gè)高斯分布的協(xié)方差矩陣。該模型的解對應(yīng)參數(shù)值αc，μc，Σc。其求解通過優(yōu)化證據(jù)下界ELBO（Evidence Lower BOund）過程完成。更新后的參數(shù)帶入式（2），求出節(jié)點(diǎn)取得各個(gè)狀態(tài)的概率，即每個(gè)超體素從屬于各狀態(tài)的概率分布表，其對應(yīng)馬爾可夫網(wǎng)絡(luò)中節(jié)點(diǎn)的初始狀態(tài)。在本文的研究中，該初始狀態(tài)給出了每一個(gè)節(jié)點(diǎn)從屬于建筑物和非建筑物的初始概率。

2.2 超體素鄰域檢索

超體素鄰域Ni指第i個(gè)超體素Si周圍的最鄰近區(qū)域，其由多個(gè)與Si相鄰的超體素，k=1，…，K構(gòu)成。本文采用由點(diǎn)鄰域到超體素鄰域拓展的方式，確定超體素Si的鄰域信息。首先采用KD 樹（Buitinck 等，2013）檢索Si中每一個(gè)點(diǎn)的鄰接點(diǎn)。若鄰接點(diǎn)的集合中含有其他超體素中的點(diǎn)，則認(rèn)為該鄰接點(diǎn)從屬的超體素位于Si鄰域，此時(shí)，記錄該鄰接點(diǎn)對之間的距離，并作為與Si之間的距離dk。以此類推，確定與Si相鄰的所有超體素，并將其集合作為初始鄰域。如圖2所示，從檢索點(diǎn)的鄰接關(guān)系出發(fā)，擴(kuò)展到超體素的鄰接關(guān)系，最終得到超體素A的鄰域包含超體素B-H。此外，本文采用與點(diǎn)鄰域檢索類似的鄰域大小限制，即根據(jù)每一個(gè)鄰接體素到Si之間的距離dk，取距離Si最近的Kmax個(gè)超體素構(gòu)成Si的最終鄰域Ni。引入該閾值Kmax，可以增強(qiáng)本文方法在不同點(diǎn)密度和超體素尺寸條件下的調(diào)節(jié)能力。例如，若超體素尺寸較大，將設(shè)置較小的Kmax，反之，則設(shè)置較大。

圖2 超體素鄰域示意圖Fig. 2 The schematic diagram of the neighborhood of a supervoxel

Kmax個(gè)鄰接超體素和當(dāng)前超體素Si共同構(gòu)成了高階馬爾可夫網(wǎng)絡(luò)中每一個(gè)高階因子，該結(jié)構(gòu)填補(bǔ)了點(diǎn)云上下文信息的缺失。

2.3 基于方差分析原理的高階因子構(gòu)建

假設(shè)預(yù)先計(jì)算了T個(gè)維度的超體素特征X={x1，x2，…，xM}，其中xi={xi1，xi2，…，xiT}T。依方差分析原理，對X進(jìn)行線性最小二乘方差分析檢驗(yàn)，其觀測值為超體素各維度特征向量。首先計(jì)算各維度（組）觀測值均值xˉi以及全部觀測值的總體均值：

式中，Kmax是當(dāng)前鄰域內(nèi)的超體素個(gè)數(shù)，Kmax+ 1表示同時(shí)考慮了當(dāng)前超體素及其鄰域，則Kmax相當(dāng)于該過程的自由度。n表示全部觀測值個(gè)數(shù)。接著，根據(jù)均值計(jì)算各誤差平方和，包括全部觀測值與總均值的誤差平方和SST（Sum of Squares for Total），組間誤差平方和SSA（Sum of Squares for Factor A）以及組內(nèi)誤差平方和SSE（Sum of Squares for Error）：

為了得到無偏估計(jì)結(jié)果，需要用各平方和除以自由度（即Kmax），該結(jié)果在方差分析中被稱作方差。相應(yīng)的，組間方差MSA（Mean Squares for Factor A）和組內(nèi)方差MSE（Mean Squares for Error）定義為

進(jìn)而，檢驗(yàn)統(tǒng)計(jì)量Fv表示為

最后，根據(jù)檢驗(yàn)統(tǒng)計(jì)量Fv，利用多元F分布檢驗(yàn)，找出Fv的伴隨P值。該統(tǒng)計(jì)量Fv和P值可以有效反映超體素之間的統(tǒng)計(jì)相關(guān)性，將用以實(shí)現(xiàn)高階因子構(gòu)建，并作為勢函數(shù)計(jì)算的依據(jù)。

P值的大小反應(yīng)了該鄰域內(nèi)所有超體素特征服從同一概率分布的可能性大小，若其較大，本文將該鄰域的超體素結(jié)構(gòu)均納入到同一個(gè)因子中。因該因子所包含的節(jié)點(diǎn)（即超體素）多于2個(gè)，所以稱作高階因子。被劃分到同一高階因子的節(jié)點(diǎn)之間，使用無向邊相連。

本文在經(jīng)典的成對馬爾可夫網(wǎng)絡(luò)基礎(chǔ)上，增加高階因子，構(gòu)成高階馬爾可夫網(wǎng)絡(luò)，以提高網(wǎng)絡(luò)上下文信息表達(dá)能力。此時(shí)，高階馬爾可夫網(wǎng)絡(luò)的聯(lián)合概率分布形式如下：

式中，X={x1，x2，…，xM}是網(wǎng)絡(luò)中所有節(jié)點(diǎn)構(gòu)成的集合，Z*是配分函數(shù)，ψF(·)、φF(·)、φF(·)分別是節(jié)點(diǎn)、邊、高階勢函數(shù)，E是邊集合，F(xiàn)SH是高階因子集合，Hl是第l個(gè)高階因子。接著，使用信念傳播算法對式（11）的聯(lián)合概率分布進(jìn)行最大后驗(yàn)MAP（Maximum A Posterior）狀態(tài)推理，得到最佳狀態(tài)分布：

2.4 P值引導(dǎo)的勢函數(shù)計(jì)算

2.4.1 節(jié)點(diǎn)勢函數(shù)

經(jīng)過貝葉斯高斯混合模型的獨(dú)立預(yù)測過程，各節(jié)點(diǎn)（超體素）從屬于各類別的概率p(xm∈Lc)被計(jì)算出來，其中Lc是第c個(gè)類別對應(yīng)的節(jié)點(diǎn)集合。則節(jié)點(diǎn)勢函數(shù)定義如下：

式中，ψF(xm∈Lc)對應(yīng)于未歸一化的概率分布表，其事件空間包含C個(gè)事件（即類別）。馬爾可夫網(wǎng)絡(luò)的歸一化概率，最后由配分函數(shù)Z*歸一化求得，本文采用信念傳播算法求取Z*。

2.4.2 邊勢函數(shù)

邊勢函數(shù)反映了相鄰節(jié)點(diǎn)所對應(yīng)超體素之間的特征相關(guān)性，相關(guān)性越高，從屬于同一類別的概率越大。本文采用前述方差分析所得到的P值來確定該先驗(yàn)相關(guān)性。具體地，P值越大，說明該鄰域內(nèi)超體素之間的同質(zhì)性越高。因此，若P值大于某一閾值，則該鄰域內(nèi)的所有超體素從屬于同一類別。若P值較小，則需要對鄰域內(nèi)的每一對超體素進(jìn)行多重對比檢驗(yàn)（Test for multiple comparisons）。本文采用Tukey HSD（Tukey’s Honestly Significant Difference）檢驗(yàn)方法，計(jì)算每一對超體素之間的HSD 檢驗(yàn)Ptukey值，以確定每一對超體素之間的相關(guān)性大小。綜上所述，每一條邊上的先驗(yàn)相關(guān)信息定義如下：

式中，TRP是P值的顯著性閾值，研究一般認(rèn)為，TRP值小于0.05為有統(tǒng)計(jì)學(xué)差異?；诖?，本文取TRP= 0.05。Ptukey是Tukey HSD 檢驗(yàn)方法計(jì)算得到的P值。

P值雖然刻畫了不同超體素服從同一分布的聯(lián)合概率，但其仍然無法在局部區(qū)域上衡量邊上相關(guān)性的大小。本文利用一種局部相關(guān)性先驗(yàn)概率加權(quán)該聯(lián)合分布，以增強(qiáng)其表達(dá)能力。假設(shè)當(dāng)前待確定的邊為E(i，j)，則其鄰域?yàn)槠溥B接的兩個(gè)節(jié)點(diǎn)鄰域集合的并集，即S(i&j) =S(i) ∪S(j)。因各節(jié)點(diǎn)存在多個(gè)狀態(tài)（對應(yīng)于分類問題的類別數(shù)），本文為每一類別計(jì)算了一個(gè)先驗(yàn)概率以衡量各對鄰接超體素的局部相關(guān)性。

式中，I(xm∈Lc)是指示函數(shù)：

進(jìn)而，邊勢函數(shù)定義如下：

式（18）表達(dá)xi和xj具有相同狀態(tài)時(shí)的勢函數(shù)，式（19）表達(dá)xi和xj具有不同狀態(tài)時(shí)的勢函數(shù)。ratee是調(diào)節(jié)系數(shù)。

2.4.3 高階勢函數(shù)

高階勢函數(shù)中包含的變量對應(yīng)多個(gè)節(jié)點(diǎn)（即超體素），參考節(jié)點(diǎn)和邊勢函數(shù)的定義，其具體形式定義如下：

式中，L(x1) =L(x2) = … =L(xKmax+1)對應(yīng)于各節(jié)點(diǎn)取得相同狀態(tài)的情況。rateh是調(diào)節(jié)系數(shù)。其中P是方差分析計(jì)算得到的P值。

2.5 基于信念傳播算法的高階馬爾可夫網(wǎng)絡(luò)推理

本文使用信念傳播算法（Hazan 和Shashua，2010）對上述高階馬爾可夫網(wǎng)絡(luò)進(jìn)行最大后驗(yàn)狀態(tài)推理。信念傳播算法首先將高階馬爾可夫網(wǎng)絡(luò)轉(zhuǎn)化成因子圖GF，其由因子節(jié)點(diǎn)和變量節(jié)點(diǎn)構(gòu)成，因子節(jié)點(diǎn)和變量節(jié)點(diǎn)之間由無向邊相連。信念傳播算法通過在因子和變量節(jié)點(diǎn)之間迭代傳遞消息的方式，不斷更新因子和勢函數(shù)直到收斂。

信念傳播算法是一種求解圖模型推理問題的高效算法，與基于能量最小化原理的圖割法相比，靈活度更高。信念傳播算法可以對任意勢函數(shù)進(jìn)行優(yōu)化，包括一些圖割算法所不能處理的非規(guī)則（不滿足次模性的）勢函數(shù)（Koller 和Friedman，2009；Le^-Huu，2019；徐勝軍 等，2013）。由于信念傳播算法在機(jī)器學(xué)習(xí)任務(wù)中被廣泛使用，本文不對其原理進(jìn)行冗余敘述。

2.6 面向建筑物提取的特征計(jì)算

本文旨在利用高階馬爾可夫網(wǎng)絡(luò)提取建筑物，因此面向建筑物提取相應(yīng)特征。本文所提取的特征包括：（1）離散度，（2）法線向量傾角，（3）法線方向方差，（4）曲率和（5）維度特征。

（1）離散度。離散度衡量超體素所包含點(diǎn)的空間分布離散程度，其計(jì)算超體素內(nèi)各點(diǎn)pi=(xi，yi，zi)到中心點(diǎn)(xm，ym，zm)距離的平均值作為離散度。即：

式中，N是超體素中點(diǎn)的總數(shù)。

（2）法線向量傾角。本文以超體素為單元，利用各超體素所包含的點(diǎn)擬合空間平面，進(jìn)而計(jì)算超體素法線向量。本文以角度的余弦值來表達(dá)法線向量傾角。設(shè)第m個(gè)超體素內(nèi)的點(diǎn)集是Sm={pi|i= 1，2，…，N}，pi=(xi，yi，zi)，那么超體素的協(xié)方差矩陣為

（3）法線方向方差。本文引用文獻(xiàn) （Du 等，2017）的定義，檢索超體素內(nèi)部每一個(gè)點(diǎn)的鄰域點(diǎn)集，利用式（22）和（23）計(jì)算協(xié)方差矩陣和法線向量傾角。再計(jì)算該超體素內(nèi)部所有點(diǎn)的法線向量傾角方差σm，用以描述超體素的平整度。

（4）曲率。對平整度的衡量，通常還可以使用曲率，其基于式（22）得到的協(xié)方差矩陣來計(jì)算。即計(jì)算Cm的3 個(gè)特征值λ0≤λ1≤λ2，則曲率為

（5）維度特征。維度特征（Demantké等，2011）使用廣泛，且效果較好，本文引用其中兩種，具體如下：

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

3 實(shí)驗(yàn)分析與對比

3.1 實(shí)驗(yàn)數(shù)據(jù)

本文使用國際攝影測量與遙感學(xué)會（ISPRS）提供的參考數(shù)據(jù)集，包括Vaihingen 和Toronto 兩組機(jī)載激光雷達(dá)點(diǎn)云（Liu 等，2020）。為了便于與其他方法對比，本文對兩個(gè)數(shù)據(jù)集中的各標(biāo)注區(qū)域采用類似的命名方式，即區(qū)域1—3（來自Vaihingen 數(shù)據(jù)集）和區(qū)域4—5（來自Toronto 數(shù)據(jù)集）。區(qū)域1—3的平均點(diǎn)密度約為4—7個(gè)點(diǎn)/m2，位于住宅區(qū)。區(qū)域4—5的平均點(diǎn)密度約為6個(gè)點(diǎn)/m2，位于商業(yè)區(qū)。圖3展示了各區(qū)域按高程著色的點(diǎn)云場景。

圖3 高程著色的測試數(shù)據(jù)Fig. 3 Testing data colored by elevation

區(qū)域1 包括37 棟建筑，主要由密集的古老建筑組成，形狀相對復(fù)雜。區(qū)域2 包括14 棟建筑，主要由屋頂水平的高層住宅組成。區(qū)域3 包括56棟建筑，主要由獨(dú)立建筑組成，屋頂結(jié)構(gòu)簡單，沿路有植被。區(qū)域4 包括58 棟建筑，由不同高度建筑組成，建筑物結(jié)構(gòu)復(fù)雜，頂部有附著物。區(qū)域5 包括38 棟建筑，建筑物結(jié)構(gòu)多樣，形狀復(fù)雜，頂部帶有附著物；其中包括一個(gè)顯著高于其他地物的高層建筑。

3.2 精度評價(jià)指標(biāo)

本文采用ISPRS 指標(biāo)體系（Rottensteiner 等，2014）衡量結(jié)果質(zhì)量，包括Completeness、Correctness、Quality和F1指數(shù)4個(gè)評價(jià)指標(biāo)，計(jì)算方法如下：

以上4 種評價(jià)指標(biāo)，在采用不同量化因子時(shí)，含義存在細(xì)微差別。其中，基于面積和基于對象的指標(biāo)體系最為常用?；诿娣e的評價(jià)指標(biāo)記為Compar、Corrar、Qar、F1ar，以投影到二維平面的面積作為量化因子，該面積通常使用像素個(gè)數(shù)衡量，精度與正確分類的區(qū)域相關(guān)。此時(shí)，TP表示正確檢測的建筑物面積，F(xiàn)P表示錯(cuò)誤檢測的建筑物面積，F(xiàn)N表示遺漏的建筑物面積?；趯ο蟮脑u價(jià)指標(biāo)記為Compobj、Corrobj、Qobj、F1obj，以建筑物個(gè)數(shù)作為量化因子，若提取對象與參考對象重疊度大于50%，則認(rèn)為識別到該對象。此時(shí)，TP表示正確檢測的建筑物個(gè)數(shù)，F(xiàn)P表示錯(cuò)誤檢測的建筑物個(gè)數(shù)，F(xiàn)N表示遺漏的建筑物個(gè)數(shù)。此外，Rottensteiner 等（2014）還提供了用于大型建筑物檢測的評價(jià)指標(biāo)，但是面積大于50 m 的建筑物提取精度普遍較高，可比性較弱，因此，不單獨(dú)對大面積建筑物進(jìn)行精度評價(jià)。為更細(xì)致地評估建筑物提取性能，本文增加了基于點(diǎn)的評價(jià)指標(biāo)，即逐點(diǎn)根據(jù)式（27）—式（30）計(jì)算精度，記為Comppoi、Corrpoi、Qpoi、F1poi。此時(shí)，TP表示正確檢測的建筑物點(diǎn)數(shù)，F(xiàn)P表示錯(cuò)誤檢測的建筑物點(diǎn)數(shù)，F(xiàn)N表示遺漏的建筑物點(diǎn)數(shù)。

3.3 參數(shù)設(shè)置

本文使用Ni 和Niu（2020）提出的SVLA 算法構(gòu)建超體素，該算法可以獲得緊致的超體素簇，并通過調(diào)整種子分辨率（Rseed）來控制超體素大小。體素的大小是影響建筑物提取結(jié)果的重要因素，若構(gòu)建的點(diǎn)云超體素較小，則意味著點(diǎn)云數(shù)據(jù)中所含超體素過多，這會增加算法的計(jì)算復(fù)雜度和運(yùn)行時(shí)間，降低建筑物提取效率。反之，若構(gòu)建的超體素較大，雖充分顧及了上下文信息，卻忽略了局部細(xì)節(jié)特征，進(jìn)而影響建筑物的提取精度。表1 記錄了區(qū)域1 中Rseed設(shè)置與建筑物提取精度的相互影響關(guān)系。為了簡化表達(dá)，本文以基于點(diǎn)的評價(jià)指標(biāo)為例進(jìn)行分析，其他基于面積和對象的精度變化情況與基于點(diǎn)的精度變化情況類似。從表中可以看到，當(dāng)0.9≤Rseed≤1.3 時(shí)，建筑物提取精度逐漸上升，Rseed取1.3 時(shí)，各項(xiàng)精度評價(jià)指標(biāo)達(dá)到最大值。當(dāng)Rseed＞1.3 時(shí)，隨著Rseed的增大，建筑物提取的精度逐漸下降。因此，構(gòu)建合適大小的超體素至關(guān)重要。本文經(jīng)過實(shí)驗(yàn)對比，最終將區(qū)域1—3的種子分辨率設(shè)置為1.3，區(qū)域4—5的種子分辨率設(shè)置為2.0。此外，在輸入數(shù)據(jù)中不可避免地會存在一部分較小的超體素，包括噪聲點(diǎn)和小型地物，一定程度上影響了結(jié)果的精度。因此，在超體素分割完成后，需要剔除較小的超體素，即超體素內(nèi)的點(diǎn)數(shù)小于t時(shí)，刪除該超體素。如表2 所示，區(qū)域4 和區(qū)域5 高層建筑物居多，建筑物體積較大，因此，設(shè)置t的值也最大。

表1 Rseed對建筑物提取精度的影響（區(qū)域1）Table 1 The relationship between Rseed and accuracy （Area 1）

鄰域檢索參數(shù)Kmax（2.2 節(jié)）控制鄰域內(nèi)超體素的個(gè)數(shù)，其對應(yīng)了高階馬爾可夫網(wǎng)絡(luò)中高階因子所包含的節(jié)點(diǎn)個(gè)數(shù)。由于對區(qū)域1—3 的超體素分割較為細(xì)碎，鄰域內(nèi)包含的超體素較多，因此，設(shè)置Kmax為9。相應(yīng)的，區(qū)域4—5 的體素塊較大，相鄰區(qū)域包含的超體素較少，設(shè)置Kmax為7。在2.4 節(jié)計(jì)算邊和團(tuán)勢函數(shù)的過程中，經(jīng)過粗對比，調(diào)節(jié)系數(shù)ratee和rateh在0.1—1.5時(shí)，建筑物提取精度較高。圖4表示不同ratee和rateh取值對應(yīng)的平均F1poi值。由圖4 可知，在區(qū)域1—3 上，ratee和rateh的值為0.3 和1.0 時(shí)，平均F1poi值最高；在區(qū)域4—5 上，ratee和rateh的值為0.5 和1.5 時(shí)，平均F1poi值最高。各區(qū)域參數(shù)設(shè)置見表2。

表2 實(shí)驗(yàn)參數(shù)設(shè)置Table 2 Parameter settings of this experiment

圖4 ratee和rateh的參數(shù)分析Fig. 4 Parameter analysis of ratee and rateh

3.4 實(shí)驗(yàn)結(jié)果

圖5展示了本文方法在不同區(qū)域的建筑物提取結(jié)果。從圖中可以看到，場景中的建筑物提取完整，屋頂內(nèi)部無缺失，建筑物邊界清晰。大型建筑物（區(qū)域4、區(qū)域5）內(nèi)部無明顯噪聲。植被覆蓋密集區(qū)域（區(qū)域2）建筑物邊界明顯，建筑物完整。小型建筑物（區(qū)域3）提取完整，周圍無明顯噪聲，建筑物與籬笆、道路的邊界明顯。圖6為本文建筑物提取的投影面積可視化結(jié)果情況，該結(jié)果對應(yīng)了精度評價(jià)指標(biāo)Compar、Corrar、Qar和F1ar；其中黃色、紅色、藍(lán)色分別代表TP、FN和FP，可以看到，本文方法正確提取了大部分建筑物。區(qū)域1 和2 中FN的面積最小，說明此區(qū)域建筑物提取最為完整。然而，提取結(jié)果中還是存在少量誤差，多余提取和提取不完整的誤差分別為FP誤差和FN誤差。其中FP誤差主要來源是場景邊緣的樹木點(diǎn)（區(qū)域1 藍(lán)框）、離散點(diǎn)（區(qū)域3、4 藍(lán)框）以及低矮灌木點(diǎn)（區(qū)域2藍(lán)框），F(xiàn)N誤差主要存在于結(jié)構(gòu)復(fù)雜的建筑物區(qū)域，其中包括部分建筑物頂部（圖6 的綠框）和部分建筑物立面（圖6黑框）。這是因?yàn)榻ㄖ镯敳扛郊^復(fù)雜的區(qū)域，不符合預(yù)設(shè)建筑物表面由平面構(gòu)成的前提；此外，部分建筑物立面點(diǎn)分布稀疏，與預(yù)設(shè)中建筑物點(diǎn)密度較高有所出入。另一方面，方差分析的引入使得不同超體素特征差異更為明顯，顯著提高了非監(jiān)督識別的整體精度。但與此同時(shí)，方差分析的引入使本文方法的識別過程對超體素特征差異的依賴度更高，進(jìn)一步增加了不符合預(yù)設(shè)物體被誤判的概率。因此，如何更好地進(jìn)行特征選擇和特征權(quán)重設(shè)置，是后續(xù)研究的重點(diǎn)。表3列出了這兩個(gè)數(shù)據(jù)集上建筑物提取結(jié)果的各項(xiàng)精度值。可以看到，所有精度值都保持在80%以上，區(qū)域1—3的精度優(yōu)于區(qū)域4—5，說明本文方法對住宅區(qū)的提取效果更佳。綜上所述，本文方法能有效提取不同場景下的建筑物且精度較高，對于低層建筑主導(dǎo)的住宅區(qū)提取效果最佳。

表3 本文方法建筑物提取結(jié)果精度情況Table 3 The accuracy values of building extraction results produced by our method

圖5 基于高階超體素馬爾可夫網(wǎng)絡(luò)的建筑物提取結(jié)果Fig. 5 Building extraction results based on the high-order supervoxel Markov network

圖6 建筑物提取投影面積可視化結(jié)果Fig. 6 Visualization of building extraction results at the per-area level

3.5 消融分析

本文提出方差分析引導(dǎo)的高階馬爾可夫網(wǎng)絡(luò)，并將其應(yīng)用于點(diǎn)云建筑物提取，該網(wǎng)絡(luò)可以分解為成對馬爾可夫網(wǎng)絡(luò)、高階因子、貝葉斯高斯混合先驗(yàn)和方差分析4個(gè)模塊，其中成對馬爾可夫網(wǎng)絡(luò)結(jié)構(gòu)是本文方法框架的基礎(chǔ)。本節(jié)選取區(qū)域2數(shù)據(jù)進(jìn)行消融實(shí)驗(yàn)，通過對比建筑物提取的可視化結(jié)果和精度，驗(yàn)證各模塊的有效性。

圖7展示了刪除特定模塊后建筑物提取結(jié)果和投影面積可視化結(jié)果，其中為了便于展示結(jié)果細(xì)節(jié)，圖中選取區(qū)域2的一個(gè)局部區(qū)域進(jìn)行分析。方法1采用成對馬爾可夫網(wǎng)絡(luò)實(shí)現(xiàn)建筑物提取，刪除了本文方法中的方差分析、高階因子和貝葉斯高斯混合先驗(yàn)?zāi)K，其建筑物提取效果較差，投影面積缺失嚴(yán)重。方法2 在方法1 的基礎(chǔ)上加入貝葉斯高斯混合先驗(yàn)?zāi)K，極大提高了檢測到的建筑物面積，但仍然存在建筑物邊界混亂、多余檢測（FP值較高）問題，這是由于貝葉斯高斯混合模型僅在獨(dú)立假設(shè)條件下進(jìn)行預(yù)測，未考慮到超體素的鄰接關(guān)系。方法3 在方法2 的基礎(chǔ)上加入高階因子模塊，高階因子聚合了相鄰的超體素，為實(shí)驗(yàn)增加了鄰域信息，其建筑物提取結(jié)果更接近于真值，但建筑物邊界依然不夠清晰，且建筑物與植被粘連區(qū)域的FN 和FP 值較大。方法4 即本文方法，在方法3基礎(chǔ)上加入了方差分析模塊，對比前3 部分實(shí)驗(yàn)結(jié)果，其建筑物邊界更加清晰，投影面積誤差最小，這說明方差分析可以有效區(qū)分不同質(zhì)的高階因子，提高建筑物提取精度，并優(yōu)化建筑物邊界。此外，表4列出了消融實(shí)驗(yàn)中各方法的建筑物提取精度。為簡化表達(dá)，本文以基于點(diǎn)的評價(jià)指標(biāo)為例，分析各模塊對建筑物提取精度的影響，基于投影面積和對象的精度變化情況與基于點(diǎn)的精度變化情況類似。從表4中可以看到，隨著各模塊的加入，各項(xiàng)精度逐步提升。其中，每增加一個(gè)模塊，表示整體建筑物提取質(zhì)量的精度評價(jià)指標(biāo)Q和F1 增幅保持在10%—29%。加入所有模塊后，建筑物提取精度最高，效果也最佳。因此，方差分析引導(dǎo)的高階馬爾可夫網(wǎng)絡(luò)中各模塊對建筑物提取均至關(guān)重要。其中，貝葉斯高斯混合先驗(yàn)?zāi)K提供了較為準(zhǔn)確的樣本先驗(yàn)狀態(tài)，為后續(xù)優(yōu)化提供了良好數(shù)據(jù)基礎(chǔ)。高階因子模塊補(bǔ)充了樣本之間的鄰接關(guān)系，提高了建筑物提取結(jié)果的平整度。方差分析模塊則充分考慮了樣本之間的統(tǒng)計(jì)差異，通過P值檢驗(yàn)，有效區(qū)分不同質(zhì)的高階因子，得到更為精確的點(diǎn)云建筑物提取結(jié)果，并優(yōu)化了建筑物邊界。

表4 消融實(shí)驗(yàn)的建筑物提取精度（區(qū)域2）Table 4 The accuracy of ablation studies on building extraction （Area 2）

圖7 消融實(shí)驗(yàn)結(jié)果可視化對比Fig. 7 Visual comparisons of the ablation studies

3.6 對比分析

為充分驗(yàn)證本文方法的優(yōu)越性和不足，本文與現(xiàn)有方法進(jìn)行了比較。由于現(xiàn)有方法的實(shí)驗(yàn)過程所使用的數(shù)據(jù)集不盡相同，本文首先與利用Vaihingen 數(shù)據(jù)集（區(qū)域1—3）驗(yàn)證的8 種方法進(jìn)行對比，再與利用Toronto 數(shù)據(jù)集（區(qū)域4—5）驗(yàn)證的5種方法進(jìn)行對比，從而充分論證本文方法的優(yōu)勢和特點(diǎn)。其中，CSU（Du 等，2017）是融合點(diǎn)特征和網(wǎng)格特征的建筑物提取方法，其采用圖割算法檢測建筑物，并通過形態(tài)學(xué)濾波優(yōu)化建筑物區(qū)域。AISeg（Chen 等，2020）是基于數(shù)據(jù)融合的建筑物提取方法，其采用自適應(yīng)迭代分割（Adaptive Iterative Segmentation） 和分層分割（Hierarchical Segmentation）原理。MDFC（Maltezos 等，2019）是基于多維特征向量（Multi-Dimensional Feature Vector）的建筑物提取方法。MPP（Liu 等，2020）是基于最小割（Min-cut）算法的建筑物提取方法，其后處理過程使用上下文信息和一致性約束消除建筑物的內(nèi)部異質(zhì)性。HANC2（Niemeyer 等，2013）是基于CRF 框架的建筑物提取方法，其納入多尺度特征改善提取結(jié)果。ITCR（Gerke 和Xiao，2014）是基于RF 算法的建筑物提取方法，其實(shí)驗(yàn)數(shù)據(jù)融合了ALS 點(diǎn)云信息和光譜圖像信息。TUM（Wei 等，2012）是基于Adaboost 分類器的建筑物提取方法。Mar1（Mongus 等，2013）是基于差分形態(tài)學(xué)的點(diǎn)云建筑物提取方法。MON 和MON2（Awrangjeb 和Fraser，2014）是基于點(diǎn)云共面特征和鄰域信息的建筑物提取方法，兩種方法均使用DEM（Digital Elevation Model）分離地面點(diǎn)和非地面點(diǎn)。

表5和表6列出了不同方法在區(qū)域1—3和區(qū)域4—5 的平均精度結(jié)果。從精度對比來看，本文提出的非監(jiān)督方法在面積和對象級別上精度普遍較高。在區(qū)域1—3上，本文方法的Compar、Qar、F1ar、Compobj、Qobj和F1obj值均高于其他方法，這意味著在區(qū)域1—3，本文方法更易識別出建筑物。在區(qū)域4—5，本文方法結(jié)果在以面積為量化因子的情況下，精度低于一部分方法（如MPP）。通過仔細(xì)分析，其原因在于區(qū)域4—5 內(nèi)高層建筑立面點(diǎn)分布稀疏，這與本文先驗(yàn)假設(shè)——建筑物離散度較低有所出入，因此，該區(qū)域檢測的建筑物立面缺失較嚴(yán)重。其次，對于部分頂部附件復(fù)雜的建筑物，其不符合建筑物點(diǎn)云由平面組成這一預(yù)設(shè)，因此這部分建筑物檢測效果不佳。但是在以對象為量化因子的情況下，本文方法在區(qū)域4—5 的整體精度處于領(lǐng)先水平。綜上所述，方差分析引導(dǎo)的超體素高階馬爾可夫網(wǎng)絡(luò)，在不需要標(biāo)注數(shù)據(jù)訓(xùn)練的情況下，可以有效提取點(diǎn)云建筑物，方法的創(chuàng)新性在理論上和實(shí)驗(yàn)上，均得到了驗(yàn)證。

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

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

4 結(jié) 論

本文提出一種方差分析引導(dǎo)的高階超體素馬爾可夫網(wǎng)絡(luò)，并將其應(yīng)用于ALS 點(diǎn)云建筑物提取任務(wù)，提高了非監(jiān)督點(diǎn)云建筑物提取精度。首先，本文構(gòu)建了一種由粗到精的無監(jiān)督建筑物提取框架。該框架在獨(dú)立假設(shè)條件下，利用貝葉斯高斯混合模型捕捉超體素類別先驗(yàn)，實(shí)現(xiàn)建筑物粗提取。接著結(jié)合高階馬爾可夫網(wǎng)絡(luò)和方差分析為超體素鄰域的相關(guān)性建模，對點(diǎn)云建筑物提取結(jié)果予以優(yōu)化。第二，本文構(gòu)建了一種基于方差分析P值檢驗(yàn)的高階無向圖勢函數(shù)計(jì)算方法，使超體素鄰域相關(guān)性建模從傳統(tǒng)的經(jīng)驗(yàn)、幾何分析，轉(zhuǎn)化成依概率統(tǒng)計(jì)屬性的精確推理過程，增強(qiáng)了勢函數(shù)對鄰接超體素之間相關(guān)性的表達(dá)能力。第三，本文將高階因子引入到點(diǎn)云超體素馬爾可夫網(wǎng)絡(luò)模型，構(gòu)建了一種高階無向圖網(wǎng)絡(luò)，實(shí)現(xiàn)了三維超體素局部鄰域相關(guān)性到全局依賴關(guān)系的拓展。并采用信念傳播算法，實(shí)現(xiàn)最大后驗(yàn)狀態(tài)推理，以精準(zhǔn)識別建筑物。

實(shí)驗(yàn)采用國際攝影測量與遙感學(xué)會（ISPRS）提供的Vaihingen 和Toronto 機(jī)載激光雷達(dá)點(diǎn)云數(shù)據(jù)集，結(jié)合4種通用精度評價(jià)指標(biāo)，對本文方法進(jìn)行性能分析，并與現(xiàn)有方法進(jìn)行對比。結(jié)果表明，本文方法能有效提取不同場景下的建筑物且精度較高，特別在建筑物面積較大的住宅區(qū)，取得了優(yōu)異的提取效果。此外，為分析各子模塊的有效性，本文進(jìn)行了消融實(shí)驗(yàn)，從定量和定性角度出發(fā)，肯定了各子模塊的積極作用。

本文方法的局限性在于尚未知如何在高階無向圖網(wǎng)絡(luò)框架下，選擇特征并確定其在網(wǎng)絡(luò)構(gòu)建過程中的權(quán)重。在后續(xù)的研究中，將聚焦該問題并尋求突破。