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

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

(1.中國林業(yè)科學(xué)研究院資源信息研究所，北京 100091；2.中南林業(yè)科技大學(xué)，湖南 長沙 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]。隨著對森林經(jīng)營認(rèn)識的深入，各國都意識到維護(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)營的熱點[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ù)動態(tài)地調(diào)整。權(quán)兵、王靈霞、張敏等人將林分生長方程和撫育間伐相結(jié)合，以林分密度為指標(biāo)研究連續(xù)動態(tài)的森林經(jīng)營，并設(shè)計了森林三維可視化系統(tǒng)[14-16]，但是僅采用林分密度指標(biāo)難以準(zhǔn)確分析出林分的綜合狀況。另外有學(xué)者對連年調(diào)查樣地研究發(fā)現(xiàn)可變生長率的生長方程在擬合單木生長方程方面有較大的優(yōu)勢[17-18]。因此，本研究基于模擬退火算法、可變生長率的生長方程并結(jié)合結(jié)構(gòu)化森林經(jīng)營思想，研究一種考慮林分結(jié)構(gòu)健康狀況和平均胸徑的多目標(biāo)森林經(jīng)營動態(tài)模擬方法，同時結(jié)合Unity3D 可視化引擎，實現(xiàn)林分經(jīng)營動態(tài)可視化模擬。

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

1.1 實驗區(qū)概況

實驗區(qū)設(shè)置在湖南省攸縣黃豐橋國有林場，林場地 處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%，是湖南省重點杉木大徑材推廣示范基地。

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

用常規(guī)的測樹學(xué)方法對選擇的樣地進(jìn)行每木檢尺和單木定位，測量其胸徑、樹高、冠高、冠幅、活枝下高、生長情況以及樹木的相對位置 (x, y,z)，按照1 年時間為間隔期，采集了2012—2017 年的連年調(diào)查數(shù)據(jù)，總結(jié)其中地位指數(shù)相同的5 塊樣地數(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 個結(jié)構(gòu)參數(shù)加以描述，分別代表了林分的分布狀況、林分的競爭情況、林分種間的隔離程度、林木的營養(yǎng)空間大小，根據(jù)前人的研究成果，健康的林分結(jié)構(gòu)在這4 個參數(shù)中均可以用林分均值定量化來表示[11]。本研究以杉木人工純林為研究對象，考慮到人工純林混交度為零，所以剔除混交度參數(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)行間伐木的判斷，按照惠剛盈等人提出的方法需要對每一株木的結(jié)構(gòu)參數(shù)進(jìn)行計算，并且需要對模擬出的多種預(yù)間伐后的方案進(jìn)行比較，雖然這種窮舉法可以得出經(jīng)營的最優(yōu)解，但是，一旦樹木的數(shù)量較大，計算起來可能會出現(xiàn)難以忍受的等待時間。考慮節(jié)省人力物力的條件下，滿足森林經(jīng)營整體需求的次優(yōu)解也是被允許的，因此本研究設(shè)計了一種基于模擬退火算法(SAA)的最優(yōu)間伐經(jīng)營方案選擇，模擬退火算法是求解組合問題的主要方法之一，它主要是基于Monte-Carlo 迭代求解策略的一種隨機(jī)尋優(yōu)算法，對于求解最優(yōu)組合快速、穩(wěn)定、有效[19-20]。基于SAA 的最優(yōu)間伐經(jīng)營方案選擇過程如下：

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

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

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

①初始化林分場景。通過讀取Excel 數(shù)據(jù)表，設(shè)杉木根部坐標(biāo)為 (x, y,z)，初始單木胸徑為D0，樹高為H0，冠幅為CW0，活枝下高為UBH0，則初始時林木的位置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)營方案選擇?；诳臻g結(jié)構(gòu)參數(shù)分析森林經(jīng)營緊迫性，根據(jù)分析結(jié)果和經(jīng)營者的實際考慮確定是否進(jìn)行間伐。一旦確定進(jìn)行間伐，首先通過競爭指數(shù)的比較初步篩選出預(yù)間伐木，建立間伐經(jīng)營方案鏈表，然后基于模擬退火算法確定最優(yōu)間伐經(jīng)營方案，被確定為間伐木對象其間伐屬性(Cut)設(shè)置為1，Tree_Position=(x,y,z),x,y,z設(shè)置為無限遠(yuǎn)，并用Unity3D 中Asset Store 下載的樹樁模型（如圖1）替換樹木原來的位置。

③林分生長動態(tài)變化?；赨nity3D 關(guān)鍵幀技術(shù)，將樹木的生長以一定的間隔化分為若干階段，每一階段用一個關(guān)鍵幀表示。在基于競爭指數(shù)的可變生長率的生長方程的模擬下，計算出第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)營。將攸縣黃豐橋國有林場5 號樣地數(shù)據(jù)導(dǎo)入，初始化渲染場景，分析樣地的空間結(jié)構(gòu)特征如表2，經(jīng)營緊迫性為適當(dāng)經(jīng)營，表明有其中一個指標(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 競爭指數(shù)的可變生長率的單木生長方程擬合與檢驗

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

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

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

表3 基于非線性擬合生長方程結(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 種樹木因子擬合值與實測值比較Fig.3 Comparison of fitting value and measured value of four tree factor

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

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

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

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

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

圖4 生長前（左）后（右）經(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)營動態(tài)過程Fig.6 Dynamic process of management

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

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

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

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

4 討論

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

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

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

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

5 結(jié)論

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