基于模擬退火算法的林分多目標(biāo)經(jīng)營動(dòng)態(tài)可視化模擬

沈 康，楊廷棟＊，張懷清，張 鴻，朱念福，劉 華

(1.中國林業(yè)科學(xué)研究院資源信息研究所，北京 100091；2.中南林業(yè)科技大學(xué)，湖南 長(zhǎng)沙 410004)

森林多目標(biāo)經(jīng)營優(yōu)化的研究由來已久，從傳統(tǒng)數(shù)學(xué)模型到人工智能算法[1-4]，主要側(cè)重于在經(jīng)營目標(biāo)和約束條件的指導(dǎo)下尋找一種合適的最優(yōu)解，多為森林某一階段的靜態(tài)經(jīng)營[5-7]。隨著對(duì)森林經(jīng)營認(rèn)識(shí)的深入，各國都意識(shí)到維護(hù)森林生態(tài)系統(tǒng)結(jié)構(gòu)與多樣性的重要性[8]，并相繼提出了保持森林生態(tài)結(jié)構(gòu)的經(jīng)營理念，如德國的近自然經(jīng)營和美國的生態(tài)系統(tǒng)管理，近年來以惠剛盈等人的結(jié)構(gòu)化森林經(jīng)營成為國內(nèi)外學(xué)者研究森林經(jīng)營的熱點(diǎn)[9-13]。以湯孟平、李建軍、曹旭鵬等人為代表，通過林分空間結(jié)構(gòu)分析構(gòu)建多目標(biāo)林空間優(yōu)化模型來進(jìn)行森林空間結(jié)構(gòu)調(diào)整[6-7,9]，但結(jié)構(gòu)調(diào)整不是連續(xù)動(dòng)態(tài)地調(diào)整。權(quán)兵、王靈霞、張敏等人將林分生長(zhǎng)方程和撫育間伐相結(jié)合，以林分密度為指標(biāo)研究連續(xù)動(dòng)態(tài)的森林經(jīng)營，并設(shè)計(jì)了森林三維可視化系統(tǒng)[14-16]，但是僅采用林分密度指標(biāo)難以準(zhǔn)確分析出林分的綜合狀況。另外有學(xué)者對(duì)連年調(diào)查樣地研究發(fā)現(xiàn)可變生長(zhǎng)率的生長(zhǎng)方程在擬合單木生長(zhǎng)方程方面有較大的優(yōu)勢(shì)[17-18]。因此，本研究基于模擬退火算法、可變生長(zhǎng)率的生長(zhǎng)方程并結(jié)合結(jié)構(gòu)化森林經(jīng)營思想，研究一種考慮林分結(jié)構(gòu)健康狀況和平均胸徑的多目標(biāo)森林經(jīng)營動(dòng)態(tài)模擬方法，同時(shí)結(jié)合Unity3D 可視化引擎，實(shí)現(xiàn)林分經(jīng)營動(dòng)態(tài)可視化模擬。

1 實(shí)驗(yàn)區(qū)概況及數(shù)據(jù)采集

1.1 實(shí)驗(yàn)區(qū)概況

實(shí)驗(yàn)區(qū)設(shè)置在湖南省攸縣黃豐橋國有林場(chǎng)，林場(chǎng)地 處113°04′～113°43′E、26°43′～27°06′N 之 間，最低海拔115 m，最高海拔1 270 m。隸屬于亞熱帶季風(fēng)氣候，平均氣溫17.8 ℃，年降水量1 410.8 mm，森林覆蓋率為90.07%，是湖南省重點(diǎn)杉木大徑材推廣示范基地。

1.2 數(shù)據(jù)采集方法

用常規(guī)的測(cè)樹學(xué)方法對(duì)選擇的樣地進(jìn)行每木檢尺和單木定位，測(cè)量其胸徑、樹高、冠高、冠幅、活枝下高、生長(zhǎng)情況以及樹木的相對(duì)位置 (x, y,z)，按照1 年時(shí)間為間隔期，采集了2012—2017 年的連年調(diào)查數(shù)據(jù)，總結(jié)其中地位指數(shù)相同的5 塊樣地?cái)?shù)據(jù)如表1 所示：

表1 樣地調(diào)查數(shù)據(jù)Table 1 Sample plot survey data

2 研究方法

2.1 空間結(jié)構(gòu)函數(shù)構(gòu)建

根據(jù)惠剛盈等人提出的結(jié)構(gòu)化森林的經(jīng)營思想，健康的林分結(jié)構(gòu)可以通過角尺度、大小比數(shù)、混交度、擁擠度4 個(gè)結(jié)構(gòu)參數(shù)加以描述，分別代表了林分的分布狀況、林分的競(jìng)爭(zhēng)情況、林分種間的隔離程度、林木的營養(yǎng)空間大小，根據(jù)前人的研究成果，健康的林分結(jié)構(gòu)在這4 個(gè)參數(shù)中均可以用林分均值定量化來表示[11]。本研究以杉木人工純林為研究對(duì)象，考慮到人工純林混交度為零，所以剔除混交度參數(shù)，建立空間結(jié)構(gòu)函數(shù)，公式如下，當(dāng)空間結(jié)構(gòu)函數(shù)Q(g)越小，林分結(jié)構(gòu)越接近健康水平。

2.2 基于模擬退火算法的最優(yōu)間伐經(jīng)營方案選擇

通過空間結(jié)構(gòu)函數(shù)來進(jìn)行間伐木的判斷，按照惠剛盈等人提出的方法需要對(duì)每一株木的結(jié)構(gòu)參數(shù)進(jìn)行計(jì)算，并且需要對(duì)模擬出的多種預(yù)間伐后的方案進(jìn)行比較，雖然這種窮舉法可以得出經(jīng)營的最優(yōu)解，但是，一旦樹木的數(shù)量較大，計(jì)算起來可能會(huì)出現(xiàn)難以忍受的等待時(shí)間。考慮節(jié)省人力物力的條件下，滿足森林經(jīng)營整體需求的次優(yōu)解也是被允許的，因此本研究設(shè)計(jì)了一種基于模擬退火算法(SAA)的最優(yōu)間伐經(jīng)營方案選擇，模擬退火算法是求解組合問題的主要方法之一，它主要是基于Monte-Carlo 迭代求解策略的一種隨機(jī)尋優(yōu)算法，對(duì)于求解最優(yōu)組合快速、穩(wěn)定、有效[19-20]。基于SAA 的最優(yōu)間伐經(jīng)營方案選擇過程如下：

初始化溫度T，設(shè)置溫度下降的速率r 和臨界值，計(jì)算出未間伐經(jīng)營條件下的空間結(jié)構(gòu)函數(shù)的值為初始解ans。根據(jù)競(jìng)爭(zhēng)指數(shù)的比較，建立間伐經(jīng)營方案鏈表，在某一溫度下隨機(jī)選取一個(gè)索引的間伐經(jīng)營方案的空間結(jié)構(gòu)函數(shù)作為新解tmp，當(dāng)Δ=tmp-ans＜0 時(shí)，接受tmp；當(dāng)Δ=tmp-ans＞0 時(shí)，根據(jù)數(shù)據(jù)量設(shè)置該溫度階段的循環(huán)次數(shù)K=10。當(dāng)循環(huán)終止時(shí)，按照exp((tmp-ans)/T)的概率接受次優(yōu)解tmp 并且將溫度下降為T=T*r，當(dāng)溫度降到臨界值時(shí)，得到間伐經(jīng)營方案的解趨向于全局最優(yōu)解。

2.3 經(jīng)營動(dòng)態(tài)可視化模擬

采用了張雄清等人[17]的基于Hegyi 競(jìng)爭(zhēng)指數(shù)的可變生長(zhǎng)率的單木生長(zhǎng)方程作為林分的生長(zhǎng)驅(qū)動(dòng)方程，結(jié)合基于模擬退火算法的最優(yōu)間伐經(jīng)營方案選擇實(shí)現(xiàn)林分的生長(zhǎng)與經(jīng)營的動(dòng)態(tài)交互，并在動(dòng)態(tài)的變化中滿足經(jīng)營者的需求。經(jīng)營動(dòng)態(tài)可視化模擬的過程如下：

①初始化林分場(chǎng)景。通過讀取Excel 數(shù)據(jù)表，設(shè)杉木根部坐標(biāo)為 (x, y,z)，初始單木胸徑為D0，樹高為H0，冠幅為CW0，活枝下高為UBH0，則初始時(shí)林木的位置Tree_Position= (x, y,z)，X-Z 平面縮放系數(shù)Tree_XZscale=D0/D(D 為原始模型的胸徑大小，如圖1 所示)，Y 方向上的縮放系數(shù)為Tree_Yscale=H0/H(H 為原始模型的樹高大小)。

圖1 初始杉木模型和樹樁模型Fig.1 Initial Chinese fir model &stump model

②林分經(jīng)營判斷與最優(yōu)經(jīng)營方案選擇。基于空間結(jié)構(gòu)參數(shù)分析森林經(jīng)營緊迫性，根據(jù)分析結(jié)果和經(jīng)營者的實(shí)際考慮確定是否進(jìn)行間伐。一旦確定進(jìn)行間伐，首先通過競(jìng)爭(zhēng)指數(shù)的比較初步篩選出預(yù)間伐木，建立間伐經(jīng)營方案鏈表，然后基于模擬退火算法確定最優(yōu)間伐經(jīng)營方案，被確定為間伐木對(duì)象其間伐屬性(Cut)設(shè)置為1，Tree_Position=(x,y,z),x,y,z設(shè)置為無限遠(yuǎn)，并用Unity3D 中Asset Store 下載的樹樁模型（如圖1）替換樹木原來的位置。

③林分生長(zhǎng)動(dòng)態(tài)變化。基于Unity3D 關(guān)鍵幀技術(shù)，將樹木的生長(zhǎng)以一定的間隔化分為若干階段，每一階段用一個(gè)關(guān)鍵幀表示。在基于競(jìng)爭(zhēng)指數(shù)的可變生長(zhǎng)率的生長(zhǎng)方程的模擬下，計(jì)算出第i 關(guān)鍵幀的杉木胸徑設(shè)為Di，樹高設(shè)為Hi，冠幅設(shè)為CWi，活枝下高設(shè)為UBHi，則此刻林木的位置Tree_Position= (x, y,z)，X-Z 平面縮放系數(shù)Tree_XZscale=Di/Di-1(i＞1)，Y 方向上的縮放系數(shù)為Tree_Yscale=Hi/Hi-1(i＞1)。

④以健康林分結(jié)構(gòu)特征和目標(biāo)林分胸徑為經(jīng)營目標(biāo)作為停止條件，將步驟②和步驟③進(jìn)行循環(huán)往復(fù)，直到滿足經(jīng)營目標(biāo)為止。

3 結(jié)果與分析

3.1 林分經(jīng)營緊迫性性分析

根據(jù)惠剛盈等人提出的健康林分的判斷標(biāo)準(zhǔn)[11]，將純林林分是否間伐劃分為4 種情況：不經(jīng)營、適當(dāng)經(jīng)營、需要經(jīng)營、必須經(jīng)營。將攸縣黃豐橋國有林場(chǎng)5 號(hào)樣地?cái)?shù)據(jù)導(dǎo)入，初始化渲染場(chǎng)景，分析樣地的空間結(jié)構(gòu)特征如表2，經(jīng)營緊迫性為適當(dāng)經(jīng)營，表明有其中一個(gè)指標(biāo)不滿足條件，由表2可以看出，平均角尺度不在[0.475～0.517]的健康林分區(qū)間。

表2 初始林分結(jié)構(gòu)參數(shù)分析結(jié)果Table 2 Analysis results of initial stand structure parameters

3.2 基于Hegyi 競(jìng)爭(zhēng)指數(shù)的可變生長(zhǎng)率的單木生長(zhǎng)方程擬合與檢驗(yàn)

3.2.1 基于Hegyi 競(jìng)爭(zhēng)指數(shù)的可變生長(zhǎng)率的單木生長(zhǎng)方程擬合 通過改進(jìn)張雄清等人的基于Hegyi 競(jìng)爭(zhēng)指數(shù)的可變生長(zhǎng)率的單木生長(zhǎng)方程，進(jìn)行胸徑生長(zhǎng)擬合，參數(shù)方程如公式（2）：

D(i,t+1)： 表示第i 棵樹在t+1 年的胸徑；D ( i,t)：表示第i 棵樹在t 年的胸徑；At：表示林分的年齡；CIi：表示第i 棵樹的Hegyi-Ci 的競(jìng)爭(zhēng)指數(shù)；B1、B2、B3 為參數(shù)。

對(duì)前4 塊樣地的連年調(diào)查數(shù)據(jù)進(jìn)行處理，剔除波動(dòng)較大的異常值數(shù)據(jù)，然后通過ForStat2.2 進(jìn)行非線性擬合，擬合胸徑的生長(zhǎng)方程結(jié)果如表3 所示，參數(shù)B1、B2、B3 取值的95%的置信區(qū)間上下限均沒有跨越0，表明參數(shù)值可用；表4 數(shù)據(jù)顯示B1、B2、B3 三個(gè)參數(shù)的相關(guān)性較低，說明函數(shù)構(gòu)造合理，故擬合出的胸徑生長(zhǎng)的方程。在得出胸徑的生長(zhǎng)方程后，通過線性回歸分析可以得出胸徑與樹高、冠幅以及樹高與枝下高之間存在明顯的線性關(guān)系，如圖2 所示，從圖中可以看出D-H、DCW、H-UBH 的R2分別是0.790 8、0.541 8、0.578 6，表明胸徑樹高相關(guān)性較高，胸徑冠幅和樹高枝下高相關(guān)性較差一些。

表3 基于非線性擬合生長(zhǎng)方程結(jié)果Table 3 Fitting growth equation results based on nonlinear method

表4 參數(shù)的漸近相關(guān)系數(shù)Table 4 Asymptotic correlation coefficient of parameters

圖2 3 種樹木因子線形回歸結(jié)果（胸徑-樹高、胸徑-冠幅、樹高-枝下高）Fig.2 Linear regression results of three tree factors (DBH-HDBH-CW, H-UBH)

圖3 4 種樹木因子擬合值與實(shí)測(cè)值比較Fig.3 Comparison of fitting value and measured value of four tree factor

3.2.2 數(shù)據(jù)檢驗(yàn) 選擇5 號(hào)樣地的連年數(shù)據(jù)對(duì)上述單木生長(zhǎng)方程結(jié)果進(jìn)行驗(yàn)證，由圖3 可知胸徑與樹高估計(jì)值與實(shí)測(cè)值偏差較小，枝下高與冠幅的估計(jì)值與實(shí)測(cè)值偏差較大，這一趨勢(shì)與上面擬合的單木生長(zhǎng)方程R2的大小比較結(jié)果一致，雖然枝下高與冠幅的擬合結(jié)果與真實(shí)值有一定的差距，但總體上具有一定的相關(guān)性，因此基本滿足生長(zhǎng)擬合的要求。

3.3 基于模擬退火算法的最優(yōu)間伐經(jīng)營方案選擇結(jié)果

根據(jù)間伐量不得大于生長(zhǎng)量的原則，在競(jìng)爭(zhēng)指數(shù)的比較中對(duì)采伐木進(jìn)行初步的篩選，本研究根據(jù)樣地實(shí)際情況，篩選出競(jìng)爭(zhēng)指數(shù)較大的前10%作為預(yù)間伐木，然后對(duì)初步篩選后的間伐木進(jìn)行組合排列，制定出所有經(jīng)營方案的集合，通過模擬退火算法對(duì)所有間伐經(jīng)營方案進(jìn)行最優(yōu)解的計(jì)算，得出最優(yōu)采伐方案，同時(shí)對(duì)最優(yōu)解與未間伐下的空間函數(shù)進(jìn)行比較，以確定是否需要實(shí)施經(jīng)營，由于初始林分（16 a）時(shí)林分經(jīng)營緊迫性分析為適當(dāng)經(jīng)營（如圖4 左），所以對(duì)該樣地進(jìn)行連續(xù)5 a 的生長(zhǎng)模擬，再對(duì)生長(zhǎng)后（21 a）的林分重新進(jìn)行經(jīng)營緊迫性分析，分析結(jié)果為需要經(jīng)營（如圖4 右），因此對(duì)其進(jìn)行經(jīng)營模擬。圖5 左為通過競(jìng)爭(zhēng)指數(shù)初步篩選的間伐木結(jié)果，從圖中可以看出初步篩選了23株預(yù)采伐木，圖5 右為基于模擬退火算法判斷的最優(yōu)間伐經(jīng)營方案選擇結(jié)果，實(shí)際間伐了4 株樹木。

3.4 經(jīng)營動(dòng)態(tài)可視化模擬過程

經(jīng)營動(dòng)態(tài)可視化模擬過程可以抽象成林分狀態(tài)場(chǎng)景、林分結(jié)構(gòu)與經(jīng)營緊迫性分析、林分生長(zhǎng)、間伐經(jīng)營方案和經(jīng)營目標(biāo)5 個(gè)模塊，它們之間相互連接與迭代構(gòu)成了林分動(dòng)態(tài)生長(zhǎng)與經(jīng)營過程。如圖6所示，林分狀態(tài)場(chǎng)景用于模擬林分不同時(shí)期的三維渲染；林分結(jié)構(gòu)與經(jīng)營緊迫性分析用于分析不同時(shí)期林分結(jié)構(gòu)的狀況和經(jīng)營緊迫性的等級(jí)，主要側(cè)重于數(shù)字化的表達(dá)，從而判斷林分是否生長(zhǎng)還是間伐；林分生長(zhǎng)通過生長(zhǎng)率的變化模擬出下一年的生長(zhǎng)量；間伐經(jīng)營方案則是在分析結(jié)構(gòu)與胸徑的基礎(chǔ)上以間伐樹木控制林分的狀態(tài)；在實(shí)施生長(zhǎng)或者間伐后，判斷經(jīng)營目標(biāo)，重新回到林分狀態(tài)，如果滿足經(jīng)營要求停止活動(dòng)，否則繼續(xù)以上步驟。

圖4 生長(zhǎng)前（左）后（右）經(jīng)營緊迫性分析Fig.4 Analysis of management urgency before (left) and after (right) growth

圖5 間伐木判斷（左）與間伐執(zhí)行結(jié)果（右）Fig.5 Judgment of thinning (left)and results of thinning execution(right)

圖6 經(jīng)營動(dòng)態(tài)過程Fig.6 Dynamic process of management

本研究以林分平均胸徑達(dá)到30 cm 和健康的林分結(jié)構(gòu)為經(jīng)營目標(biāo)，對(duì)5 號(hào)樣地進(jìn)行經(jīng)營動(dòng)態(tài)可視化模擬，模擬結(jié)果如圖7 所示。5 號(hào)樣地初始林分場(chǎng)景時(shí)（如表5 所示），林分平均年齡為16 a，平均胸徑為20.21 cm，平均角尺度為0.64，即分布狀態(tài)為團(tuán)狀分布，經(jīng)營緊迫性為適當(dāng)經(jīng)營，所以對(duì)5 號(hào)樣地進(jìn)行采取生長(zhǎng)模擬。本研究將間伐間隔期設(shè)定為5 a，模擬出5 號(hào)樣地21 年生的林分狀況，此時(shí)林分平均胸徑為23.79 cm，平均角尺度沒有變化，因?yàn)闆]有對(duì)林分進(jìn)行間伐模擬，林分樹木相對(duì)位置沒有發(fā)生變化，林分擁擠度為0.81，表明林分平均冠幅變大，林分逐漸擁擠，經(jīng)營緊迫性為需要經(jīng)營，對(duì)林分進(jìn)行間伐模擬，模擬結(jié)果為圖7 中上所示，間伐了其中的4 株樹木，間伐后重新對(duì)林分進(jìn)行分析，結(jié)果表明林分經(jīng)營緊迫性仍為需要經(jīng)營，但間伐后如圖8 所示密集度分布向右偏移，表明擁擠程度有所下降，大小比數(shù)向左偏移，表明垂直競(jìng)爭(zhēng)下降，所以經(jīng)營仍然有效。由于平均胸徑未達(dá)到目標(biāo)，因此繼續(xù)對(duì)林分進(jìn)行生長(zhǎng)模擬。在連續(xù)的生長(zhǎng)模擬與間伐模擬后，當(dāng)林分平均年齡模擬到31 a 時(shí)，林分平均胸徑達(dá)到30.02 cm，林分經(jīng)營緊迫性為必須經(jīng)營，表明仍需要對(duì)林分進(jìn)行經(jīng)營模擬，對(duì)林分進(jìn)行間伐模擬后，經(jīng)營緊迫性等級(jí)下降一級(jí)，間伐后平均胸徑為30.10 cm，故結(jié)束此次林分模擬。

圖7 經(jīng)營動(dòng)態(tài)可視化模擬過程Fig.7 Management dynamic visualization simulation process

表5 經(jīng)營統(tǒng)計(jì)表Table 5 Management statistics

圖8 經(jīng)營前后結(jié)構(gòu)參數(shù)等級(jí)差值分布(左21a，右26a)Fig.8 Distribution of structural parameter grades before and after management(left 21a, right26a)

4 討論

本研究以湖南省攸縣黃豐橋國有林場(chǎng)5 塊長(zhǎng)期固定樣地為數(shù)據(jù)源，通過林分結(jié)構(gòu)參數(shù)，建立了空間結(jié)構(gòu)函數(shù)，并采用模擬退火算法進(jìn)行間伐經(jīng)營方案的最優(yōu)選擇，用以指導(dǎo)經(jīng)營采伐；同時(shí)，研建了基于競(jìng)爭(zhēng)指數(shù)的可變生長(zhǎng)率的生長(zhǎng)方程，得到了胸徑與樹高、胸徑與冠幅、樹高與枝下高之間具有相關(guān)性性，且胸徑與樹高之間相關(guān)性最強(qiáng)；并基于Unity3D 三維渲染引擎，模擬了林分動(dòng)態(tài)經(jīng)營采伐過程，效果形象逼真，運(yùn)行流暢，但也存在一些不足之處。

1）本研究模擬的為同齡純林，未涉及混交異齡的經(jīng)營，因此空間參數(shù)選擇較少，對(duì)于混交度、開敞度等其他空間結(jié)構(gòu)參數(shù)等沒有引入，因此經(jīng)營緊迫性劃分等級(jí)較為粗放，另外沒有考慮補(bǔ)植經(jīng)營措施的影響，因此需要在以后的研究中進(jìn)行深入。

2）基于模擬退火算法的最優(yōu)間伐經(jīng)營方案選擇與模擬退火溫度的上下限以及下降速率有關(guān)，本研究考慮了計(jì)算機(jī)的性能，適當(dāng)?shù)慕档土藴囟鹊纳舷孪藓吞岣吡讼陆邓俾省Ａ硗猓紤]到樣地林木數(shù)量較少，故在競(jìng)爭(zhēng)指數(shù)初步篩選間伐木時(shí)，僅選擇了競(jìng)爭(zhēng)指數(shù)較大的前10%的樹木作為預(yù)間伐木，因此可能會(huì)出現(xiàn)間伐后的效果不是太明顯，間伐結(jié)果不一定是最優(yōu)經(jīng)營解的情況。

3）對(duì)于可視化模擬大面積動(dòng)態(tài)經(jīng)營的能力還有所欠缺，主要包括兩個(gè)方面的原因：一方面大面積樹木的渲染需要消耗大量GPU 性能，雖然可以通過層次細(xì)節(jié)（LOD）和射線碰撞檢測(cè)進(jìn)行緩解，但是問題依舊存在；另一方面前期的空間結(jié)構(gòu)參數(shù)的計(jì)算以及組合經(jīng)營方案集的遍歷也會(huì)隨著樹木數(shù)量的上升，出現(xiàn)等待的情況，在以后的研究中仍需要研究有效的算法來進(jìn)行組織。

5 結(jié)論

模擬退火算法是解決大量組合排列問題的有效方法之一，通過合理的構(gòu)建目標(biāo)函數(shù)，對(duì)于解決森林經(jīng)營中存在的大量組合排列問題也同樣具有很好的效果。本研究通過可視化的手段形象地描述了模擬退火算法在森林經(jīng)營中應(yīng)用的可行性，生動(dòng)地展示了在多目標(biāo)需求下的森林生長(zhǎng)與經(jīng)營的整個(gè)過程，對(duì)于進(jìn)一步提高經(jīng)營水平和經(jīng)營的精細(xì)化程度具有一定的作用。