張衛(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所示。

圖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 圖像處理

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 系統(tǒng)實現(xiàn)與分析

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)提出的基于暗通道的去霧算法進行圖像增強，從而提高對冠層大孔和小孔的識別效果。

4 結(jié)論

為解決現(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)測分析，具有一定的實用價值。