基于模擬退火算法的林分多目標經營動態可視化模擬

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

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

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

1 實驗區概況及數據采集

1.1 實驗區概況

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

1.2 數據采集方法

用常規的測樹學方法對選擇的樣地進行每木檢尺和單木定位，測量其胸徑、樹高、冠高、冠幅、活枝下高、生長情況以及樹木的相對位置 (x, y,z)，按照1 年時間為間隔期，采集了2012—2017 年的連年調查數據，總結其中地位指數相同的5 塊樣地數據如表1 所示：

表1 樣地調查數據Table 1 Sample plot survey data

2 研究方法

2.1 空間結構函數構建

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

2.2 基于模擬退火算法的最優間伐經營方案選擇

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

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

2.3 經營動態可視化模擬

采用了張雄清等人[17]的基于Hegyi 競爭指數的可變生長率的單木生長方程作為林分的生長驅動方程，結合基于模擬退火算法的最優間伐經營方案選擇實現林分的生長與經營的動態交互，并在動態的變化中滿足經營者的需求。經營動態可視化模擬的過程如下：

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

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

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

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

④以健康林分結構特征和目標林分胸徑為經營目標作為停止條件，將步驟②和步驟③進行循環往復，直到滿足經營目標為止。

3 結果與分析

3.1 林分經營緊迫性性分析

根據惠剛盈等人提出的健康林分的判斷標準[11]，將純林林分是否間伐劃分為4 種情況：不經營、適當經營、需要經營、必須經營。將攸縣黃豐橋國有林場5 號樣地數據導入，初始化渲染場景，分析樣地的空間結構特征如表2，經營緊迫性為適當經營，表明有其中一個指標不滿足條件，由表2可以看出，平均角尺度不在[0.475～0.517]的健康林分區間。

表2 初始林分結構參數分析結果Table 2 Analysis results of initial stand structure parameters

3.2 基于Hegyi 競爭指數的可變生長率的單木生長方程擬合與檢驗

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

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

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

表3 基于非線性擬合生長方程結果Table 3 Fitting growth equation results based on nonlinear method

表4 參數的漸近相關系數Table 4 Asymptotic correlation coefficient of parameters

圖2 3 種樹木因子線形回歸結果（胸徑-樹高、胸徑-冠幅、樹高-枝下高）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 數據檢驗 選擇5 號樣地的連年數據對上述單木生長方程結果進行驗證，由圖3 可知胸徑與樹高估計值與實測值偏差較小，枝下高與冠幅的估計值與實測值偏差較大，這一趨勢與上面擬合的單木生長方程R2的大小比較結果一致，雖然枝下高與冠幅的擬合結果與真實值有一定的差距，但總體上具有一定的相關性，因此基本滿足生長擬合的要求。

3.3 基于模擬退火算法的最優間伐經營方案選擇結果

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

3.4 經營動態可視化模擬過程

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

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

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

圖6 經營動態過程Fig.6 Dynamic process of management

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

圖7 經營動態可視化模擬過程Fig.7 Management dynamic visualization simulation process

表5 經營統計表Table 5 Management statistics

圖8 經營前后結構參數等級差值分布(左21a，右26a)Fig.8 Distribution of structural parameter grades before and after management(left 21a, right26a)

4 討論

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

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

2）基于模擬退火算法的最優間伐經營方案選擇與模擬退火溫度的上下限以及下降速率有關，本研究考慮了計算機的性能，適當的降低了溫度的上下限和提高了下降速率。另外，考慮到樣地林木數量較少，故在競爭指數初步篩選間伐木時，僅選擇了競爭指數較大的前10%的樹木作為預間伐木，因此可能會出現間伐后的效果不是太明顯，間伐結果不一定是最優經營解的情況。

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

5 結論

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