基于L-系統的3D虛擬植物冠層光合作用模擬模型

陳 剛 陳 斌 林郁欣 錢基德

(1.中國科學院成都計算機應用研究所， 成都 610041； 2.中國科學院大學， 北京 100049； 3.中國科學院廣州電子技術研究所， 廣州 510070; 4.福建網龍網絡公司研發中心， 福州 350001)

0 引言

植物冠層光合作用理論是作物生長動力學模型建立的基礎。由于植物冠層內葉片遮擋，導致光輻射(PAR)分布難以直接計算。JOHNSON等[1]采用Monsi-Saeki理論對冠層內葉片對光截獲進行計算，結合光合環境因子完成冠層內光合作用逐日同化物累積量估算；于強等[2]、李艷大等[3]采用經驗統計公式對作物的冠層進行描述；李緒孟等[4]基于數學方法進行水稻虛擬切層的構建，建立了水稻群體冠層內光分布及光合作用速率模型。

隨著虛擬植物技術的發展，基于L-系統[5]、參考軸Amap[6]、雙尺度自動機[7]等能夠直觀生成3D植物模型。在準確模擬植物葉片數量、葉片傾角和方位角分布規律基礎上，可構建特定種類、大小的植物虛擬冠層，將其用于植物冠層PAR分布及光合速率研究。一般采用計算機光線跟蹤和輻射度算法模擬PAR在冠層內的分布。中國農業大學研究人員馬韞韜、勞彩蓮[8]、王錫平[9]在植株測量的3D結構體中，基于光線跟蹤、天空可見率的方法進行冠層內的PAR傳輸模擬，并在儀器研制、算法設計、實驗驗證等方面做了較為完整的工作。近年來，曹秉帥[10]通過構建水稻功能-結構耦合模型G-ORY，研究了兩個不同地理位置的3個水稻株型冠層的PAR分布差異；石佩[11]利用三維數字化儀器構建不同冠層形狀的桃樹模型，分析測定了紡錘形和開心形桃樹冠層內不同區位的PAR、冠層凈光合作用速率，得到若干桃樹冠層光合作用規律；李艷大等[3]以雙季稻為研究對象，假設冠層內同一水平面輻射強度相同，構建水稻冠層光合作用速率與干物質積累模型，并基于實測數據驗證了模型精度。當前，基于虛擬植物光合作用模型的研究現狀及不足在于：①以虛擬植物技術代替三維數字化儀，在滿足結構模型精度條件下，可通過快速改變植物結構模型，進行冠層光合作用能力變動估測。②已有研究關注虛擬植物PAR分布模擬較多，較少研究關注對植物進行長時間內多次冠層PAR分布模擬，結合實測的光合環境因子(CO2濃度、大氣濕度)計算冠層光合作用速率，以及周期內干物質積累模型。③基于虛擬冠層的光合作用模型，較少研究將真實環境的實時光照強度、光合環境因子(氣溫、濕度、CO2濃度)作為輸入，進行冠層光合作用速率計算。

本文采用L-系統構建3D虛擬植物拓撲結構，基于概率分布曲線動態調整冠層內葉片分布。基于OpenGL圖形庫繪制虛擬植物冠層模型，利用GPU加速算法對冠層內的PAR輻射傳輸及分布模擬進行快速計算，同時進行太陽光角度、CO2濃度、空氣濕度等光合因子的逐時模擬，計算虛擬冠層光合作用速率，進行模擬周期內生物量積累的估算。

1 3D虛擬冠層構建

植物冠層結構是影響冠層內太陽光合有效輻射分布的最重要因素。虛擬植物技術能較好地刻畫重建植物3D虛擬冠層結構，劉剛等[12]對農林作物三維重建方法進行了研究綜述。虛擬植物結構模型包括拓撲模型和幾何結構模型。拓撲結構模型主要用于描述植物離散結構單元之間的連接關系，并可以通過抽象的語言反映出植物的生長階段，經典方法包括L-系統、參考軸、雙尺度自動機等方法。幾何結構模型是對虛擬植物模型基本組成單元的三維信息描述，包括冠層葉片的密度及分布(傾角、方位角)。中國科學院自動化所中法實驗室、中國農業大學資源環境學院開展了大量虛擬植物建模、植物環境交互模擬、植物功能-結構模型的研究，研制出了GreenLab[13]，開展了多種林木、農作物的模擬研究[14-16]。構建一個高精度的3D植物冠層，主要是枝干拓撲構建和葉片分布模擬。本文選取針葉型杉木冠層作為模擬對象。

1.1 枝干模擬

依據杉木枝干特征，利用L-系統規則構建各個生長階段的幼齡杉木型枝干模型。基于L-系統的拓撲結構模型構建流程為：①依據杉木生長過程中拓撲變化規律，設計公理和產生式算法。②設計各分枝產生式符號、分枝角度、長度、半徑、產生概率等。③設計產生式約束條件。④程序對規則進行解析、字符迭代、繪制，生成三維樹木模型。

圖1 用于幼齡杉木虛擬冠層構建的L-系統規則Fig.1 L-system rules for generating virtual canopy of young Chinese fir

幼齡杉木3D虛擬冠層生成對應的L-系統規則如圖1所示。規則解釋為：由產生式A(l,w)開始推導，主軸A的頂點在每一步的推導中產生1個節間F和4個側枝B。4個側枝為杉木的一輪模擬枝條。產生式B(l,w)表示側軸B的頂點產生1個節間F和兩個側枝C。兩個側枝之間有一個節間F，為模擬側枝上的不嚴格輪生(輪生枝條的著生點不在主干的同一高度)。側枝C為一個節間F，迭代步長設置12。

1.2 葉片分布模擬

葉子是光合作用的主要器官，葉子的位置和方向是影響光分布的主要因素。杉木葉片呈線狀披針形。為了真實地表示葉片，采用3DS MAX建模軟件進行單葉曲面建模，獲得葉片骨架點的坐標和點位間的構網索引。在L-系統解析軟件中讀取點位坐標和索引，繪制出精細單葉模型。單葉大小可在規則中交互輸入數值進行控制。葉片在枝干上的分布參數包括葉片的密度、葉片傾角、葉片方位角。根據觀測數據，確定杉木的葉傾角分布范圍，使葉片隨機分布于該角度范圍內。單葉和葉片在枝干上的分布模擬如圖2所示。

圖2 幼齡杉木針形葉片及其在枝干上的分布模擬Fig.2 Simulation of young Chinese fir needles and their distribution on branches

1.3 3D杉木虛擬冠層模型描述

以圖3b中5年生杉木為研究對象，按照1.1節與1.2節的虛擬冠層生成方法，模擬生成3D虛擬杉木冠層結構：模型高度為3.30 m，冠幅為1.73 m，模型中單葉面積為2.25 cm2，整個虛擬冠層葉片個數為32 081，葉總面積為7.21 m2。

圖3 杉木模型與杉木照片對比Fig.3 Comparison between Chinese fir model and actual picture

2 基于虛擬冠層模型的PAR輻射傳輸模擬

2.1 虛擬冠層頂部太陽光直射方向與強度模擬

真實世界植物冠層頂部太陽直射方向和太陽光強度隨地理位置、時間、天氣狀況而變化。根據天文學、大氣相關公式和模型，分別對冠頂太陽直射光的方向、冠頂PAR強度進行模擬。

2.1.1太陽光直射方向計算

根據文獻[17]，經緯度地理位置不同時刻太陽光直射方向(高度角和方位角)的計算流程如圖 4所示，三維可視化模擬結果如圖5所示。

圖4 太陽光直射方向模擬流程Fig.4 Flow chart of sun direct radiation

2.1.2太陽直射和散射光強計算

到達植物冠頂的輻射包括太陽直射和天空散射兩部分。

2.1.2.1直射光強模擬

太陽光的傳輸過程分為2個階段：從太陽到達地球大氣層上界面和太陽光在大氣層中的傳輸。

(1)太陽光到達地球大氣層上界面的傳輸過程

太陽到達地球大氣層上界面的太陽光強與日地距離成反比。對于垂直于太陽光線的平面，其輻射強度計算公式為

圖5 不同時刻太陽直射三維可視化模擬結果Fig.5 Direct sun direction simulation results

(1)

式中E0——太陽輻射強度常數，世界氣象組織推薦值為1 367 W/m2

(r0/r)2——日地距離訂正系數

r0——日地平均距離

r——當前日地距離

(2)太陽光在大氣層中傳輸

垂直于太陽光線方向的太陽直射輻射強度與其穿過大氣層的路徑和大氣透明度有關。大氣量與太陽高度角和當地的地形高度以及大氣壓有關。一定地形高度下的大氣量計算公式與海平面大氣量的計算公式為

Mh=M0ph/p0

(2)

其中

M0=[1 229+(614sinδ)2]1/2-614sinδ

(3)

ph/p0=[(288-0.006 5h)/288]5.256

(4)

式中Mh——一定地形高度下的大氣量

M0——海平面上大氣量

δ——太陽高度角

ph/p0——大氣壓修正系數

ph——一定海拔高度下大氣壓

p0——標準大氣壓

h——地形高度

根據文獻[17]提出的晴朗無云條件下的太陽直射大氣透明系數計算公式為

τb=0.56(e-0.56Mh+e-0.095Mh)

(5)

式中τb——太陽直射大氣透明系數

在考慮大氣衰減情況下，垂直到達地表的太陽直射光強計算公式為

Ib=E0(r0/r)2τbsinδ

(6)

式中Ib——太陽輻射到達地表的太陽直射光強

2.1.2.2散射光強模擬

在晴朗無云條件下天空散射是一個均勻散射，在直接輻射和散射輻射之間存在線性關系[18]。散射光強計算公式為

Is=E(r0/r)2τd

(7)

其中

τd=0.271-0.294τb

(8)

式中τd——散射輻射透明系數

Is——太陽輻射到達地表的天空散射光強

2.1.2.3光合有效輻射PAR計算

光合有效輻射(PAR)部分的光能約占太陽總輻射的40%，單位為μmol/(m2·s)，其同太陽光強(單位：W/m2)的換算系數隨光強不同而改變。根據文獻[19]將單位W/m2轉換到單位μmol/(m2·s)的轉換系數設為4.6。

2.2 虛擬冠層光合有效輻射傳輸模擬

采用光線跟蹤算法和天空可見率算法分別計算太陽直射PAR和散射PAR在虛擬冠層內的空間分布，精確獲取某一時刻虛擬冠層內葉片上每一個三角面元的光強，用于植物葉片光合作用速率的計算。

2.2.1太陽直射PAR在虛擬冠層內的分布模擬

采用正向光線跟蹤算法，從光源向場景發射足夠數量的光線，讓光線與虛擬植物冠層相交，并繼續跟蹤光線與冠層組分相交后產生的反射和投射光線，直至光線的光強足夠小或者逸出場景為止。為提高光線跟蹤算法效率，將虛擬冠層進行空間剖分，建立幾何包圍盒包圍葉片，使用八叉樹進行空間索引。根據太陽方位角、高度角、太陽直射能量等確定光線跟蹤投射的光線方向，單條光線攜帶的能量等；采用遞歸方法實現跟蹤光線同葉片是否相交的計算。根據直射光線與葉片所在包圍盒相交后，按設定的投射和反射能量分量，進行正向光線跟蹤，直到光線能量衰減到設定閾值內。技術流程圖如圖 6示，3D模擬結果如圖7所示。

圖6 太陽直射光在虛擬冠層內的輻射傳輸模擬流程圖Fig.6 Simulation flowchart of radiation transmission of direct radiation in virtual canopy

圖7 虛擬冠層模型空間剖分及光線投射位置Fig.7 Space segmentation of virtual canopy and ray plane location

2.2.2太陽散射PAR在冠層內分布模擬

采用龜形算法來計算虛擬冠層中每個葉片對天空散射的透過率。葉片的天空透過率可表示為以葉片上的三角面元的中心點向該點上半球空間望去，不受虛擬冠層內部其他結構組分遮擋，而可以直接可視到天空的立體角范圍所占整個上半球天穹的比例。對虛擬冠層中的一個葉片，計算其天空散射透過率的流程圖如圖8所示。

圖8 基于虛擬冠層的PAR散射模擬技術流程圖Fig.8 Flow chart of PAR scattering simulation based on virtual canopy

虛擬半球天穹被投影三角形遮擋的計算流程為：計算出投影三角形3個頂點在半球穹頂上所處的高度角和方位角，這3組高度角和方位角可確定一個經度方向和緯度方向上的跨度范圍，可縮小覆蓋格網的范圍，加快遍歷。其中，若已知投影三角形某個頂點坐標為(x,y,z)，則可以計算該頂點所處位置的高度角α為

(9)

方位角β計算公式為

(10)

依據式(10)求得β位于-90°～90°范圍內，根據x和z的正負確定β所在的象限，計算出方位角。求出實際的高度角α和方位角β后，即可確定3對高度角、方位角中的最大、最小值，將這組最大值和最小值作為查詢的區間對相應高度角和方位角范圍內的格網進行遮擋判斷。大大減少了不必要的遮擋判斷。圖9為某面元天空可見率的模擬結果示意圖。

圖9 天空可見率模擬結果示意圖Fig.9 Sketch of sky visibility calculation

圖9中亮黃色三角形為計算天空可見率的三角面元，中間紅色三角形為虛擬冠層中某個遮擋三角形，頂部紅色三角形為遮擋三角形在穹頂上的投影。半球天穹上被遮擋區域的天穹格網點不予以繪制，白色格網點代表天空立體角可見部分。由此，計算天空可見立體角格網占總格網數量的百分比，即天空散射PAR的三角形面元上的天空可見率。

某一地理位置的經緯度根據某一時刻可計算冠層頂部太陽光方向和PAR值，經上述模擬流程可知虛擬冠層內每個三角面元的PAR值。使用顏色梯度代表PAR值。對上文所建的幼齡杉木進行PAR輻射傳輸模擬，模擬結果如圖10所示。

3 虛擬冠層光合作用模型

圖10 5年生幼齡杉木太陽直射和散射輻射結果Fig.10 Direct and scattered radiation results from five-year-old Chinese fir

植物光合作用模型根據尺度可分為單葉尺度模型、冠層尺度模型、群體和區域尺度模型。上文模型可模擬計算冠層內每一個三角面元的PAR值，根據其他光合因子，可以累積冠層內每一個三角面元光合速率得到冠層光合作用模型。采用LEUNING[20]建立的光合速率與光強之間的非直角雙曲線關系為

(11)

式中 ?i——初始量子效率

P——最大光合速率

θ——曲線凸度

Pg——總光合速率

PAR——光合有限輻射強度

Pmax——最大光合速率

冠層光合作用計算流程圖如圖11所示。

圖11 逐時光合速率計算流程圖Fig.11 Flow chart of time-lapse photosynthetic environment factor simulation and canopy photosynthetic rate calculation

虛擬冠層的光合作用速率可視為冠層內所有葉片光合作用速率之和。冠層平均光合速率計算公式為

(12)

式中Pn——冠層光合作用速率

Pi——第i片葉光合作用速率

Si——第i片葉的葉表面積

n——葉片數量

量子效率是植物光響應特征的重要指標。其受細胞間隙的CO2濃度的影響比較明顯，其計算公式為

?i=?i0(Ci-T)/(Ci+2T)

(13)

式中 ?i0——植物葉片CO2同化的內稟量子效率，其與植物葉片自身光合特性相關

Ci——葉肉細胞CO2濃度，本文將葉肉細胞CO2濃度等同于葉片空氣中的CO2濃度

T——溫度

根據文獻[21]，最大光合速率計算公式為

Pmax=Vm0(Pi-G)/[Pi+Kc(1+PO2/Ko)]

(14)

式中Vm0——單位面積上Rubisco酶最大羧化能力

Kc、Ko——Rubisco酶反應中的CO2和O2米氏反應常數

PO2——葉內O2分壓

4 模擬實驗

4.1 光合同化物積累估算

對幼齡杉木模型分別進行太陽光直射和天空散射模擬計算，獲取冠層內葉片三角面元的PAR值。設模擬周期為7月1日—12月30日。根據2.1節所述虛擬冠層太陽輻射時刻點計算方法求得各輻射點的時刻值，模擬尺度為日尺度內每隔6 h，在中間正點時刻設置一個輻射點；在周期尺度內將每月的第15天設置為該月的一個輻射模擬日期。冠層光合作用生產力含義為1 h內虛擬冠層進行光合作用同化的CO2質量。模擬實驗結果如表1所示。

4.2 單葉光合特性參數選取

光合作用曲線方程中曲線凸度θ的值位于[0,1]之間，當θ為0時，光合作用曲線為一條直線，當θ為1時，光合作用曲線為一條直角曲線。根據文獻[22]中對不同葉齡針葉光響應曲線的描述，將θ設為0.3。內稟量子效率?i0是模型中表示植物葉片光合特性的參數。根據文獻[22]中對內稟量子效率的測定，將幼齡杉木葉片內稟量子效率設為0.03。CO2補償點為當植物葉片光合作用吸收CO2和呼吸作用釋放CO2相等，即凈光合速率為零時，外界CO2的濃度。根據文獻[23]中對不同種源杉木光合特性的比較中列出的CO2補償點的描述，將幼齡杉木的CO2補償點確定為150 μL/L。根據文獻[24]中相關描述，當溫度為25℃時，Rubisco酶的最大羧化能力Vm0為100 μmol/(m2·s)，Rubisco反應中，CO2米氏反應常數Kc為63 Pa，O2米氏反應常數Ko為36 000 Pa。

4.3 生物量估算結果分析

已知虛擬冠層光合作用速率p，根據生物量累積算法，該正點時刻的光合作用速率被視作前后各3 h的平均光合速率，生物累積量為

表1 部分輻射點時刻值、環境因子、冠層光合生產力、凈光合生產力模擬結果Tab.1 Simulation results of radiation, environmental factors, canopy photosynthetic productivity and net photosynthetic productivity at some radiation point

q=21 600kp

(15)

式中k——CO2量轉換為生物量的質量轉化系數，取30

由葡萄糖到碳水化合物的轉化系數取為0.72，生物量中的礦物質含量設為5%。按模擬周期內計算的各個輻射時刻，依次進行輻射模擬計算瞬時冠層光合速率，經系數轉換可得到周期內生物累積量。每月的生物累積量、生長周期內的生物累積量和模擬周期內總的生物累積量如表2所示。

表3中列出了不同種源杉木的光合生產力[23]。從表中可知，光合作用生產力同單株葉面積呈顯著正相關，各種種源的杉木光合生產力同單株葉面積之比位于[367，639]之間。其中福建沙縣種源杉木的比例為428.15。由表3結合輻射幼齡杉木模型的幾何結構參數可知，7—10月，各個輻射點上的虛擬冠層光合作用生產力同單株葉面積之比位于[397，454]之間。由于表3中數值測定的季節、光照條件、杉木冠層形狀、立地條件等因素均不統一，故無法作出精確的比較，但估算值接近于實測值，具有一定的可信度。其具體的精度，有待進一步的實測驗證。其次，表1模擬結果顯示了冠層光合作用速率隨季節變化，從夏季到冬季，冠層生產力呈遞減趨勢，與實際情況相符。

表2 幼齡杉木冠層生物量累積量估算Tab.2 Estimation of biomass accumulation for young Chinese fir g

表3 不同種源杉木光合生產力的比較Tab.3 Comparison of photosynthetic productivity for different provenances Chinese fir

5 模擬軟件程序

模擬結果的準確性取決于PAR強度、CO2濃度等光合因子的瞬時數值。該模型支持多種參數交互輸入，軟件界面如圖12所示。

圖12 模擬程序軟件界面Fig.12 Simulation program software interface

以計算虛擬杉木冠層瞬時光合速率為例，設定時間為2018年4月23日14時，地理位置為119°18′E，26°5′N，海拔為100 m。其他參數不變，改變CO2濃度，圖13中CO2濃度從3×10-4改為4.5×10-4。虛擬冠層光合作用生產力從2 188.39 mg/(h·株)提高為2 566.82 mg/(h·株)，符合植物光合作用特性。

6 結束語

光輻射(PAR)傳輸模型可以準確模擬某一位置、某一時刻冠層頂部太陽直射光照的入射方向；根據大氣透明系數，模擬冠頂的太陽直射、散射光強；結合CO2濃度、空氣濕度等因子，進行冠層尺度光合速率的計算。該模型可以準確模擬冠層內的光照輻射傳輸，在已知其他光合作用環境因子情況下，理論上給出一個準確冠層尺度光合作用模型。對5年生幼齡杉木進行三維建模，再進行光合作用環境因子模擬，在此基礎上完成太陽直射和天空散射在虛擬冠層內傳輸模擬與光合速率計算，在考慮呼吸作用的基礎上，完成一定周期內生物量累積量估算。估算結果在一定范圍內與實際相符。

圖13 模擬程序示例Fig.13 Simulator usage example