劉昌振,馬 威,馬 紅,魏世軒
1. 重慶市測繪科學(xué)技術(shù)研究院,重慶 401121; 2. 自然資源部智能城市時空信息與裝備工程技術(shù)創(chuàng)新中心,重慶 401121; 3. 自然資源部國土空間規(guī)劃監(jiān)測評估預(yù)警重點實驗室,重慶 401147; 4. 重慶交通大學(xué)智慧城市學(xué)院,重慶 400074
建筑物作為構(gòu)筑城市的主體,是重要的地理空間要素。利用高分辨率遙感影像[1]、LiDAR點云[2]等各類數(shù)據(jù)提取建筑物輪廓是攝影測量與遙感領(lǐng)域的研究熱點,也是地圖制圖[3]、智慧城市建設(shè)[4]等應(yīng)用的重要數(shù)據(jù)源。然而,在原始數(shù)據(jù)獲取和建筑輪廓提取過程中往往存在變形和誤差,導(dǎo)致其幾何特征無法滿足相應(yīng)的規(guī)則和應(yīng)用要求。因此,在建筑物輪廓提取后,需要進一步對其幾何進行規(guī)則化處理,構(gòu)建其規(guī)則特征[5-6]。
目前,規(guī)則化方法主要分為主方向法和最小二乘法。主方向法以計算出的建筑物主方向為約束,將滿足條件的輪廓邊進行旋轉(zhuǎn),使其和主方向垂直或平行來實現(xiàn)直角化。常用的建筑物主方向計算方法有最小外接矩形法[7-10]、統(tǒng)計加權(quán)法[11-13]和墻均值法[14-15]。其中,最小外接矩形法是通過計算建筑物的最小外接矩形,并將其作為建筑物的方向。如文獻[7]使用最小外接矩形估計概略主方向,然后利用Hough變換檢測主方向,從而提取建筑物輪廓邊。文獻[10]以最小外接矩形為主方向,通過柵格填充的方式實現(xiàn)規(guī)則化。統(tǒng)計加權(quán)法是先確定一系列候選方向,然后計算建筑物各邊對每個候選方向的貢獻值,取貢獻值最大的方向作為建筑物的方向。文獻[2]通過加權(quán)平均值計算建筑物主方向,將輪廓邊分為平行和垂直兩類,以輪廓邊中心點為軸心,旋轉(zhuǎn)輪廓邊實現(xiàn)規(guī)則化。墻均值法是以建筑物輪廓每條邊方向與權(quán)重(即長度)乘積之和的平均值作為建筑物的方向。文獻[16]使用墻均值法加強開放街道地圖中建筑物的直角特征。另外,文獻[17]結(jié)合Radon變換和主軸分析確定建筑物的主方向,然后將輪廓邊分類并進行精確定位實現(xiàn)規(guī)則化建筑物輪廓。文獻[18]使用主方向法優(yōu)化遙感影像提取的建筑物輪廓。基于主方向法的建筑輪廓提取算法相對簡單,且計算效率高,但現(xiàn)有方法在計算過程中通常只能得到一個方向結(jié)果,無法應(yīng)對具有多個方向特征的建筑。
最小二乘法是以建筑物輪廓節(jié)點或輪廓邊上的點作為觀測值,將建筑物輪廓相鄰邊的垂直關(guān)系作為限制條件,使用總體最小二乘方法求解,獲得直角化后滿足限制條件的節(jié)點坐標(biāo),從而實現(xiàn)建筑物輪廓的直角化[6,19-20],其改進方法有混合最小二乘-總體最小二乘直接解法[21]、加權(quán)總體最小二乘法[22]和限制條件偏最小二乘法[23]。此外,也有研究將最小二乘方法應(yīng)用在點云數(shù)據(jù)提取建筑物輪廓規(guī)則化[24-25]、室內(nèi)地圖邊界校正[26]等方面。相比主方向法,基于最小二乘法的直角化結(jié)果誤差更小,但也存在建模、推導(dǎo)和解算困難的問題[21]。
現(xiàn)有主方向法和最小二乘法解決了具有一個主方向的簡單建筑物直角化問題,然而對于復(fù)雜建筑物及建筑群的規(guī)則化尚沒有較好的應(yīng)對方法。建筑物的規(guī)則特征主要表現(xiàn)為建筑物的方向及相鄰輪廓邊的夾角,當(dāng)建筑物輪廓邊的夾角全部為直角,則輪廓邊的方向可由兩個垂直的方向表示。然而,現(xiàn)實中許多建筑物與其相鄰輪廓邊夾角中存在非直角的情況,需要用多個方向來表達。目前,相關(guān)的研究尚缺乏對這種復(fù)雜建筑輪廓規(guī)則化的解決方法。此外,現(xiàn)有方法一般針對單個建筑物進行規(guī)則化,未考慮相鄰建筑之間的群體關(guān)系,導(dǎo)致處理建筑群時整體規(guī)則化效果不佳。
針對上述問題,本文提出一種基于向量重組的建筑物輪廓方向計算和規(guī)則化方法。首先將建筑物的輪廓邊轉(zhuǎn)化為一組按照一定順序組成的向量集合,根據(jù)向量的方向特征進行分組和變換,對向量組計算方向向量,突破現(xiàn)有方法只能描述建筑物單一主方向的限制。在此基礎(chǔ)上,使用直角改正加強方向的直角特征,進一步修正建筑規(guī)則化的幾何特征。最后,以方向向量為約束條件調(diào)整建筑物輪廓邊,實現(xiàn)建筑物輪廓的規(guī)則化。
基于向量重組的建筑物輪廓規(guī)則化算法流程如圖1所示,主要分為2個部分:方向向量計算和規(guī)則化。首先輸入建筑物原始輪廓,將其轉(zhuǎn)換為一組按照一定順序組成的向量集合,并計算得到初始向量組。然后,根據(jù)向量的方向特征進行分組和變換,將每組向量描述為一個方向特征,進一步通過直角改正加強方向的直角特征,從而得到描述建筑物輪廓方向特征的方向向量。最后,以方向向量為約束條件,計算輪廓邊旋轉(zhuǎn)基點,對建筑物輪廓邊進行調(diào)整,得到建筑物規(guī)則化結(jié)果。
圖1 建筑物輪廓方向計算和規(guī)則化方法流程Fig.1 Flowchart of building footprints orientation calculation and regularization
基于向量重組的建筑物輪廓方向計算方法主要思想是將建筑物輪廓邊視為向量,并設(shè)計一定規(guī)則對每個建筑物輪廓邊進行重組。然后,依據(jù)方向特征對重組邊進行分組,將同一組向量邊變換為同一方向。進而計算每個分組方向,判斷分組方向是否滿足垂直關(guān)系,如果滿足則進行直角改正,從而獲得描述建筑物輪廓方向的一組方向向量。在向量重組和分組的過程中,相對平行的邊分在同一組,計算得到的方向向量表示了同一方向的所有平行邊,因此本文提出的向量重組方法得到的方向向量可以從整體上描述建筑物輪廓的方向特征。
1.2.1 向量重組
向量重組是將建筑物的輪廓邊轉(zhuǎn)化為向量,則建筑物輪廓可以看作為按照一定順序組成的向量集合[27],如圖2(a)所示,建筑多邊形由P1,P2,…,P8共8個點組成,將多邊形的邊表示為向量,則多邊形由向量P1P2,P2P3,…,P8P1組成的集合。然后,將所有向量的起點平移至同一點,因原始建筑物輪廓的規(guī)則特征,重組后的向量具有一定的方向特征。如圖2(b)所示,將向量的起點平移到同一個點O,組成新的向量,如向量P1P2重組為向量OP2,向量P2P3重組為向量OP3。
圖2 向量重組Fig.2 Vectors reconstruction
1.2.2 向量分組和變換
通過向量重組后,因原始建筑物的規(guī)則特征,重組邊具有明顯的方向特征??紤]到原始數(shù)據(jù)或提取誤差的存在,方向接近向量邊的方向值不完全相等。因此需設(shè)定條件閾值,將方位角小于條件閾值或旋轉(zhuǎn)角度π后小于條件閾值的邊分為同一組。根據(jù)建筑物的復(fù)雜程度,可得到多組向量,每一組向量表示建筑物輪廓的一個方向。如圖3所示,重組后向量方位角取值為[0,2π),根據(jù)向量方位角進行分組,分組規(guī)則如下:設(shè)定一個角度閾值θ,向量的方位角α1、α2,若|α1-α2|<θ、|α1-α2-π|<θ或|α2-α1-π|<θ,則將向量分為一組。θ的取值和原始數(shù)據(jù)精度有較大的相關(guān)性,誤差越大,則同一方向邊的角度差值越大,相應(yīng)θ的取值應(yīng)大一些,向量重組后,同一方向的向量有一定的差值,角度閾值θ可依據(jù)該差值確定。根據(jù)分組規(guī)則及閾值設(shè)定,圖3中的原始向量可分為3組:{OP2,OP5,OP8}、{OP1,OP4,OP6}、{OP3,OP7}。對分組后向量的進行變換,使同一組的向量方向相同,如第一組中,OP5、OP8的方向和OP2的方向相反,對OP5、OP8旋轉(zhuǎn)角度π,分別得到OP5′、OP8′。對其余向量用同樣的方法變換,最終向量分為3組:G1={OP2,OP5′,OP8′},G2={OP1,OP4′,OP6},G3={OP3′,OP7}。
圖3 向量分組和變換Fig.3 Vector grouping and transformation
1.2.3 分組方向計算
利用向量重組方法獲得每一組的向量邊后,需要計算每一組的方向值代表該組的整體方向。合理的分組方向值要考慮每個向量邊的長度和方向,向量邊的長度不同,對分組方向的貢獻應(yīng)該也不同[14-15]。本文采用總體最小二乘方法獲得每個分組的方向值:計算每個向量邊終點到分組方向的距離,當(dāng)距離的平方和最小時,則認為該方向為分組方向。
如圖4所示,向量α1,α2,…,αn為一組向量,起點均為點O,假設(shè)其分組方向經(jīng)過點p0,則其向量可表示為p0,向量v為分組方向p0的單位法向量。當(dāng)分組方向滿足總體最小二乘時,向量α1-p0,α2-p0,…,αn-p0在向量v投影長度的平方最小,即向量α1-p0,α2-p0,…,αn-p0和向量v點乘的平方和最小,可得目標(biāo)函數(shù)
圖4 分組方向計算Fig.4 Orientation calculation of group
(1)
表示為矩陣形式
(2)
對列向量p0求導(dǎo)可得
(3)
(4)
1.2.4 直角改正
原建筑物輪廓邊本應(yīng)垂直的方向向量夾角不一定為π/2,在獲得每組方向向量后,仍需對分組方向進一步處理,使應(yīng)垂直的方向向量嚴(yán)格垂直,本文稱該過程為分組方向的直角改正。若兩分組方向向量的方位角|α1-α2-π/2|<θ,則對其進行直角改正,對方位角較大的向量,順時針旋轉(zhuǎn)π/2,對旋轉(zhuǎn)后的向量和方位角較小的向量求和,求和后向量的方向為直角改正后方位角較小向量的方向,逆時針旋轉(zhuǎn)π/2后的方向為方位角較大的向量方向,分組向量長度取原向量在改正后方向向量的投影長度,最終得到直角改正后的分組向量。由式(1)—式(4)的推導(dǎo)可知,直角改正過程滿足總體最小二乘。
對圖5(a)所示的建筑物輪廓計算方向,計算分組向量和:OG1=OP2+OP5′+OP8′,OG2=OP1+OP4′+OP6,OG3=OP3′+OP7。分組向量OG1和OG2滿足直角改正的條件,將向量OG2順時針旋轉(zhuǎn)π/2后得到向量OR2,對向量OG1和向量OR2求和,取求和后向量的方向和OG1在求和向量上的投影長度,得到OG1直角改正后的向量OGR1,將OGR1逆時針旋轉(zhuǎn)π/2,并取OG2的在該旋轉(zhuǎn)向量的投影長度,得到OG2直角改正后的向量OGR2,得到建筑物輪廓的最終方向向量:OGR1、OGR2和OG3。
圖5 方向計算和直角改正Fig.5 Orientation calculation and right angle correction
在計算得到建筑物方向向量的基礎(chǔ)上,構(gòu)建建筑物規(guī)則化通用公式,基于方向向量的旋轉(zhuǎn)基點進行節(jié)點坐標(biāo)的調(diào)整,從而實現(xiàn)建筑物規(guī)則化。如圖6所示,l為待調(diào)整邊,l1和l2為其相鄰邊,l′為其通過向量重組計算出的方向,對建筑物輪廓規(guī)則化,就是求解一點P,將待調(diào)整邊l繞點P旋轉(zhuǎn)到和l′相同的方向,P為旋轉(zhuǎn)基點。
圖6 輪廓邊旋轉(zhuǎn)基點Fig.6 The rotation base point of footprints edge
S=(x1-x)2+[y1-(k1x+b1)]2+
(x2-x)2+[y2-(k2x+b2)]2
(5)
將y1=k1x1+b1和y2=k2x2+b2代入式(5),可得
(6)
對x求導(dǎo)可得
(7)
令導(dǎo)數(shù)等于0,可求得
(8)
設(shè)l′和l1的夾為α1,設(shè)l′和l2的夾為α2,l′的方位角為π/2,k1和α1的關(guān)系為k1=cotα1,k2和α2的關(guān)系為k2=cotα2,代入式(8)可得
(9)
式(9)表示出了求解點P橫坐標(biāo)的公式,以l′的方向為橫坐標(biāo)軸方向,同樣的推導(dǎo)方法可得縱坐標(biāo)公式
(10)
式(10)是假設(shè)l1和l2的斜率都存在情況下得到的,當(dāng)斜率不存在時,夾角α=π/2,此時有cotα=0,代入式(10)即可得到斜率不存在時縱坐標(biāo)的計算公式。
通過式(9)和式(10)即可求解出旋轉(zhuǎn)基點P的坐標(biāo)。式(9)和式(10)表明,當(dāng)P1和P2坐標(biāo)固定,則點P坐標(biāo)只和相鄰兩條邊的夾角有關(guān)。當(dāng)坐標(biāo)系發(fā)生旋轉(zhuǎn)時,如圖7所示,相鄰邊的夾角不會改變,點P相對于P1和P2的位置關(guān)系不會改變,式(9)和式(10)對旋轉(zhuǎn)后的坐標(biāo)系同樣適用,因此式(9)和式(10)是解算旋轉(zhuǎn)基點P的一般公式。
圖7 坐標(biāo)系旋轉(zhuǎn)Fig.7 Coordinate system rotation
特別地,當(dāng)α1=α2時,即待調(diào)整邊和相鄰兩個邊的夾角相同時,有
(11)
(12)
式(11)和式(12)表示,當(dāng)待調(diào)整邊和相鄰兩個邊夾角相同時,待調(diào)整邊的中點即為旋轉(zhuǎn)基點。
圖8 建筑物規(guī)則化方法Fig.8 Building regularization method
為驗證本文方法的有效性和效率,本文首先對復(fù)雜建筑和建筑群進行規(guī)則化,然后和主方向法進行對比分析,最后對算法效率分析。
現(xiàn)有方法難以實現(xiàn)具有多方向特征復(fù)雜建筑物的規(guī)則化,選取形狀較復(fù)雜的建筑物輪廓進行方向計算和規(guī)則化試驗。圖9(a)為原始建筑輪廓經(jīng)過向量重組變換后的重組邊和方向向量對比(為便于展示,方向向量的長度進行了縮放),原始重組邊的方向有一定的差別,經(jīng)過方向計算和直角改正后,最終得到3組相互垂直的方向向量。以方向向量為控制條件,對建筑輪廓進行規(guī)則化處理,結(jié)果如圖9(b)所示,試驗結(jié)果證明方法對于復(fù)雜建筑同樣適用,在加強建筑物規(guī)則特征的同時保持了原始形狀特征。
圖9 復(fù)雜建筑物方向向量和規(guī)則化結(jié)果Fig.9 Complex building orientation vectors and regularization results
鄰近的建筑物方向通常具有相似性,尤其是多個建筑組成的建筑群,其中部分建筑物往往具有相同的方向。如圖10(a)所示為多個建筑組成的建筑群和單個建筑物的方向向量,可將建筑群作為一個整體,得到建筑群的方向向量,結(jié)果如圖10(b)所示,最終得到3組垂直的方向向量,建筑方向向量的角度值見表1。
表1 建筑群方向向量的角度值
圖10 建筑群方向向量Fig.10 Orientation vectors of building group
從建筑方向向量中可得到建筑物群的方向關(guān)系,如建筑1、2、3具有近似的方向特征,然而,由于原始采集誤差的存在,逐一計算單個建筑物方向會導(dǎo)致相同方向建筑物的方向向量存在一定的差別。將建筑群作為整體考慮,計算得到的方向向量更能代表建筑群的方向特征,以方向3為例,有方向3的建筑包括建筑4—建筑10共7個建筑,每個建筑方向3的具體方向值不相等,因此,考慮更多建筑的整體方向值更加準(zhǔn)確。
以得到的整體方向向量為控制條件,對每個建筑輪廓進行規(guī)則化,規(guī)則化結(jié)果如圖11所示。相比于對單個建筑規(guī)則化,以建筑群為整體得到的規(guī)則化結(jié)果,具有接近方向特征的建筑輪廓方向完全相同,因此規(guī)則化結(jié)果更加合理。
圖11 建筑群規(guī)則化結(jié)果Fig.11 Regularization results of building group
為進一步驗證本文方法,將本文方法與主方向法進行對比分析?;谧钚⊥饨泳匦畏ㄓ嬎愕玫降姆较蜉^為穩(wěn)健且應(yīng)用較多[7-10,18],本文使用最小外接矩形獲得建筑物主方向。主方向法只能實現(xiàn)建筑物輪廓的直角化,以建筑物輪廓的直角化來對比兩種方法,在直角化的過程中,應(yīng)盡量減少節(jié)點的移動距離,因此通過節(jié)點移動距離的均方根誤差評價兩種方法。試驗建筑物輪廓有6個節(jié)點且相鄰輪廓邊相互垂直,節(jié)點坐標(biāo)見表2,以該建筑物的節(jié)點坐標(biāo)作為真實值,通過加入高斯誤差的方式獲得模擬數(shù)據(jù)。為充分對比兩種方法,共進行1000次模擬試驗,加入均值為0、標(biāo)準(zhǔn)差為0.5 m的高斯誤差。加入高斯誤差后,相鄰輪廓邊不再滿足垂直條件,對模擬結(jié)果進行直角化處理,角度閾值取值為π/12,計算直角化前后節(jié)點移動距離的均方根誤差。原始模擬結(jié)果節(jié)點均方根誤差分布如圖12(a)所示,為直觀對比兩種方法誤差的分布規(guī)律,以主方向法均方根誤差排序結(jié)果如圖12(b)所示,以本文方法均方根誤差排序結(jié)果如圖12(c)所示。
表2 試驗建筑物節(jié)點坐標(biāo)
圖12 主方向法和本文方法1000次模擬節(jié)點均方根誤差Fig.12 Coordinate RMSE of main orientation method and the proposed method for 1000 simulations
在1000次模擬的直角化結(jié)果中,主方向法節(jié)點移動最大均方根誤差為0.972 3 m,平均值為0.324 9 m,本文方法最大均方根誤差為0.736 8 m,平均值為0.305 1 m。由圖12(b)可知,以主方向法均方根誤差從小到大排序,本文方法的誤差曲線主要在主方向法誤差曲線的下方波動;由圖12(c)可知,以本文方法均方根誤差從小到大排序,主方向的誤差曲線在本文方法誤差曲線的上方波動。
進一步和真實值比較,計算原始模擬數(shù)據(jù)、主方向法和本文方法直角化結(jié)果節(jié)點坐標(biāo)相對于真實值節(jié)點坐標(biāo)的均方根誤差,原始模擬結(jié)果節(jié)點均方根誤差分布如圖13(a)所示,為展示原始數(shù)據(jù)和直角化前后誤差的變化分布規(guī)律,對原始模擬數(shù)據(jù)誤差排序結(jié)果如圖13(b)所示。
圖13 相對于建筑物節(jié)點真實值均方根誤差Fig.13 RMSE relative to real building coordinates
在本文1000次模擬試驗中,相對于真實數(shù)據(jù),原始模擬數(shù)據(jù)節(jié)點均方根誤差最大值為1.028 2 m,平均值為0.481 2 m;主方向法直角化后節(jié)點均方根誤差最大值為0.892 1 m,平均值為0.385 1 m;本文方法直角化后節(jié)點均方根誤差最大值為0.822 1 m,平均值為0.366 9 m。從1000次模擬結(jié)果節(jié)點的均方根誤差平均值來看,主方向法將原始模擬數(shù)據(jù)的均方根誤差減小了19.97%,本文方法將原始模擬數(shù)據(jù)的均方根誤差減小了23.75%。進一步對誤差分布圖分析可以看出,主方向法直角化后大部分建筑的均方根誤差得到了減小,但出現(xiàn)部分結(jié)果均方根誤差增大的情況,在1000次模擬試驗中,共出現(xiàn)80次,增大的最大值為0.159 6 m。本文方法直角化的均方根誤差分布曲線主要在原始模擬數(shù)據(jù)均方根誤差曲線的下方波動,只有極個別出現(xiàn)直角化后均方根誤差增大的情況,在1000次模擬試驗中,共出現(xiàn)2次,增大的最大值為0.001 2 m。因此,本文方法通過直角化方法對建筑輪廓直角化不僅可以從視覺上得到更好的結(jié)果,還可以進一步提高原始數(shù)據(jù)的整體精度。在精度提高和整體穩(wěn)定性方面,本文方法優(yōu)于主方向法。
本文方法的時間復(fù)雜度主要與建筑物輪廓的方向個數(shù)有關(guān),具體每個步驟的時間復(fù)雜度及說明見表3。
表3 算法的時間復(fù)雜度
試驗對約16 km2的重慶市中心城區(qū)建筑物進行規(guī)則化,建筑物輪廓如圖14所示。試驗區(qū)建筑物輪廓數(shù)據(jù)取自1∶2000比例尺地形圖,采用立體航測法采集,立體航片分辨率為0.2 m,建筑物采集精度為0.2 m,其中包含建筑物11 214個,共有139 658條邊,平均每個建筑物12.5條邊。測試計算機配置為Intel(R) Core(TM) i5-10400F CPU @2.90 GHz,內(nèi)存8 GB,角度閾值取值為π/20。
圖14 試驗區(qū)建筑物輪廓Fig.14 Building footprints of experiment area
對每個建筑物計算方向向量和規(guī)則化,處理總時間為2 min 38 s,平均每個建筑物耗時為0.014 s,算法效率較高。建筑物方向向量個數(shù)統(tǒng)計見表4,方向向量總個數(shù)為31 290個,平均每個建筑物2.8個方向向量,方向向量為2個的建筑物共8111個,占建筑物總數(shù)的72.33%,方向向量總個數(shù)為4個及以下的占建筑物總數(shù)的88.67%。
表4 試驗數(shù)據(jù)方向向量個數(shù)統(tǒng)計
本文針對復(fù)雜建筑物和建筑群輪廓的方向定量描述和規(guī)則化等問題,從建筑物輪廓的幾何規(guī)則特征分析出發(fā),提出了一種基于向量重組的建筑物輪廓方向計算和規(guī)則化方法。在方向計算方面,通過向量重組獲得描述建筑物輪廓方向特征的重組邊,經(jīng)過向量變換和方向計算,實現(xiàn)了復(fù)雜建筑物方向的定量描述。在此基礎(chǔ)上,以方向向量為控制條件,通過計算旋轉(zhuǎn)基點,進一步調(diào)整節(jié)點坐標(biāo),實現(xiàn)了復(fù)雜建筑物和建筑群的規(guī)則化。
通過推導(dǎo)驗證,本文方法在方向計算和規(guī)則化過程中均滿足總體最小二乘,有效降低了直接構(gòu)建總體最小二乘模型的復(fù)雜解算困難。通過與主方向法對比試驗可知,本文提出的向量重組方法在節(jié)點整體移動距離和穩(wěn)定性方面更具有優(yōu)勢,通過規(guī)則化算法可進一步提高數(shù)據(jù)的整體精度。通過復(fù)雜建筑物和建筑群的試驗,進一步驗證了本文方法在復(fù)雜建筑物和建筑群規(guī)則化的適用性,為建筑群空間分布特征識別、智慧城市構(gòu)建、城市建筑變化檢測等應(yīng)用提供一種解決思路。