張衛(wèi)正 張偉偉 李燦林 萬瀚文 張秋聞 劉 巖 金保華
(1.鄭州輕工業(yè)大學計算機與通信工程學院 鄭州 450002;2.鄭州大學國際學院 鄭州 450052)
樹冠通透度指基于樹木冠形的、透過樹木冠層可見天空區(qū)域的系數(shù),是評估森林健康和樹木健康的重要指標之一(Svein,2009;Spiecker,1996;宋文龍等,2015)?,F(xiàn)階段,樹冠通透度缺乏估測的標準方法,在實際作業(yè)中,觀察者通常手工繪制樹木冠層的二維形態(tài)和輪廓,估計透過冠層的可見天空區(qū)域占樹冠區(qū)域的百分比,從而實現(xiàn)等級評估(Bond,2012);也有研究人員通過參考不同通透度等級的照片(卡片),人工判讀通透度的大概程度或等級(Schomakeretal.,2007)。但這些方法均依賴樹冠形態(tài)和葉量的人工視覺判斷,受觀察者視角和主觀性影響較大,具有高度的可變性,重復性差;而且,由于每株樹木的樹冠形態(tài)和葉片分布各不相同,難以建立統(tǒng)一的判別標準。
將數(shù)字圖像處理技術用于定量評估樹冠通透度,具有快速準確、自動化程度高等優(yōu)勢,能夠克服人工判別可變性大、重復性差及易受主觀影響等局限。Zhang等(2007)采用分形維數(shù)對樹冠的空間格局進行了研究;Córica(2009)分析了陽光透過樹木冠層后在外墻上產(chǎn)生光斑的分布模式;Barthélémy等(2007)將樹冠孔隙度指標與樹冠通透度和樹木發(fā)育階段聯(lián)系起來,但是該指標依賴特定環(huán)境,并需要根據(jù)樹種進行校準;Clark等(2003)探討樹冠通透度估算方法,以“收縮包裹”的視角定義樹冠光滑邊界,提出了廣義的冠層區(qū)域和梯度變化概念,但沒有詳細說明如何設置閾值以及閾值對樹木分枝、嫩枝等的影響。
本研究基于樹木數(shù)字圖像進行單木結(jié)構特征識別與整合,不僅對樹冠所形成內(nèi)部區(qū)域的通透度進行計算,采用K-均值聚類算法將樹冠中足夠大的可見天空區(qū)域設定為大孔(通常是由于主要樹枝的缺失造成的),將較小的可見天空區(qū)域設定為小孔(大多是由于水分脅迫、葉子病變或過早退化造成的),并計算大孔、小孔的分布密度,而且還對樹冠外部區(qū)域的通透度進行研究,當可見天空區(qū)域位于樹冠輪廓外部時,探索可見天空區(qū)域凹陷進入樹冠的深度,對樹冠邊界輪廓進行平滑處理,建立樹冠的對稱凸包,計算平滑輪廓距凸包的距離和可見天空區(qū)域的凹陷深度,通過對距離設置閾值,計算深度凹陷密度?;跇涔谥械拇罂缀托】酌芏取涔谳喞纳疃劝枷菝芏?個可量化指標計算樹冠通透度。
本研究的創(chuàng)新之處在于:1) 將2個原始輪廓定義為參考形狀,即用于評估樹冠深度凹陷的對稱凸包,以便于計算小孔和大孔的樹冠平滑輪廓;2) 采用K-均值聚類算法,將計算得到的平滑輪廓上各點到凸包最近距離的集合自動分為2類,計算類間閾值(查找小的類中的最大值與大的類中的最小值,計算二者之和的均值作為閾值),通過閾值區(qū)分輕度和深度凹陷,并統(tǒng)計樹冠中的大孔和小孔;3) 根據(jù)對樹冠通透度的貢獻,賦予深度凹陷密度、大孔和小孔密度3個參數(shù)不同權重?;谝陨蟿?chuàng)新之處,計算單木通透度,以期為單木的健康狀況監(jiān)測和生長狀態(tài)分析提供技術支持。
圖像采集與預處理流程如圖1所示。
圖1 單木通透度計算流程Fig.1 Flowchart for calculating the single tree permeability
1.1 圖像采集 采用微軟公司開發(fā)的便攜式平板電腦Surface Pro 4獲取單木圖像。該設備觸控顯示屏尺寸為12.3英寸,分辨率為2 736×1 824像素,采用第6代Intel酷睿處理器,系統(tǒng)內(nèi)存為16 G,具有800萬像素后置自動對焦攝像頭,配備1 024級壓感觸控筆,電池續(xù)航時間長達9 h,機身厚度8.45 mm,質(zhì)量786 g,性能強大又不失輕盈,攜帶方便,便于移動操作,滿足本研究所需圖像拍攝和圖像處理要求。
以雪松(Cedrusdeodara)為研究對象,該樹種為常綠喬木,樹冠尖塔形,大枝平展,小枝略下垂,在我國多地栽培引種,除了可用作庭園觀賞外,也是一種重要的建筑用材。采用設備自帶的高清攝像頭獲取單木圖像,圖像中可清晰辨識枯枝、新枝等反映樹木生長狀態(tài)的細節(jié)信息。在采集圖像時,確保觀察者處于合適的距離,以相機可拍攝樹木全貌且圖像顯示清晰為準,同時盡量保持仰視角小于60°、俯視角小于30°,水平距離大于樹高,如圖2所示。標記觀察者獲取圖像時站立的位置、樹木種類、天氣情況、拍攝日期等信息,便于后續(xù)不同時段進行重復觀測,分析和評估單木通透度的變化。
圖2 拍攝單木圖像示意Fig.2 Schematic diagram of taking an image of single tree
1.2 單株樹木提取 考慮到實際采集環(huán)境,目標樹木兩側(cè)和背景中可能有建筑、車輛、相鄰粘連交叉樹木等干擾,會對單木圖像分割和提取產(chǎn)生影響,而常用的圖像分割和提取方法難以達到理想效果,故本研究利用Surface Pro 4自帶的壓感觸控筆手工圈存圖像中冠層區(qū)域,以消除圖像中的無關信息,增強單木冠層的可檢測性并最大限度簡化處理數(shù)據(jù)的復雜度,從而改進特征提取、圖像分割、匹配和識別的可靠性(劉博峰,2016)。樹冠大都包含樹枝和樹葉等詮釋樹木主要形態(tài)特征和生物學特性的組件,本研究針對樹冠進行通透度分析,排除了樹干。
采用Matlab中的函數(shù)get(axesHandle,‘CurrentPoint’)獲取鼠標點擊處的圖像坐標,并將這些坐標組成封閉多邊形。為了提取多邊形所包圍的感興趣區(qū)域(即單木區(qū)域),使用函數(shù)inpolygon(m1,n1,im,in)判斷圖像中各像素點是否在多邊形內(nèi),其中im、in為構成多邊形邊界的頂點坐標,m1、n1為圖像中的像素點坐標。get函數(shù)返回結(jié)果為邏輯類型的0或1,如果該點在多邊形內(nèi)則返回1,否則為0。運用imwrite()函數(shù)將多邊形包圍的圖像區(qū)域存儲,區(qū)域外的圖像背景填充為黑色,可實現(xiàn)手動點擊產(chǎn)生多邊形,并圈存多邊形區(qū)域,完成單木區(qū)域提取,如圖3所示。
圖3 手工圈存目標樹木Fig.3 Artificial circles target treea.原始圖像Original image;b.人工圈存軌跡Manually tracing trajectory;c.存儲樹木區(qū)域格Storage tree area grid.
1.3 感興趣區(qū)域的灰度化和二值化 將圈存的單木區(qū)域作為感興趣區(qū)域進行圖像灰度化和二值化,以便于后續(xù)處理。采用平均值法(馮志新等,2013)將RGB圖像轉(zhuǎn)換為灰度圖像:
Gray=(R+G+B)/3。
(1)
運用最大類間差法(Ostu算法)確定的閾值進行圖像二值化,該算法不受圖像亮度和對比度的影響,按圖像的灰度特性,將圖像分成背景和前景2部分(廉寧等,2014)。設灰度圖像灰度級為L,則灰度范圍為[0,L-1],運用Ostu算法計算圖像的最佳閾值為:
σ=Max{w0(t)×[u0(t)-u]2+w1(t)×
[u1(t)-u]2)}。
(2)
式中:w0為背景比例;u0為背景均值;w1為前景比例;u1為前景均值;u為整幅圖像的均值;σ取最大值時的分割閾值t即為最佳閾值T。
使用閾值T對原圖像f(x,y)進行分割,將圖像分為2部分,g(x,y)即為通過閾值處理得到的二值化圖像:
(3)
通過對圈存區(qū)域進行圖像灰度化和二值化,得到樹木二值化圖像,如圖4a、b所示。
對圖3c中圈存樹冠區(qū)域的綠色通道進行二值化,得到圖4c,并將4b與c合并,最大限度保留樹冠形態(tài),得到圖4d。
圖4 感興趣區(qū)域二值化Fig.4 Region of interest binarizationa.圖像灰度化Grayscale image;b.圖像二值化Image binarization;c.綠色通道二值化Green channel binarization;d.合并b與c Merge b and c.
將圖4d中的圈存區(qū)域取反,得到圖5a。采用Matlab中的imclose函數(shù)對圖5a進行閉運算,得到圖5b,其中閉運算采用直徑為2的圓形結(jié)構元素D1。閉運算對圖像先膨脹后腐蝕,具有填充物體內(nèi)細小空洞、連接鄰近物體和平滑邊界的作用。
采用Matlab中的bwlabel函數(shù)對圖5b進行連通區(qū)域標記,并保留面積最大的連通區(qū)域,得到圖5c。
圖5 單株樹木的二值化圖像Fig.5 Binarized images of single tree a.取反Negate;b.閉運算Closed operation;c.保留面積最大的連通區(qū)域Keep the connected area with the largest area.
1.4 樹冠平滑輪廓 首先,采用Matlab中的imfill函數(shù)對樹冠二值化圖像進行孔洞填充,如圖6a所示。然后,采用半徑為2的圓形結(jié)構元素對圖6a進行閉運算,實現(xiàn)樹冠區(qū)域的邊界平滑,運用bwboundaries函數(shù)提取樹冠邊界點,得到樹冠平滑輪廓,如圖6b、c所示。
圖6 樹冠平滑輪廓Fig.6 Smooth outline of the canopy a.孔洞填充Hole filling;b.平滑輪廓Smooth outline;c.顯示平滑輪廓Display smooth outline.
2.1 對稱軸確定 正常生長樹木的冠形在豎直方向基本保持平衡和對稱,其對稱軸通常是樹木主干所在的直線。確定樹冠在豎直方向的對稱軸,建立樹冠基于對稱軸的鏡像,構建對稱樹冠凸包,從而有助于后續(xù)與通透度相關參數(shù)的計算。
針對已經(jīng)手工圈存的單木圖像,其大小為m×n(m為行數(shù),n為列數(shù))。設定對稱軸所在的列為p,則對稱軸左側(cè)的樹木像素[1~ (p-1)列的所有樹木像素]之和與對稱軸右側(cè)的樹木像素[(p+1)~n列的所有樹木像素]之和的差值最小。在Matlab中計算出p,并在圖像中畫出對稱軸,如圖7所示。
圖7 樹冠的對稱軸Fig.7 Symmetry axis of the canopya.樹冠的對稱軸The symmetry axis of the canopy;b.在全景圖像中顯示對稱軸Show the symmetry axis in the panoramic image.
2.2 樹冠對稱鏡像和對稱凸包確立 計算樹冠通透度時,考慮樹冠中的分枝、嫩枝、新枝、葉片以及可能缺失的枝葉,以對稱軸為中心左右兩側(cè)的樹冠分別做鏡像,然后合并成為新的樹冠,如圖8b所示,白色像素為原始樹冠,灰色像素為鏡像。
圖8 對稱樹冠與凸包Fig.8 Symmetrical crown and convex hulla.樹冠Canopy;b.對稱樹冠Symmetrical canopy;c.對稱樹冠的Delaunay三角網(wǎng)Delaunay triangulation of symmetrical canopy; d.對稱樹冠凸包Convex hull of symmetrical canopy.
將圖像中表示對稱樹冠的每個像素當成一個點,每個點都具有行坐標和列坐標,所有這些點構成集合X。先用DelaunayTri函數(shù)得到集合X的三角網(wǎng),后用Convex hull函數(shù)得到樹冠凸包。
采用Matlab中的DelaunayTri函數(shù)建立約束三角網(wǎng),將集合X中的所有點進行規(guī)則化處理,使得集合中每個點位于三角形的頂點且三角形構成的網(wǎng)絡具有空圓特性和最大最小角特性(余杰等,2010;Zhanetal.,2017)。采用Matlab中的Convex hull函數(shù)獲得構成凸包的邊界點(Yangetal.,2014),依次連接凸包邊界點構成樹冠凸包,如圖8d所示。
2.3 深度凹陷密度計算 為了衡量樹冠伸展情況,采用凹陷深度表示樹冠填充凸包的緊湊程度。深度凹陷密度指樹冠平滑輪廓外部可見天空區(qū)域凹陷進入樹冠平滑輪廓較深的區(qū)域面積與凸包面積的比值。
將樹冠平滑輪廓和對稱樹冠凸包添加到樹冠區(qū)域,如圖9a、b所示。圖9c顯示了構成樹冠平滑輪廓和凸包的點,并將其中一部分(黑色矩形框)放大顯示。從平滑輪廓中黑色圓點處開始沿平滑輪廓順時針方向行走,計算平滑輪廓上各點到凸包的最短距離,并畫出距離圖,如圖10a所示。
圖9 樹冠的對稱凸包及平滑輪廓Fig.9 Symmetrical convex hull and smooth outline of canopya.平滑輪廓和凸包Smooth outline and convex hulls;b.原始圖像上顯示平滑輪廓和凸包Show smooth outline and convex hull on the original image;c.放大顯示局部的平滑輪廓和凸包Zoom in to show local smooth outline and convex hull.
采用K-均值聚類算法(顏佩等,2017)確定深度和輕度凹陷閾值,當凹陷深度超過此閾值即設定為深度凹陷。對距離圖進行K-均值聚類(K=2)處理獲得閾值Ddepression,如圖10a中紅色線段所示。
將位于樹冠平滑輪廓以外、凸包以內(nèi)的區(qū)域作為凹陷區(qū)域,如圖10b中的白色區(qū)域所示。計算該凹陷區(qū)域內(nèi)像素點到凸包邊界的最短距離,如果最短距離大于Ddepression,則該點為深度凹陷點,所有深度凹陷點構成深度凹陷區(qū)域。遍歷凹陷區(qū)域得到所有深度凹陷點,如圖10c中灰色區(qū)域、圖10d中紅色區(qū)域所示。
圖10 深度凹陷區(qū)域確定Fig.10 Establishment of the deep depression areaa.平滑輪廓上各點到凸包的最短距離The shortest distance from each point on the smooth outline to the convex hull;b.白色區(qū)域為凹陷The white area is the depression;c.灰色區(qū)域為深度凹陷The gray area is the depth depression;d.紅色區(qū)域為深度凹陷The red area is the depth depression.
2.4 大小孔檢測及密度計算 將樹冠平滑輪廓內(nèi)的區(qū)域(狹義的樹冠區(qū)域)進行連通區(qū)域標記,采用K-均值聚類算法(設定為2類)將標記的連通區(qū)域分為大孔和小孔,所有標記為小孔的連通區(qū)域像素總數(shù)除以樹冠平滑輪廓所包圍的像素總數(shù)得到小孔密度,然后計算大孔密度。
采用Matlab中的連通區(qū)域標記函數(shù)bwlabel對樹冠中所有可見天空區(qū)域進行標記,統(tǒng)計其像素數(shù),并對所有連通區(qū)域的像素數(shù)進行K-均值聚類,其中K為2(樹冠中所有可見天空區(qū)域分為小孔和大孔,其中小孔對應類別1,大孔對應類別2)。
樹冠中樹干、樹枝和樹葉的值設為0,小孔的值設為0.333 3,大孔的值設為0.666 6,樹冠輪廓外部像素的值設為1,如圖11所示。
圖11 樹冠的小孔和大孔Fig.11 Small and large holes of canopy
根據(jù)所提取的樹冠、大孔和小孔,計算得大孔密度Dmacro和小孔密度Dmicro分別為0.160 5和0.179 2。
2.5 樹冠通透度計算 參考林業(yè)專家經(jīng)驗和樹木生長規(guī)律,采用深度凹陷密度、大孔和小孔密度3個系數(shù)定量評估單木通透度Tc:
(4)
深度凹陷和大孔是因樹枝缺失造成的,考慮深度凹陷對冠形和通透度的影響高于大孔,對深度凹陷賦予更大加權。此外,式(4)中也加入了小孔對通透度的貢獻。
2.6 多視圖聚合 從0°、30°、60°、90°、120°和150°共6個角度獲取6個樹冠圖像,分別計算通透度系數(shù)Tc1、Tc2、Tc3、Tc4、Tc5和Tc6,以平均值Tca作為單木通透度系數(shù),減小因視角變化引起的通透度波動,盡可能精確反映單木真實狀況。如果難以從多個視角獲取樹木圖像,則根據(jù)實際獲取數(shù)量計算平均值作為通透度系數(shù):
(5)
3.1 系統(tǒng)實現(xiàn) 由Matlab GUI (graphical user interface,為圖形用戶界面,又稱圖形用戶接口)設計開發(fā)的系統(tǒng)界面如圖12所示,依次點擊打開圖像、提取樹木、確立對稱軸、凸包與平滑輪廓、深度凹陷、大小孔檢測等計算單幅圖像通透度,分別顯示在右下方的表格中;也可在提取冠層后,直接點擊“顯示”按鈕計算通透度。
圖12 系統(tǒng)界面Fig.12 System interface
本研究對其他雪松進行圖像采集和通透度分析,如圖13所示。
圖13 實例Fig.13 Example
3.2 樹冠對稱鏡像和對稱凸包確立 為了對本研究所提出方法進行試驗驗證和精度分析,設計一個單木驗證模型,如圖14a所示。
圖14 單木驗證模型Fig.14 Single tree verification model a.本研究所設計的單木驗證模型The single-wood verification model designed by this research institute;b.人工判讀識別驗證模型并標示大孔、小孔及深度凹陷區(qū)域Manual interpretation to identify the verification model and mark the large holes,small holes and deep depression areas;c.將驗證模型顯示在液晶平面顯示器 上Display the verification model on the LCD flat panel display.
本研究所提出方法對尺度不敏感(目標樹等比例縮放對研究結(jié)果無影響),計算出的樹冠通透度為一比值,故采用AutoCAD 2010結(jié)合Photoshop CC 2017設計單木驗證模型(尺寸無需單位)。組成樹冠3個等腰直角三角形的腰長為4,從上到下垂直排列;樹冠中大孔是邊長為1的正方形,小孔是腰長為1的等腰直角三角形。采用人工方法確定該模型凸包,并人工計算深度凹陷閾值,畫出深度凹陷區(qū)域,如圖14b所示。計算樹冠中小孔、大孔、樹冠區(qū)域、凸包和深度凹陷區(qū)域的面積分別為3、3、24、40和6.514 8,理想情況下的小孔密度、大孔密度、深度凹陷密度和通透度分別為0.125 0、0.125 0、0.162 9和0.264 6。
將單木驗證模型顯示在液晶平面顯示器上,屏幕尺寸為23.8英寸,分辨率為1 920×1 080像素,垂直和水平可視角度均為178°。采用Surface Pro 4拍攝該顯示器,獲得單木模型圖像。利用本研究所提出方法計算的小孔密度、大孔密度、深度凹陷密度和通透度分別為0.117 8、0.124 1、0.164 0和0.258 6,由此得到本研究所提出方法的精度為97.73%(0.258 6 / 0.264 6)。
3.3 深度凹陷密度計算 采用K-均值聚類算法自動實現(xiàn)大孔和小孔分類以及深度和輕度凹陷分類。比較2個類的相似度,滿足以下條件,則認定2個類可以合并:
Mhigher<1.4×Mlower。
(6)
式中:Mhigher為較大的類均值;Mlower為較小的類均值。
當進行凹陷深度計算時,如果出現(xiàn)2個類合并,則認為樹冠只呈現(xiàn)輕度凹陷,此時深度凹陷對通透度系數(shù)的貢獻為零。同樣,大孔和小孔均值滿足式(6),則全部顯示為小孔,此時大孔對通透度系數(shù)的貢獻為零。
如果多角度采集單木圖像時背景存在黏連、遮擋等問題,導致難以差距6個角度的單木圖像,則式(5)中Tca應以實際采集的圖像數(shù)量為準。采用手工圈存樹冠區(qū)域方法,未能做到全自動分析,但可減小樹冠周圍復雜背景的干擾,且人工工作量小。如果被檢測樹木的背景為其他樹木或暗色干擾物時,在人工圈存樹冠區(qū)域后,可采用He等(2011)提出的基于暗通道的去霧算法進行圖像增強,從而提高對冠層大孔和小孔的識別效果。
為解決現(xiàn)有單木通透度評估受觀察者視角和主觀性影響較大、難以建立統(tǒng)一判別標準的問題,本研究采用Surface Pro 4獲取單木圖像,進行單木結(jié)構特征識別與整合,通過建立樹冠對稱凸包和樹冠平滑輪廓,不僅對樹冠所形成的內(nèi)部區(qū)域進行研究,標記內(nèi)部的可見天空區(qū)域,并統(tǒng)計各連通區(qū)的像素數(shù),采用K-均值聚類算法將各連通區(qū)域分為大孔和小孔,而且還對樹冠輪廓外部可見天空區(qū)域凹陷進入樹冠的深度進行探索,采用K-均值聚類算法判別輕度和深度凹陷。通過自動量化大孔和小孔密度、深度凹陷密度3個指標計算單木通透度系數(shù),并對過程和結(jié)果進行詳細分析。結(jié)果表明,本研究所提出方法的精度高達97.73%,測量速度快、人工工作量小,可為單木健康狀況監(jiān)測和生長狀態(tài)分析提供技術支持,同時該研究思路和方法也可以推廣應用到其他樹木和作物的監(jiān)測分析,具有一定的實用價值。