陳 剛 陳 斌 林郁欣 錢基德

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

0 引言

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

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

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

1 3D虛擬冠層構(gòu)建

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

1.1 枝干模擬

依據(jù)杉木枝干特征，利用L-系統(tǒng)規(guī)則構(gòu)建各個生長階段的幼齡杉木型枝干模型?；贚-系統(tǒng)的拓撲結(jié)構(gòu)模型構(gòu)建流程為：①依據(jù)杉木生長過程中拓撲變化規(guī)律，設(shè)計公理和產(chǎn)生式算法。②設(shè)計各分枝產(chǎn)生式符號、分枝角度、長度、半徑、產(chǎn)生概率等。③設(shè)計產(chǎn)生式約束條件。④程序?qū)σ?guī)則進行解析、字符迭代、繪制，生成三維樹木模型。

圖1 用于幼齡杉木虛擬冠層構(gòu)建的L-系統(tǒng)規(guī)則Fig.1 L-system rules for generating virtual canopy of young Chinese fir

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

1.2 葉片分布模擬

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

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

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

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

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

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

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

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

2.1.1太陽光直射方向計算

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

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

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

到達植物冠頂?shù)妮椛浒ㄌ栔鄙浜吞炜丈⑸鋬刹糠帧?/p>

2.1.2.1直射光強模擬

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

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

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

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

(1)

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

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

r0——日地平均距離

r——當前日地距離

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

垂直于太陽光線方向的太陽直射輻射強度與其穿過大氣層的路徑和大氣透明度有關(guān)。大氣量與太陽高度角和當?shù)氐牡匦胃叨纫约按髿鈮河嘘P(guān)。一定地形高度下的大氣量計算公式與海平面大氣量的計算公式為

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——大氣壓修正系數(shù)

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

p0——標準大氣壓

h——地形高度

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

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

(5)

式中τb——太陽直射大氣透明系數(shù)

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

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

(6)

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

2.1.2.2散射光強模擬

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

Is=E(r0/r)2τd

(7)

其中

τd=0.271-0.294τb

(8)

式中τd——散射輻射透明系數(shù)

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

2.1.2.3光合有效輻射PAR計算

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

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

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

2.2.1太陽直射PAR在虛擬冠層內(nèi)的分布模擬

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

圖6 太陽直射光在虛擬冠層內(nèi)的輻射傳輸模擬流程圖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在冠層內(nèi)分布模擬

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

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

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

(9)

方位角β計算公式為

(10)

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

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

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

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

3 虛擬冠層光合作用模型

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

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

(11)

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

P——最大光合速率

θ——曲線凸度

Pg——總光合速率

PAR——光合有限輻射強度

Pmax——最大光合速率

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

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

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

(12)

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

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

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

n——葉片數(shù)量

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

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

(13)

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

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

T——溫度

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

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

(14)

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

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

PO2——葉內(nèi)O2分壓

4 模擬實驗

4.1 光合同化物積累估算

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

4.2 單葉光合特性參數(shù)選取

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

4.3 生物量估算結(jié)果分析

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

表1 部分輻射點時刻值、環(huán)境因子、冠層光合生產(chǎn)力、凈光合生產(chǎn)力模擬結(jié)果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量轉(zhuǎn)換為生物量的質(zhì)量轉(zhuǎn)化系數(shù)，取30

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

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

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

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

5 模擬軟件程序

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

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

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

6 結(jié)束語

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

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