王玉堂,王佳,2*,牛利偉,常書萍,孫露
(1. 北京林業(yè)大學(xué)林學(xué)院,北京 100083; 2. 北京林業(yè)大學(xué)精準(zhǔn)林業(yè)北京市重點實驗室,北京 100083)
樹冠體積和表面積的傳統(tǒng)計算方法主要是通過伐倒木解析或使用測繪工具測量樹木的冠幅、冠高等參數(shù),根據(jù)樹冠形態(tài)將其擬合成近似的幾何圖形,根據(jù)相應(yīng)的體積公式進行求解[1-2]。傳統(tǒng)的計算方法費時費力,破壞性強,并且由于樹冠結(jié)構(gòu)的不規(guī)則性,傳統(tǒng)方法計算精度不高,很難滿足實際需求[3]。三維激光掃描技術(shù)具有高精度、穿透性、不接觸性等特點,為了提高測量精度,許多學(xué)者將其引入樹冠的測量中[4-5]。
近年來,無人機航空攝影測量技術(shù)興起,廣泛應(yīng)用于地籍測繪、城市三維建模、土方量計算等多個領(lǐng)域。由于無人機遙感具有機動、靈活、低成本、安全性高等優(yōu)勢,不斷發(fā)展成為林業(yè)調(diào)查中的主要手段[6-7]。其中正射影像結(jié)合樹木三維點云可以精確地提取樹木高度信息[8],傾斜攝影測量能夠快速、高效地獲取樹木不同角度的特征信息,真實地反映樹木的三維結(jié)構(gòu)[9]?;谝曈X的三維重建技術(shù)不受物體形狀及場景的限制,可快速實現(xiàn)全自動或半自動建模。其中,基于多視點的運動恢復(fù)結(jié)構(gòu)(structure from motion, SFM)算法,利用不同角度并且具有一定重疊度的一系列影像數(shù)據(jù),通過特征點匹配來求解相機姿態(tài)參數(shù)和變換矩陣,從而實現(xiàn)物體三維幾何信息的恢復(fù)[10]。正是由于SFM算法的通用性,其不依賴于某一特定場景,所以無人機傾斜攝影測量的多角度特征與SFM算法相結(jié)合為林業(yè)調(diào)查提供了一種新思路。
曾健等[11]證明了使用無人機傾斜攝影測量和SFM算法生成的點云數(shù)據(jù)來提取落葉松人工林地形信息的可行性;Wallace等[12]對比分析了使用機載激光雷達(ALS)和SFM算法兩種方法生成的樹木點云數(shù)據(jù),兩種方法都可以提取樹木特征信息,但ALS方法生成的樹木點云點密度更高,林下信息更豐富,提取精度也更高。但由于在實際應(yīng)用中三維激光掃描儀成本高操作復(fù)雜等缺點,并且在低郁閉度的林區(qū)兩者林分信息的提取精度較為接近[13],所以使用無人機傾斜攝影測量的多角度影像數(shù)據(jù)結(jié)合SFM算法進行的三維重建,是林分調(diào)查時的最佳低成本替代方法。
基于高精度點云數(shù)據(jù)提取樹冠體積和表面積的算法有很多,應(yīng)用較為廣泛的主要有體元法[14]和數(shù)字高程模型法[15]。這兩種算法都有其優(yōu)缺點,體元法充分考慮了樹冠內(nèi)部空隙,但該方法無法區(qū)分樹冠真實空隙和因遮擋而形成的偽空隙,而且在樹冠表面規(guī)則性較差時不能很好地體現(xiàn)樹冠的邊緣特征,并且操作復(fù)雜計算過程較慢;數(shù)字高程模型法利用了點云數(shù)據(jù)的高程信息,并且結(jié)合數(shù)字化測繪數(shù)據(jù)處理系統(tǒng)生成立體三角網(wǎng),操作方便快捷,還可以很好地反映樹冠的邊緣特征,但其進行計算時忽略了樹冠內(nèi)部的空隙,而且三角網(wǎng)容易出現(xiàn)過擬合現(xiàn)象。
本研究使用多旋翼航測無人機對北京市懷柔區(qū)主要道路的行道樹進行傾斜攝影測量和樹木模型的三維重建,在獲取單木點云數(shù)據(jù)的基礎(chǔ)上,對比分析兩種樹冠體積和表面積的提取算法在不同冠型、不同點云密度等方面的差異,以期為無人機傾斜攝影測量提取樹冠體積和表面積的算法研究提供參考。
試驗地位于北京市懷柔區(qū)(116°17′~116°63′E,40°41′~41°04′N)中心主干道,屬于典型的北溫帶半濕潤大陸性季風(fēng)氣候,夏季高溫多雨,冬季寒冷干燥,春、秋短促,年均氣溫9~13 ℃,年均降雨量600~700 mm。試驗區(qū)的行道樹主要樹種為洋白蠟(Fraxinuspennsylvanica)、側(cè)柏(Platycladusorientalis)、銀杏(GinkgobilobaL.)、美桐(PlatanusoccidentalisL.)、刺槐(RobiniapseudoacaciaLinn.)、毛白楊(Populustomentosa)等代表性樹種。
1.2.1 試驗地選取
無人機飛行作業(yè)需遠(yuǎn)離強干擾物及高大建筑物等干擾物,并且由于研究精力有限,所選行道樹要具有代表性,樹木間遮擋較少,以便單木點云數(shù)據(jù)的準(zhǔn)確提取。由于市區(qū)環(huán)境復(fù)雜,經(jīng)過飛行前的實地考察,綜合考慮以上因素后,在符合條件的區(qū)域選擇了5段50 m長的道路作為實驗地,路旁的行道樹作為實驗樹種。
1.2.2 無人機數(shù)據(jù)獲取
無人機飛行時間為2019年10月19日和20日,這兩天天氣晴朗風(fēng)力較小。為避免光照條件對實驗結(jié)果的影響,每天的飛行時間選為10:00和13:00 兩個光照較好的時間段。為保證測量結(jié)果的準(zhǔn)確性,防止傾斜攝影測量精度的不同對研究結(jié)果產(chǎn)生影響,本研究在實際測量每塊樣地時控制了變量的統(tǒng)一,保證了每次傾斜攝影測量的高度、傾斜角度以及飛行速度等參數(shù)不變。具體參數(shù)為:飛行高度30 m,航向重疊度80%,旁向重疊度70%,飛行速度為1 m/s。數(shù)據(jù)獲取平臺為大疆精靈4 RTK(公司:大疆創(chuàng)新;產(chǎn)地:中國深圳)行業(yè)級測繪無人機,無人機飛行航線示意圖見圖1(以樣地4為例)。
圖1 無人機飛行航線Fig. 1 Schematic diagram of UAV flight routes
圖3 單木點云圖Fig. 3 Single tree point cloud
利用攝影測量軟件將航拍影像生成點云模型,并對點云數(shù)據(jù)進行去噪、單木分割、樹冠提取等數(shù)據(jù)預(yù)處理操作,對處理好的點云數(shù)據(jù)分別使用體元法和數(shù)字高程模型法計算樹冠體積和表面積,實驗流程見圖2。
圖2 實驗流程Fig. 2 Experimental flow chart
1.3.1 數(shù)據(jù)預(yù)處理
1)選擇Pix4D-mapper(版本號為2.0.104;公司為Pix4D;產(chǎn)地為瑞士洛桑)作為航拍影像處理軟件,導(dǎo)入影像數(shù)據(jù),軟件自動讀取照片位置信息。Taddia等[16]在實地實驗中測試了DJI Phantom 4 RTK傾斜攝影測量數(shù)據(jù)在有像控點和沒有像控點情況下的均方根誤差(RMSE),結(jié)果表明兩者僅相差0.003 m??梢钥闯觯瑹o人機在RTK GNSS下不使用像控點和使用像控點時測量結(jié)果的誤差非常小。同時,為了減小缺少像控點帶來的誤差,將飛行高度設(shè)為30 m。選擇傾斜攝影測量建模方式,設(shè)置坐標(biāo)系和投影,默認(rèn)為WGS 84 / UTM zone 50N投影坐標(biāo)系。軟件自動進行初始化,提取并匹配特征點生成樹木三維點云數(shù)據(jù)。
2)在Pix4D-mapper軟件中的點云模式下,使用點云編輯工具人工手動進行單木分割(圖3)和樹冠部分點云的選取,并將分割好的樹冠點云數(shù)據(jù)以.las格式進行輸出。由于圖片分辨率和航線重疊度較高,傾斜攝影測量生成的點云具有較高的質(zhì)量,本研究生成的樹木點云平均密度為5 043.8點/m3。
3)使用編程語言Python選擇基于統(tǒng)計的方法,根據(jù)相鄰點之間的距離來去除噪聲點。
1.3.2 體元法計算樹冠體積和表面積
體元法計算樹冠體積和表面積的基本思想是將樹冠劃分成無數(shù)個小立方體,選擇適合大小的體元就能很好地模擬出樹冠形態(tài)及內(nèi)部結(jié)構(gòu)。先給樹冠分層,以樹冠的最高點和最低點連線作為Z軸,以k為間隔將樹冠劃分為n層,將每一層的點云投影到垂直于Z軸的XY平面,同樣沿X軸和Y軸以k為間隔劃分體元。根據(jù)韋雪花等[14]的實驗可知,當(dāng)體元邊長k≤冠徑/10時計算結(jié)果較好,計算用C#語言編程實現(xiàn),體元法計算示意圖見圖4。
圖4 體元法原理Fig. 4 The principle of voxel method
樹冠體積計算時,在每一橫截面上逐個體元判斷是否含有點云,含有點云的體元為有效體元,記為1,相反,不含點云的像元記為0,統(tǒng)計每一層有效像元個數(shù)T,則樹冠體積計算公式為:
(1)
式中:V為樹冠體積,m3;k為體元邊長,m;n為樹冠層數(shù);Ti為第i層中的有效體元個數(shù)。
樹冠表面積計算時,首先建立每層中各個體元間的拓?fù)潢P(guān)系,從頭到尾遍歷各個體元記錄當(dāng)前體元與其他體元的鄰接個數(shù),如果與4個體元鄰接就說明有0條獨立邊,與3個體元鄰接就說明有1條獨立邊,與2個體元相鄰接則有2條獨立邊,與1個體元鄰接則有3條獨立邊,最后統(tǒng)計所有體元的鄰接個數(shù)得出每一層的獨立邊個數(shù)。設(shè)第i層獨立邊的個數(shù)為Ji, 統(tǒng)計所有層的獨立邊的個數(shù),則樹冠表面積計算公式如下:
(2)
式中:S為樹冠表面積,m2;Ji為每一層的獨立邊個數(shù)。
1.3.3 數(shù)字高程模型法計算樹冠體積和表面積
數(shù)字高程模型法計算樹冠體積和表面積充分利用了點云數(shù)據(jù)所具有的高程信息,計算原理為:根據(jù)樹冠點云數(shù)據(jù)的高程值,將具有相同高度的樹冠表面的點用平滑的曲線連接,形成等高線,再根據(jù)等高線繪制方格網(wǎng),其中角點處的體積為方格體積的1/4,邊點體積為方格體積的2/4,拐點體積為方格體積的3/4,中心點體積為方格體積,所有點的體積之和為樹冠體積。計算表面積時應(yīng)選取樹冠表面點,利用其高程值連接生成不規(guī)則三角網(wǎng),分別計算每個三角形的面積,所有三角形面積之和即為樹冠表面積。
數(shù)字高程建模工具選擇的是基于CAD成圖軟件二次開發(fā)的CASS軟件(版本號:9.0;公司:廣東南方數(shù)碼科技股份有限公司;產(chǎn)地:中國廣州),具體計算流程為:①樹冠體積計算。首先將樹冠點云數(shù)據(jù)轉(zhuǎn)換為數(shù)字化測繪數(shù)據(jù)采集系統(tǒng)可以識別的.dxf格式,打開點云數(shù)據(jù)并提取高程保存為dat格式的高程數(shù)據(jù)。在軟件中展開高程數(shù)據(jù)建立數(shù)字高程模型如圖4所示,立體三角網(wǎng),用閉合復(fù)合線圈出要計算樹冠體積的區(qū)域,利用CASS軟件的DTM法土方計算功能計算出點云體積。②樹冠表面積計算。首先利用體元法選擇具有1,2,3條獨立邊的體元,將具有獨立邊的體元內(nèi)的點云作為樹冠表面點云。再根據(jù)體積計算時高程數(shù)據(jù)提取的方法提取高程數(shù)據(jù),建立三角網(wǎng),選取表面積計算區(qū)域,利用CASS軟件表面積計算功能,根據(jù)圖上高程點來計算樹冠表面積。
1.3.4 幾何法計算樹冠體積和表面積
由于樹冠體積和表面積的真實值無法準(zhǔn)確測量,所以在傳統(tǒng)的樹冠體積和表面積的計算中,根據(jù)樹冠的冠型將其模擬成規(guī)則的幾何體,以樹冠的冠幅和冠高為參數(shù),計算幾何體的體積和表面積作為樹冠的體積和表面積。本研究中根據(jù)實驗樹種的冠型將其分別近似為橢圓形、半球形和圓錐形,每種冠型的體積和表面積計算公式見表1。
表1 幾何法樹冠體積和表面積計算公式Table1 Geometrical formulas for calculating crown volumes and surface areas
對傾斜攝影測量得到的點云數(shù)據(jù)預(yù)處理后得到了50棵單木樹冠的點云數(shù)據(jù),共6個樹種。分別使用體元法、數(shù)字高程模型法和傳統(tǒng)幾何法計算樹冠體積和表面積。由于樹冠形態(tài)的不規(guī)則性以及實際操作的困難性,在不損傷樹木的情況下樹冠體積和表面積的真實值是無法測量的,所以不對算法的計算精度進行比較。對比了3種方法的計算結(jié)果,從不同樹種間的計算結(jié)果來看,美桐的3種方法計算結(jié)果相差最大,側(cè)柏最小(圖5)。從不同樹種的冠型分析,側(cè)柏的冠型更接近模擬的圓錐形,所以傳統(tǒng)方法與體元法和數(shù)字高程法的計算結(jié)果間的誤差較小。而美桐的冠型規(guī)則性較差,通過與規(guī)則幾何體對比,傳統(tǒng)幾何法只能將其近似模擬為圓錐形,與實際冠型間仍有較大差距,所以計算結(jié)果誤差較大。從樹冠體積計算結(jié)果來看,傳統(tǒng)幾何法最大,數(shù)字高程法次之,體元法最小。從3種方法的原理分析,傳統(tǒng)幾何法沒有考慮樹冠表面不規(guī)則性和樹冠內(nèi)部的空隙,所以計算結(jié)果較大;而數(shù)字高程法僅考慮了樹冠表面的不規(guī)則性,但沒有考慮樹冠內(nèi)部的空隙,所以其計算結(jié)果小于傳統(tǒng)幾何法大于體元法。樹冠表面積的計算結(jié)果中,體元法和傳統(tǒng)幾何法較為接近,這是因為這兩種方法將樹冠表面模擬得較為規(guī)則,而數(shù)字高程法則根據(jù)樹冠表面真實形狀進行計算。
圖5 樹冠體積和表面積計算結(jié)果Fig. 5 Comparisons of tree crown volume and surface area calculation results
通過分析可知傳統(tǒng)幾何法的計算結(jié)果誤差較大,所以本研究只分析體元法和數(shù)字高程法兩種算法的差異及影響因素。使用兩種算法計算結(jié)果間的相對誤差和絕對誤差來定量表示它們間的差異,具體計算公式為:
(3)
Eb=|r1-r2|
(4)
式中:Er為相對誤差;Eb為絕對誤差;r1為數(shù)字高程法計算結(jié)果;r2為體元法計算結(jié)果。
樹冠體積和表面積計算結(jié)果見表2。由表2可知,除側(cè)柏的樹冠體積外,體元法計算的樹冠體積和表面積均小于數(shù)字高程法的計算結(jié)果,兩種算法計算的樹冠體積相對誤差最大為銀杏的42.87%,最小為側(cè)柏的5.60%;絕對誤差最大為美桐的40.76 m3,最小為側(cè)柏的1.23 m3。樹冠表面積相對誤差最大為洋白蠟的66.22%,最小為側(cè)柏的33.47%;絕對誤差最大為美桐的109.98 m2,最小為側(cè)柏的29.38 m2。表2計算結(jié)果表明,銀杏等樹種的樹冠內(nèi)部空隙較多,側(cè)柏的內(nèi)部空隙較少,所以兩種算法在樹冠體積的計算上銀杏的相對誤差最大,側(cè)柏的相對誤差最小,而由于樹冠表面不規(guī)則性的存在,導(dǎo)致兩種算法間的相對誤差均在30%以上。
表2 樹冠體積和表面積計算結(jié)果Table 2 Calculation results of canopy volumes and surface areas
樹冠的內(nèi)部結(jié)構(gòu)主要體現(xiàn)在樹冠的點云密度,點云密度的大小主要由單位體積的點云數(shù)量來體現(xiàn),為便于分析對其進行均值標(biāo)準(zhǔn)化,計算方法如下:
(5)
(6)
式中:i為樹冠編號,取值為1,2,…,n;n為樹冠個數(shù),本研究取50;Pi為點云密度;Pm為標(biāo)準(zhǔn)化后的點云密度;N為點云點數(shù);V1為數(shù)字高程法計算的體積,m3;V2為體元法計算的體積,m3。
計算各樹冠的點云密度,分別對點云密度與兩種算法計算的樹冠體積和表面積的相對差值進行擬合(圖6)。從圖6可以看出,點云密度與兩種樹冠體積算法的差值之間具有較好的相關(guān)性,決定系數(shù)為0.443 1,呈指數(shù)負(fù)相關(guān)關(guān)系,即點云密度越大,兩種樹冠體積算法的計算結(jié)果越接近。由兩種樹冠體積的計算原理可知,體元法計算樹冠體積時充分考慮了樹冠的內(nèi)部空隙,而數(shù)字高程模型法在計算時只考慮了樹冠的外部形態(tài),忽略了樹冠的內(nèi)部結(jié)構(gòu),所以當(dāng)樹冠的點云密度越大時,樹冠內(nèi)部的空隙越小,樹冠越接近一個實體,兩種算法的計算結(jié)果差別也就越小。然而點云密度與兩種表面積算法之間的決定系數(shù)為0.039 2,兩者之間基本沒有相關(guān)性,這主要是由于兩種表面積算法在計算時都只考慮了樹冠的表面點云,與點云密度不相關(guān),只與樹冠的外形及樹冠表面特征點的選取有關(guān)。
圖6 點云密度對不同算法間相對誤差的影響Fig. 6 Influence of branch and leaf density on relative errors between different algorithms
同一樹種的樹冠形態(tài)大致相同,通過觀察將實驗地的6種樹木依據(jù)其外形劃分為6種不同冠型,不同樹種間相對誤差的統(tǒng)計結(jié)果見表3。由統(tǒng)計結(jié)果可知,兩種算法在計算不同樹種間結(jié)果的相對誤差不同,并且體積計算時不同樹種間的差異大于表面積。在樹冠體積的計算中,不同樹種間相對誤差平均值最小的是窄圓錐形的側(cè)柏樹冠,為14.95%;最大的是橢圓形的銀杏樹冠,為56.55%。在樹冠表面積的計算中,不同樹種間相對誤差平均值最小的也是窄圓錐形的側(cè)柏樹冠,為34.92%;最大是橢圓形的洋白蠟樹冠,為55.69%。
兩種樹冠體積算法在同一樹種間相對誤差最大的為橢圓形的銀杏樹冠,差異最小的為寬圓錐形的美桐樹冠;兩種樹冠表面積算法在同一樹種間相對誤差差異較大的為圓錐形的毛白楊樹冠,差異較小的為寬圓錐形的美桐樹冠。由于樹冠的冠型并不是規(guī)則的幾何體,在不同高度處的冠幅長短不一,為研究樹冠表面規(guī)則性對兩種樹冠體積和表面積算法的影響,本研究以0.5 m為間隔在垂直方向上將樹冠分層,計算每一層平均冠幅之間的標(biāo)準(zhǔn)差作為樹冠的形狀參數(shù)(K),定量表示樹冠表面的規(guī)則性,具體計算方法如下:
表3 不同樹種計算結(jié)果相對誤差Table 3 Statistical table of relative errors of calculation results of different tree species
(7)
式中:K為樹冠的形狀參數(shù),m;n為樹冠的層數(shù);Ci為i層樹冠的平均冠幅,m;Ca為樹冠各層平均冠幅的均值,m。
分析樹冠表面的規(guī)則性與兩種體積和表面積計算結(jié)果間的關(guān)系,分別將冠型參數(shù)與兩種樹冠體積和表面積算法計算結(jié)果間的相對差值進行擬合(圖7)。從圖7可以看出,樹冠表面的規(guī)則性對樹冠體積和表面積的計算都有影響,兩者都呈正相關(guān),且與體積相對差值間的決定系數(shù)為0.349 6,與表面積相對差值間的決定系數(shù)為0.100 1,樹冠表面的規(guī)則性對樹冠體積計算的影響大于對樹冠表面積計算的影響。在樹冠體積和表面積計算時,體元法將樹冠劃分為一個個的體元,雖能反映樹冠表面結(jié)構(gòu),但樹冠表面越不規(guī)則其計算精度越差,而數(shù)字高程法能很好地反映樹冠表面的規(guī)則性,計算結(jié)果受樹冠表面影響較小。
圖7 樹冠形狀參數(shù)對不同算法間相對誤差的影響Fig. 7 The influence of crown shape parameters on the relative errors between different algorithms
本研究利用傾斜攝影測量得到的樹冠點云數(shù)據(jù)對兩種計算樹冠體積和表面積的算法進行了對比分析,根據(jù)計算結(jié)果及兩種算法的原理分析了兩種算法在計算結(jié)果上產(chǎn)生差異的原因,結(jié)果表明:
1)數(shù)字高程法計算的樹冠體積和表面積明顯大于體元法的計算結(jié)果。樹冠內(nèi)部空隙是影響樹冠體積計算結(jié)果的重要因素,洋白蠟和側(cè)柏等樹冠內(nèi)部空隙小的樹種計算結(jié)果較為接近,內(nèi)部空隙大的銀杏樹種的計算相對誤差高達42.87%,而樹冠表面的不規(guī)則性導(dǎo)致兩種算法在樹冠表面積的計算中相對誤差均在30%以上。
2)點云密度對兩種樹冠表面積算法的結(jié)果影響較小,對樹冠體積算法的結(jié)果影響較大。樹冠點云密度與兩種樹冠體積計算結(jié)果相對差值間的相關(guān)系數(shù)為0.443 1,點云密度越大,樹冠內(nèi)部空隙越小,數(shù)字高程法計算的結(jié)果越接近體元法,兩者間相對差值越小。
3)兩種算法在不同樹種間計算結(jié)果差異較大。在計算側(cè)柏樹冠時,兩種體積和表面積算法間的差異均小于其他樹種,而在樹冠體積計算時銀杏樹冠的差異最大,在樹冠表面積計算時洋白蠟最大。
4)樹冠表面的規(guī)則性與兩種體積和表面積計算結(jié)果均具有相關(guān)性。冠型參數(shù)與體積相對差值間的決定系數(shù)為0.349 6,但與表面積相對差值間的相關(guān)性較差,決定系數(shù)為0.100 1。冠型參數(shù)越大,樹冠表面越不規(guī)則,體元法不能很好地反映冠型變化,兩種算法間的差值變大。
本研究分別從樹種冠型、樹冠表面規(guī)則性和點云密度等方面對兩種算法進行了對比分析,為實際樹冠體積和表面積的計算提供理論支持。雖然從兩種算法的理論上來說,體元法能夠更真實反映精度,但是計算過程較慢,另外需要自行二次開發(fā)計算工具,相對麻煩,而利用數(shù)字高程法,可以依賴現(xiàn)有的軟件提供的體積表面積計算工具,快速得到樹冠結(jié)構(gòu)因子。所以,對于樹木內(nèi)部空間較小、冠型較規(guī)則的樹木(例如側(cè)柏),或?qū)τ诓豢紤]樹木內(nèi)部結(jié)構(gòu)的研究(例如三維綠量等),數(shù)字高程法仍有較好的應(yīng)用前景。如果要求計算結(jié)果準(zhǔn)確,或樹冠內(nèi)部空隙較多,則體元法較為適合。
另外本研究只選擇了北京地區(qū)6個常見樹種夏季的測量數(shù)據(jù),并沒有本地區(qū)所有樹種以及同一樹種其他季節(jié)的數(shù)據(jù),在算法中只選擇了兩種算法來對比分析。在今后的研究中將進一步增加實驗數(shù)據(jù),并且進一步提高算法的速度和精度。