利用輻射度模擬虛擬枇杷果園冠層光分布

楊怡斐，唐麗玉，崔 磊，陳舒煒

YANG Yifei1，2,TANG Liyu1，2,CUI Lei1，2,CHEN Shuwei1，2

1.福州大學 空間數據挖掘與信息共享教育部重點實驗室，福州 350116

2.福州大學 地理空間信息技術國家地方聯合工程研究中心，福州 350116

1.Key Laboratory of Spatial Data Mining&Information Sharing of MOE,Fuzhou University,Fuzhou 350116,China

2.National Engineering Research Center of Geospatial Information Technology,Fuzhou University,Fuzhou 350116,China

1 引言

果樹形態結構決定了冠層內的太陽光的分布情況，而光分布影響了光合作用效率和冠層內微環境，進一步影響果實產量和質量。目前冠層內光環境量化分析方法主要有地面實測法[1]、數學模型法[2]和三維植物結構模擬法[3]。地面實測法受實際客觀條件限制，空間分辨率和時間分辨率有限；數學模擬法一般對真實冠層進行簡化或某種假設，采用一些函數描述葉傾角分布和方位角[4]，在遙感反演得到較多應用，但與真實冠層還有較大誤差。基于三維冠層結構的模擬方法已受到學者們的青睞，三維冠層結構可以更精確地描述冠層組分之間的拓撲結構，可以采用更精細面元表達冠層的三維信息。

輻射度算法和光線跟蹤算法是三維植物冠層模擬輻射傳輸的主要方法。例如，文獻[5]把冠層三維結構離散化為小立方體集合，同時使用高斯迭代法為每個立方體賦予光屬性的方法模擬離散的光線集合。文獻[6]利用嵌套輻射度方法，使用蒙特卡洛框架模擬自然光在植物冠層的三維分布。以上模擬已達到植物器官級別或更小單元，因此這類模型可以模擬冠層內任一點的冠層光分布或光強。文獻[7]提出了RAPID模型，可以模擬復雜場景中任意的植被冠層，進而快速計算二向反射因子（BRF），適用于不同分辨率的衛星數據。文獻[8]從葉片、枝梢、樹枝以及整株四個尺度上研究樹形結構對冠層光截獲的影響，并引入評價光截獲效率的STAR（Silhouette to Total Area Ratio）值，即冠層或枝梢葉片被太陽照射的葉面積與總葉面積之比。文獻[9-11]分別以黃瓜、玉米、煙草冠層為例分析了植物冠層內的光傳輸。近年來輻射度模型應用于植物冠層內光分布模擬多見于單株尺度，而果園尺度的研究對象為果樹群體，由單株走向群體是模型應用中的必然步驟。在植物群體尺度上，一些學者采用光線跟蹤方法[12-13]模擬光分布，文獻[12]未考慮冠層內部相互反射（如枝條、果實等）、葉片透射等情況，模擬蘋果樹群體冠層光分布；文獻[13]利用了開源光線跟蹤模擬平臺GroIMP，且用遞歸方法計算光線在櫸木群體冠層內的多次反射。以上兩篇文獻的不足之處在于未能考慮植株所處地形的反射光對冠層光分布的影響。另一些學者在此基礎上，加入了植物周圍環境對群體冠層內光分布的影響，如文獻[14]以玉米為例，利用輻射度-圖形學結合模型（RGM）計算玉米群體冠層內光分布情況，且考慮了玉米葉片、莖稈面元的反射率與透射率以及地形的影響。不足在于假定地形為平面，與實際情況不符。文獻[15]在溫室內人工光環境條件下，綜合考慮太陽直射、天空散射、人工光源模擬了現代月季群體冠層內光分布與光截獲，并計算了溫室內地板、墻壁與屋頂的反射光對群體冠層光分布的影響，但其研究對象冠層結構簡單，未能應用于復雜的果樹冠層。

目前植物群體光環境模擬，大多未考慮實際植物所處地形表面的起伏對光分布的影響。相關研究表明，地形因子對BRF分布有著非常明顯的影響[16]。地形是果樹的主要載體，地面反射光在冠層內的分布情況是研究虛擬枇杷果園冠層光分布必不可少的組成部分。構建地形格網，剖分為三角形面元參與輻射計算，可進一步提高光輻射模擬的準確度。

本文將采用輻射度算法模擬顧及地形影響的枇杷果園的光分布。

2 枇杷果園場景構建

2.1 數據獲取

構建精細枇杷幾何模型是進行冠層光分布模擬的基礎。本研究團隊于2016年1月在福建省云霄縣溪口村（117°7′，24°8′）選取若干枇杷樹作為調查對象，樹齡為3年。調查與測量的主要內容有：（1）使用桿尺、皮尺等進行枇杷形態結構參數測量，具體參數有樹高、葉幕層厚度、干高、冠幅、各級枝條分枝長度、周長、角度、分枝數等；（2）使用Artec Eva手持三維激光掃描儀掃描枇杷葉簇獲取葉器官點云數據；（3）使用RIGEL VZ400地面三維激光掃描儀獲取枇杷果園地形點云數據。

2.2 地形網格生成

使用VZ400地面三維激光掃描儀獲取的地形點云數據密度較大，且分布均勻，適合于反距離權重插值法（Inverse Distance Weighted，IDW）的應用條件[17]。假設未知地形點z0受近距離已知地形點的影響比遠距離已知地形點的影響更大，z0的高程為：

式中zi為已知點i的高程；di是點i與點0之間的距離；s為用到的已知點的數量；k是確定的冪，在本例中k選值為2。IDW插值的一個重要特性是所有預測點高程值都位于已知點高程值最大最小值之間，插值后地形點高程不會超出原始地形點云覆蓋的高程范圍，且分布均勻，便于構造三角格網。地形格網生成過程如圖1所示。

圖1 地形格網生成過程

2.3 單株木冠層模型構建

利用葉器官點云數據，使用Geomagic Studio與3ds Max軟件構建枇杷葉片精細三維模型。根據葉片大小不同，每個葉片約由80～200個三角面元構成。為了模擬枇杷樹葉真實的組織狀態，由葉簇掛接于枝干系統，形成虛擬枇杷樹模型，每個葉簇約由3～6片葉片構成，包含的三角面元數約為500～1 000個，如圖2所示。

采用研究團隊自主研發的ParaTree參數化單樹建模軟件[18-19]，輸入調查獲得的各項枇杷形態結構特征參數，建立枝干系統，在枝系統上掛接葉器官模型，構建精細枇杷單樹幾何模型，如圖3所示，模型總面元數為141 217個。

圖2 葉簇構建過程

圖3 單樹枇杷幾何模型

2.4 虛擬果園場景生成

虛擬果園場景抽象為由三維地形和三維冠層模型組成。本研究樣地的地形呈階梯狀，果樹栽種于每一壟。根據實地調查數據，可以計算每株果樹種植的地理坐標。精細單樹枇杷幾何模型數據結構中包含了組成樹模型不同器官（包含枝干、葉片等）的大量三角形面元的索引數組，索引數組對應于坐標數組，實現了一一對應關系。因此在虛擬枇杷果園場景中只需使用實例化技術調用樹模型對應的類的構造函數，再結合枇杷樹種植點位，即可在不同種植點位“種植”枇杷樹模型，如圖1（4）所示。此時，樹模型的坐標數組并未根據“種植”點位不同而存儲多次，僅在繪制階段調用三角形面元索引，獲取對應的坐標值，做相應的矩陣平移操作，節約了大量的內存開銷，保證了三維果園場景的流暢顯示。

3 枇杷果園冠層內光三維分布模擬

3.1 技術流程與實現

果園的光源主要包括太陽直射光和天空散射光。虛擬枇杷果園冠層內光分布模擬總體技術流程如圖4所示。

首先根據研究樣地的地理經緯度以及模擬的日期、時刻，利用天文計算公式[20]估算樣地所在地的太陽高度角和方位角，進而確定光源位置、光線方向、光源面片初始攜帶能量等參數。對于太陽直射光源，在光源位置構建光投放平行四邊形并剖分為多個三角形面元，模擬平行直射光入射植物冠層。而散射光源則構建一個包圍整個果園場景的天空球，將此天空球按照經緯度方向剖分為多個三角形面元，從而模擬均質散射光從多個方向離散射入植物冠層。根據枇杷幾何模型與光源幾何信息，采用三維離散視角因子方法[21]（3D-Discrete View Factor，3D-DVF）計算輻射度模型形狀因子。求解輻射度模型，從而計算冠層內太陽直射光分布、天空散射光分布以及地形網格反射初始能量值。將地形網格作為反射光源，與枇杷幾何模型再次構建輻射度模型并進行求解，計算枇杷果園冠層內地面反射光分布。將每個面元獲得的直射光、散射光、地面反射光的輻射強度進行累加，即得到枇杷果園冠層內每個面元的輻射強度。

圖4 虛擬枇杷果園冠層光分布模擬總體技術流程

3.2 虛擬果園場景組織

實地調查樣地枇杷樹總量為12棵，按每棵樹約由14萬個三角形面元構成，場景數據量約為170萬個面元。以計算地面反射光在冠層內分布為例，為保證精度，將地形面元剖分為1 250個三角形面元，總計算量超過21億次以上，系統內存消耗巨大，因此高效組織管理虛擬場景非常重要。為了實現樹面元與光源面元間的遮擋判斷，采用均勻體素剖分場景。根據場景樹模型，建立場景最小包圍盒，將包圍盒劃分為64×64×64均勻體素網格。針對場景中每個樹面元做歸屬判定，若某樹面元三角形中心點位于某體素內，則認為該面元與此體素具有唯一相關性。將平面直線的Bresenham算法推廣至三維條件下，設某空間直線段兩端點坐標為(x1,y1,z1)與 (x2,y2,z2)，令 dx=|x1-x2|，dy=|y1-y2|，dz=|z1-z2|，將該空間直線段投影至xoy平面后即可根據平面直線的Bresenham算法求出x坐標與y坐標；同理可求出投影至xoz平面上的x坐標與z坐標，據此可求出該直線與空間體素相交的任意三維點坐標。

光源面元發射出的光線在包圍盒內遮擋判斷，光線原點為光源三角面元中心點，方向對應于樹面元中心點，除每條光線與包圍盒相交的第一個體素外，其余相交體素內的樹面元均視為被遮擋。CUDA編程模型提供了基于線程網格與線程塊的雙層結構。采用一維線程格網與一維的線程塊，一個線程內計算單個樹面元與單個光源面元間的遮擋情況，一個線程塊內計算單個樹面元與所有光源面元間的遮擋情況。基于三維體素遍歷的并行遮擋判斷流程如圖5所示。

圖5 遮擋判斷流程圖

3.3 形狀因子與輻射度求解

形狀因子又稱視角因子，是僅與面元在場景中所處位置有關的無量綱常數，面元Si與面元Sj間的形狀因子可表示為：

式中HID(dAi,dAj)定義為遮擋系數；θi和θj是面元Si和面元Sj法向量與兩面元間夾角；rij為面元Si與面元Sj間的距離。采用3D-DVF法求解輻射度模型，具有普適性高、可并行實現的特點。該方法可表示為：

使用CUDA編程模型實現形狀因子與輻射度求解，需要迭代3次式（3），第一次將地形面元視為果園場景一部分，與光源面元進行輻射計算。采用一維線程格網與一維線程塊雙層結構，一個GPU線程計算一對樹面元（包含地形面元）與光源面元間的形狀因子能量，一個線程塊計算一個樹面元（包含地形面元）與所有光源面元間形狀因子與能量，據此求出太陽直射輻射強度[22]。同理可求出天空散射輻射強度。將地形面元接收到的直射與散射輻射能量乘以地面反射率，取值為0.20[23]，作為地面反射光源，與所有樹面元進行輻射計算，CUDA線程結構與前兩次類似。

4 分析與討論

以2016年1月20日于福建省云霄縣溪口村樣地實地調查所獲數據為依據，枇杷樹葉無凋落現象，不考慮葉片光透射的情況。估算冠頂輻射強度時，假設研究樣地天氣晴朗無云。模擬從早上8時開始，間隔為1 h，選取一天中的10個時刻，依次模擬果園群體冠層的光分布。本案例模擬太陽光在400～700 nm波段的輻射，即為光合有效輻射（Photosynthetically Active Radiation，PAR）分布。

果園場景模擬以三角形面元為最小單位，模擬結果可以輸出果園冠層每個葉片的每個三角面元截獲的PAR能量值，從而可以計算冠層在某時刻獲得的PAR能量，亦可進行三維可視化輸出。如圖6（1）所示為樣地正午12時冠層光分布模擬效果圖，在接受太陽直射、天空散射與地面反射光源輻射后，模型中顏色越亮的部分代表到達該面元的光輻射能量越多，PAR強度越強。

根據研究樣地地形數據和樹木的分布情況，構建12株與16株兩種不同種植密度的虛擬枇杷果園，其中12株情況為樣地真實果樹位置，16株情況為增加種植密度后的對照組。分別研究太陽直射、天空散射、地面反射以及總PAR輻射強度。單位面積葉獲得PAR瞬時能量值代表該時刻枇杷果園冠層對光的截獲能力，以及枇杷植株對光的利用率，進而可判斷該種植密度是否合理。

圖6（2）中早8時12株密度下獲得平均太陽直射PAR強度高于16株密度，是因為早晨太陽高度角較低，光線從植株側面進入冠層，植株間遮擋明顯，更大的種植密度的情況下會有更多的葉片相互遮擋；在正午時刻，光線從冠頂進入冠層，這兩種種植密度的冠層內部葉片相互遮擋情況類似，且3年生枇杷還未郁閉，植株間的遮擋不明顯，因此，此時太陽直射瞬時PAR能量接近。圖6（3）中12株密度下天空散射瞬時能量值在各個時刻均高于16株密度，且差值相近，是因為不同時刻下散射光源天穹發出的光線方向一致，僅有輻射能量大小不同，枇杷果園冠層內葉片相互遮擋情況幾乎一致。圖6（4）和圖6（5）表明果園冠層葉片獲得經過地面反射后的光輻射能量約為地形網格初始反射能量的1/3，且正午時刻最高截獲能量速率為35 μmol·m-2·s-1。結合圖6（6）分析，同一時刻冠層光截獲地面反射約占冠層總光截獲的1/20～1/10。兩種不同種植密度下的虛擬枇杷果園在一天10個時刻內的單位面積的PAR變化趨勢一致，早晨和傍晚太陽輻射強度弱，因此冠層獲得的PAR值相應變小。對一天10個時刻的模擬值進行平均，12株密度下果園獲得的平均PAR值為236.4 μmol·m-2·s-1，16 株密度下則為 224.8 μmol·m-2·s-1，兩者相差約為11.6 μmol·m-2·s-1。模擬結果表明，對于3年生的枇杷植株，在樣地的生境下，這兩種種植密度均未達到郁閉，株間遮擋對光分布影響不是很大。

圖6 果園冠層內光分布分析

5 結束語

目前果樹冠層光分布模擬多數停留在單株尺度，針對這種研究現狀，同時考慮地面反射對冠層內光截獲的影響，設計并實現了利用輻射度的枇杷果園冠層光分布模擬方法。利用地面激光雷達掃描獲取的地形點云數據構建果園地形網格，手持激光掃描獲取的枇杷葉器官點云數據構建葉簇實體三維模型，然后采用單樹建模軟件ParaTree建立枇杷冠層模型，最后根據植株的分布規則，形成三維果園場景。使用3D-DVF方法求解輻射度模型，開展不同種植密度下果園冠層光分布模擬，可定量化計算冠層內每個組分獲得的光的強度，且可直觀表達光在冠層內空間分布情況。模擬結果表明，果園冠層內光強度分布日變化呈近似拋物線，早晚較小而正午到達最大值；果園冠層截獲的地面反射PAR能量主要作用于冠層內部與下層，這部分光能對于果園內枇杷植株生長發育的作用不容忽視。

參考文獻：

[1]高照全，趙晨霞，李志強，等.我國4種主要蘋果樹光合能力差異研究[J].中國生態農業學報，2013，21（7）：853-859.

[2]Baldocchi D，Hutchison B.On estimating canopy photosynthesis and stomatal conductance in a deciduous forest with clumped foliage[J].Tree Physiology，1986，2：155-168.

[3]Tang L，Hou C，Huang H，et al.Light interception efficiency analysis based on three-dimensional peach canopy models[J].Ecological Informatics，2015，30：60-67.

[4]Ross J.The radiation regime and architecture of plant stands[J].The Quarterly Review of Biology，1982，71（3）：344.

[5]Gastellu-Etchegorry J P，Demarez V，Pinel V，et al.Modeling radiative transfer in heterogeneous 3-D vegetation canopies[J].Remote Sensing of Environment，1996，58（2）：131-156.

[6]Chelle M，Saint-Jean S.Taking into account the 3D canopy structure to study the physical environment of plants：The Monte Carlo solution[C]//4th International Workshop on Functional-Structural Plant Models.Montpellier，France：CIRAD-AMAP，2004：176-180.

[7]Huang H，Qin W，Liu Q.RAPID：A Radiosity Applicable to Porous IndiviDual objects for directional reflectance over complex vegetated scenes[J].Remote Sensing of Environment，2013，132（10）：221-237.

[8]Silva D D，Han L，Costes E.Light interception efficiency of apple trees：A multiscale computational study based on MAppleT[J].Ecological Modelling，2014，290（10）：45-53.

[9]Wiechers D，Kahlen K，Stützel H.Evaluation of a radiosity based light model for greenhouse cucumber canopies[J].Agricultural and Forest Meteorology，2011，151（7）：906-915.

[10]孔婭，勞彩蓮，曹素云.利用3D模型模擬天空與葉面散射對玉米冠層截光率的影響[J].農業工程學報，2011，27（5）：248-252.

[11]黃帆，勞彩蓮，肖翠霞.基于光線跟蹤的冠層光分布模型參數研究[J].中國農業大學學報，2013，18（6）：96-101.

[12]Kang H，Fiser M，Shi B，et al.IMapple-functional structural model of apple trees[C]//2016 IEEE International Conference on Functional-Structural Plant Growth Modeling，Simulation，Visualization and Applications，China，2016：90-97.

[13]Reinhard H，Ole K，Dirk L，et al.The rule-based language XL and the modelling environment GroIMP illustrated with simulated tree competition[J].Functional Plant Biology，2008，35（10）：739-750.

[14]溫維亮，孟軍，郭新宇，等.基于輻射照度的作物冠層光分布計算系統設計[J].農業機械學報，2009，40（S1）：190-193.

[15]Buck-Sorlin G，de Visser P H，Henke M，et al.Towards a functional-structural plant model of cut-rose：Simulation of light environment，light absorption，photosynthesis and interference with the plant structure[J].Annals of Botany，2011，108（6）：1121-1134.

[16]王新云，郭藝歌，過志峰，等.地形因子對森林冠層BRF模擬的影響分析[J].遙感信息，2012（4）：82-85.

[17]劉光孟，汪云甲，張海榮，等.空間分析中幾種插值方法的比較研究[J].地理信息世界，2011（3）：41-45.

[18]Tang L，Chen C，Zou J，et al.OntoPlant：An integrated virtual plant software package for different scale applications[C]//Spatial Data Mining and Geographical Knowledge Services（ICSDM），2011：308-314.

[19]林定，陳崇成，唐麗玉，等.基于顏色編碼的虛擬樹木交互式修剪技術及其實現[J].計算機輔助設計與圖形學學報，2011，23（11）：1799-1807.

[20]王炳忠.太陽輻射計算講座：第一講 太陽能中天文參數的計算[J].太陽能，1999（2）：8-10.

[21]趙權，黃運保，孫宇航.CUDA架構下的靶丸輻射能流并行計算[J].計算機輔助設計與圖形學學報，2013，25（7）：937-945.

[22]侯璨，唐麗玉，陳崇成，等.基于并行輻射度的虛擬植物冠層內光分布模擬[J].系統仿真學報，2015（10）：2337-2343.

[23]林正云.福建省太陽總輻射和地面輻射平衡的分布[J].太陽能學報，1994（3）：248-256.