松材線蟲病害擴散預測模型的設計

姚 緣

上海市嘉定區林業站，上海 201822

借助計算機技術，林業工作人員嘗試構建松材線蟲病害擴散預測模型。一方面，根據松材線蟲病害聚集分布特點，建立基于區域的病蟲害擴散預測模型。與傳統的數值測算方法相比，基于區域的預測模型對病蟲害擴散的推演全部在二維平面上完成，通過對樣條曲線進行擬合，得到病蟲害區域參數，并結合能量函數預估一定時間內的病蟲害擴散區域?；谏植∠x害影響力，使用回歸方法預測病蟲害在二維平面上的傳播方向與距離，經過參數校正之后得到病蟲害擴散預測輪廓。與元胞自動機方法相比，基于區域病蟲害擴散預測模式能夠更為準確地預估病蟲害擴散機制，能夠以較為直觀的方式幫助林業從業者了解松材線蟲病害擴散情況。

另一方面，根據松材線蟲病害高傳染率的特點，創建基于傳染病動力學的病蟲害擴散預測模型。該模型通過對比自然狀態與施藥狀態下，松材線蟲病害的擴散差異，以動態化的方式展示SEI模型（自然狀態下松材線蟲病害擴散模型）與SEIR（施藥狀態下松材線蟲病害擴散模型）模擬擴散情況，為松材線蟲病害防治提供數據參考。以下針對這2種病害擴散預測模型進行詳細分析。

1 基于區域的病蟲害擴散預測模型

1.1 圖像擬合

構建預測模型之后，相關工作人員需要使用專業的信息提取設備獲得目標區域的矢量數據，通過分批處理病蟲害擴散區域圖像對矢量圖像進行擬合，進而得到病蟲害區域坐標點集合。該過程主要分為3個步驟：（1）調節圖像尺寸，根據具體的設計需求，將此次模型設計中的圖像尺寸調節為W=500，H=500；（2）將矢量圖像轉變為二值圖像，其中深色的區域代表未受到松材線蟲病害影響，淺色區域代表已經受到了松材線蟲病害影響；（3）提取二值圖像上的病蟲害擴散區域輪廓點坐標[1]。

提取輪廓點坐標之后，采用均勻間隔提取模式，選取輪廓指標點并將其擬合為曲線型值點，并根據實際需要，調整設計型值點的精度。需要注意的是，閉合曲線首尾需要重合，因此，設計人員假定病蟲害輪廓提取點的數量為Ncdp-1個(此處的Ndtp-1數值為100)，則共需要提取Ndtp個型值點。為了保障采樣點在圖像呈均勻分布狀態，工作人員要確定每一個型值點在二維平面上的間隔，其公式如下：

上述2個公式中，floor表示向下取整函數，Ncdp代表坐標點的數量，Nitv則表示型值點取值間隔，Nitv則表示雙精度浮點型。先利用浮點數確定二維圖像中的坐標值，在確定坐標序數的情況下對大型值數據進行四舍五入，進而得到間隔數據（圖1）。

圖1 Nitv與Nitv提取型值點的區別

分析圖1可以發現，2種數值提取方式在輪廓點末端存在巨大的差異。第2種數值提取方式中，目標區域（Ω區）提取的型值點比第1種方式提取的型值點更加均勻。因此，此次模型設計決定使用第22種型值提取方式。

1.2 病蟲害擴散模擬

1.2.1 確定病蟲害擴散因素 在正式開始模擬之前，需要確定影響松材線蟲病害擴散的主要因素。經過研究發現，一方面，病蟲害的擴散與林地的坡度呈正相關關系，即坡度越大，病蟲害發病概率越高；另一方面，林地的海拔與病蟲害擴散呈負相關關系，即林地的海拔越高，松材線蟲病害的發病概率越低?；谶@種特點，相關工作人員將上述2種因素作為影響松材線蟲病害擴散的特征。

1.2.2 計算病蟲害無差別擴散輪廓 假設某林區發生松材線蟲病害，設計人員將不受海拔與坡度影響的病蟲害擴散輪廓與受海拔與坡度影響的病蟲害擴散輪廓進行對比。在對比過程中，相關工作人員將只由于松材線蟲繁殖而產生的病蟲害擴散稱為無差別擴散。在計算無差別擴散時，工作人員利用3次B樣條閉合曲線面積的擬合，得到無差別擴散輪廓，利用面積值計算公式得到無差別擴散面積與質心數據，其公式如下：

式（3）中，m00代表封閉3次B樣條曲線零階距，即擴散面積。xi與yi為B樣條曲線控制點的橫縱坐標。在得到面積之后，利用1階距m10，m01計算質心，其公式為：

式（4）中，xc表示型心的橫坐標，yc表示型心的縱坐標。

在確定病蟲害擴散輪廓面積和質心之后，可以開始對無差別擴散情況進行模擬。先計算2次病蟲害擴散范圍的面積倍數，假設St為在t時刻病蟲害的面積，則t時刻與t+1時刻病蟲害面積的倍數為通過計算得到t+1時刻病蟲害無差別擴散輪廓Pnui(t+1)，其計算公式為：

式（5）、（6）中，Pnui(t+1)（0）與Pnui(t+1)（1）表示在t+1時擴散輪廓的橫坐標與縱坐標，而Pt（0）與Pt（1）表示在t時病蟲害實際擴散的橫坐標與縱坐標[2]。

1.3 預測模擬計算

假設t+1時刻無差別擴散輪廓的型值 點 為Pnui(t+1):(x0,x1,x2,…,xNdtp－1)，考 慮坡度和海拔因素的病蟲害擴散輪廓型值點為Pt+1:(x0,x1,x2,…,xNdtp－1)，對比二者之間的差異。

在確定Emin數值的基礎上，病蟲害擴散輪扣型值點之間最小距離Dist與對應點角度Angle可視為模型數據集中的Y值，計算對應點連接后形成的直線與水平方向射線之間的夾角，計算公式如下：

在得到θ的具體數據之后，工作人員就能夠預測t+1時刻的病蟲害擴散輪廓（表1）。

表1 松材線蟲病害擴散預測流程

2 基于傳染病動力學松材線蟲病蟲害擴散的預測模型

構建基于區域的松材線蟲病蟲害擴散預測模型，能夠使相關工作人員以較為直觀的方式觀測病蟲害擴散情況，但該模型無法準確描述病蟲害擴散流行規律。針對這一問題，工作人員嘗試利用傳染病動力學理論構建SEI模型（自然狀態下松材線蟲病害擴散模型）和SEIR模型（施藥狀態下松材線蟲病害擴散模型）。通過對比，幫助工作人員了解松材線蟲病害擴散動態化過程。

2.1 建立模型

作為一種外來入侵型病原生物，松材線蟲對林業的影響巨大，這種生物將松褐天牛作為主要宿主，其繁殖過程會堵塞松樹內部的導管，導致松樹枯萎。在正式構建模型之前，工作人員假定松褐天牛在林區之中均勻分布，并且松材線蟲僅在病樹中存在。假設松材線蟲的傳播方式與人類傳染病的傳播方式相同，分別設計SEI與SEIR模型[3]。

2.1.1 SEI模型 如果林業管理者并未發現林區中存在松材線蟲，則在自然狀態下，只需要五年左右的時間，松材線蟲就能摧毀一整片松林。由于松材線蟲為外來入侵物種，因此，中國本土松樹并不具備針對松材線蟲的抵抗能力。一旦松樹被感染，如果不進行人工干預，松樹幾乎沒有可能自行恢復健康。基于這種情況，在設計模型時，不考慮非人工干預狀態下病樹轉變為移出者R的可能。

當松樹枯萎后，松材線蟲還可以繼續在松樹內存活，以松樹內的真菌為食。同時，松材線蟲的主要傳播媒介松褐天牛習慣在枯萎的松樹中產卵繁殖，因此，在設計SEI模型時仍然將枯萎的松樹視為感染者I。

此外，松材線蟲病有著很長的潛伏期，松樹在發病初期其表面并無明顯變化。通常情況下，松褐天牛并不會從健康的松樹中飛出，因此可以將處于染病初期的松樹視作潛伏者E。

基于上述背景，工作人員將自然狀態下的松材線蟲病害擴散模型劃分為3個艙室，再利用SEI模型在不同的艙室中進行模擬（圖2）。

圖2 SEI模型

SEI模型中，S表示健康但不具備免疫力的松樹；E表示處于染病初期但不具備傳染性的松樹；I表示已經受到感染的松樹，其微分方式模型如下：

該公式中，系數β代表染病松樹將松材線蟲傳染給其他松樹的有效擴散率，其中，潛伏者E有一定的概率（μ）表現出病癥，松材線蟲病害的潛伏期即1/μ。

2.1.2 SEIR模型 目前，松材線蟲防治方式主要包括樹干注藥、砍伐疫木以及誘捕器誘殺等方式，橫向對比之下，樹干注藥的防治效果較為優秀。因此，此次設計以樹干注藥為基礎建立SEIR模型。

在松樹的樹干中注入藥液，利用松樹的蒸騰作用將藥液輸送至松樹的各個部位。需要注意的是，由于健康的松樹內部含有松脂，不具備脂溶性特點的藥液無法融入松脂，使得藥液無法被輸送至松樹的每一個部位[4]。因此，健康的松樹有一定的概率成為移出者R。

一些松樹年代久遠，具有很高的觀賞價值和研究價值，因此相關工作人員研制出用于急救的藥液，但此類藥品的價格較高，無法普遍使用。因此，在建立SEIR模型時，排除松樹因得到了急救藥而轉好的可能性，其微分方程模型如下：

該公式中，系數σ表示健康松樹在注射藥液之后獲得免疫能力并轉變為移出者R的概率。

2.2 可視化設計

為了能夠讓相關工作人員以較為直觀的方式觀察松材線蟲病擴散過程，設計人員對上述模型進行可視化升級。利用計算機技術設計SEI與SEIR模型可視化界面（圖3）。

圖3 SEI與SEIR模型可視化界面

分析圖3可以發現，此次模型設計使用的參數中，C（林區松樹總數量）為29 753，點擊“All Green”按鈕，UI頁面實時反饋的深色部分為健康松樹分布位置，點擊“SEIimulating”和“SEIRimulating”按鈕，可以發現UI界面中出現黃色區域，代表受到松材線蟲病害影響的松樹分布范圍。點擊“Intitalization”按鈕使得UI頁面恢復初始狀態[5]。

此時，工作人員可以根據實際需求設定迭代間隔時間，并點擊橙色的“開始按鈕”，UI頁面會按照設定好的迭代時間，展示健康松樹與受到松材線蟲病害影響的松樹分布變化情況，通過這種方式幫助工作人員了解該林業松材線蟲病害擴散情況。

3 結束語

松材線蟲病害對中國林業資源開發工作產生了很大的影響，如何解決松材線蟲病害問題已經成為擺在林業從業者面前的一個難題。想要提高病蟲害防治效率，就要對松材線蟲病害擴散情況進行精準預測，并以此為基礎，制定切實可行的預防措施。因此，相關工作人員嘗試基于病害區域和傳染病動力學建立松材線蟲病害擴散預測模型，并對該模型進行可視化設計，方便林業從業人員對松材線蟲病害的擴散進行動態觀察，為更好地保護林業資源提供技術支撐。