陳洪立，屠貴林，喬 欣，彭繼宇

(1.浙江理工大學(xué) 機械與自動控制學(xué)院，浙江 杭州 310018；2.浙江工業(yè)大學(xué) 特種裝備制造與先進加工技術(shù)教育部/浙江省重點實驗室，浙江 杭州 310023)

茶是我國最重要的飲品，茶葉的質(zhì)量直接影響其產(chǎn)品出口。目前，隨著消費者群體對身體健康的重視，其對產(chǎn)品質(zhì)量的要求越來越高，西方發(fā)達國家茶葉農(nóng)殘的進口標(biāo)準(zhǔn)日趨嚴(yán)格，茶葉的出口阻力加大。農(nóng)藥殘留直接影響到茶葉整體質(zhì)量安全[1-2]，因此茶葉農(nóng)藥殘留問題是當(dāng)下需要迫切解決的問題之一。目前，茶葉農(nóng)藥殘留檢測的方法普遍存在費時、費力和檢測精度不高等問題。通過檢索國內(nèi)外文獻，發(fā)現(xiàn)從分子角度解析茶葉農(nóng)藥降解規(guī)律的相關(guān)研究較少，合理的分子數(shù)值模擬方法結(jié)果可以作為實驗驗證的依據(jù)，彌補實驗工作的不足[3]。

近年來許多茶葉農(nóng)藥殘留的研究是基于偏微分方程的模型，較為典型的有指數(shù)負(fù)增長函數(shù)模型、多項式回歸模型灰色預(yù)測GM(1，1)模型等[4]。元胞自動機是一種離散的演化系統(tǒng)，它由許多小的組建構(gòu)成，具有模擬復(fù)雜系統(tǒng)的能力。于乃功等[5]在形態(tài)學(xué)結(jié)構(gòu)動力學(xué)模型的基礎(chǔ)上，建立了青霉素發(fā)酵的元胞自動機模型。因此，利用元胞自動機模擬農(nóng)藥殘留降解過程有著巨大優(yōu)勢。筆者采用仿真實驗與數(shù)值模擬結(jié)果比較的研究方法，探索元胞自動機算法模擬農(nóng)藥殘留降解過程的可行性。本研究以茶葉農(nóng)藥殘留的機理和推導(dǎo)出的動力學(xué)微分方程為基礎(chǔ)，建立農(nóng)藥殘留降解的元胞自動機模型，并根據(jù)設(shè)計的規(guī)則編程設(shè)計，進行仿真實驗，使得演化結(jié)果可視化。

1 農(nóng)殘降解的動力學(xué)微分方程

農(nóng)藥在植物上的衰減速度與農(nóng)藥的質(zhì)量分?jǐn)?shù)相關(guān)，質(zhì)量分?jǐn)?shù)越高，衰減速度越快，可用一級反動力學(xué)公式來表示，該公式是不考慮外在條件的理想化偏微分方程模型[4]，即

(1)

式中：y為t時刻的農(nóng)藥質(zhì)量分?jǐn)?shù)；t為施藥后時間；k為農(nóng)藥降解速率常數(shù)；a為農(nóng)藥的初始質(zhì)量分?jǐn)?shù)。解微分方程(1)可得

y=ae-kt

(2)

當(dāng)t取足夠小，則有

(3)

式中：t=kT0，T0為離散時間間隔。由式(2～3)推導(dǎo)出方程的離散形式，即

y(t+T0)-y(t)=T0ae-kt

(4)

式(4)為元胞自動機演化的設(shè)計基礎(chǔ)。

2 元胞自動機理論在農(nóng)藥降解中的應(yīng)用

2.1 農(nóng)藥殘留降解中的元胞

農(nóng)藥降解的元胞自動機模型中，在每個時刻kT，元胞都具有不同的狀態(tài)：存在或已降解。在該模型中，農(nóng)藥分子的殘留降解過程其實就是元胞自動機的各個元胞在不同狀態(tài)值之間的轉(zhuǎn)變過程[6-7]。演化規(guī)則作用于每個元胞節(jié)點，它是一個狀態(tài)轉(zhuǎn)移函數(shù)，下一時刻元胞的狀態(tài)僅與當(dāng)前自身狀態(tài)和鄰居狀態(tài)有關(guān)[8]，若當(dāng)前元胞存在，則掃描其周圍鄰居的狀態(tài)，如果周圍鄰居狀態(tài)滿足該元胞降解的條件，那么該元胞在下一時刻有一定概率降解。

2.2 元胞自動機模型建模的基本思想

農(nóng)藥殘留降解過程與施藥量、風(fēng)雨和時間等多種因素有關(guān)，其消解過程同放射性物質(zhì)衰變一樣，符合一級動力學(xué)方程。為了更好地模擬農(nóng)藥殘留降解過程，利用數(shù)值計算的方法對降解場景進行仿真[9-10]，對該模型作出以下假設(shè)：1) 每個元胞代表農(nóng)藥分子的集合；2) 根據(jù)建模需要，只考慮農(nóng)藥降解，不考慮中間產(chǎn)物的生成；3) 忽略外部因素的影響。在此假設(shè)條件下，分別對元胞、元胞空間、狀態(tài)集合、元胞鄰域和演化規(guī)則進行定義。最后利用計算機語言直觀地模擬該降解過程。

2.3 元胞自動機的應(yīng)用

根據(jù)農(nóng)藥殘留降解過程的動力學(xué)特性，建立了相關(guān)的元胞自動機模型。采用二維結(jié)構(gòu)作為農(nóng)藥降解的空間，演化規(guī)則根據(jù)農(nóng)藥降解機理以及一級動力學(xué)方程確定。模型中的每個元胞都代表特定數(shù)量的農(nóng)藥分子集合，它們具有不同的狀態(tài)，為了研究農(nóng)藥殘留降解速率與其質(zhì)量分?jǐn)?shù)的關(guān)系，采用鄰居較多的鄰域類型。

3 元胞自動機設(shè)計

本次模擬采用二維元胞自動機模型，模擬茶葉農(nóng)藥降解的二維元胞自動機可以定義為一個4 元組[7]，即

CA={t,Cells,CellSpace,Neiborhoods,Rules}

(5)

式中：t為離散時間；Cells為CA的基本元素，即元胞；Cellspace為元胞空間；Neighborhoods為CA的鄰域，即演化規(guī)則的定義域；Rules為CA的演化規(guī)則。

3.1 CA的離散時間

對農(nóng)藥降解的整個過程進行時間離散化，對于每一個離散時刻t=kT0，元胞自動機模型都會統(tǒng)計當(dāng)前農(nóng)藥的質(zhì)量分?jǐn)?shù)，隨后根據(jù)所有離散時間對應(yīng)的質(zhì)量分?jǐn)?shù)，繪制農(nóng)藥降解過程的農(nóng)殘曲線。

3.2 CA的元胞和元胞空間

元胞是元胞自動機的基本組成要素，每個元胞代表一群農(nóng)藥分子的集合，其結(jié)構(gòu)簡單，只有0和1兩個狀態(tài)。建立二維元胞空間，如圖1所示。元胞總數(shù)N為55×55，在農(nóng)藥降解的元胞自動機模型中，農(nóng)藥的質(zhì)量分?jǐn)?shù)=初始質(zhì)量分?jǐn)?shù)×(農(nóng)藥元胞個數(shù)/元胞總數(shù))。

圖1 二維元胞空間Fig.1 Two dimensional cell space

3.3 CA的元胞鄰域

在二維元胞自動機模型中，元胞自動機的鄰居類型主要有Moore型、擴展Moore型或者標(biāo)準(zhǔn)的Von.Neumann型等[11]，如圖2所示。

圖2 元胞自動機鄰域Fig.2 Neighbors of cellular automata

元胞鄰域是其演化規(guī)則的演化作用范圍，也就是元胞演化規(guī)則的定義域[12]。由于該模型需要研究農(nóng)藥分子與其質(zhì)量分?jǐn)?shù)之間的關(guān)系，因此這里采用鄰居較多的Moore型鄰域，即每個農(nóng)藥分子有8 個鄰居。

3.4 CA的元胞狀態(tài)

為了用數(shù)學(xué)模型再現(xiàn)農(nóng)藥分子降解動態(tài)，用數(shù)值計算的方法實現(xiàn)對整個復(fù)雜反應(yīng)過程的模擬[13-14]，并將最后結(jié)果顯示出來，該過程需要將不同時刻的元胞定義成不同的狀態(tài)。自動機中每一個元胞代表一群農(nóng)藥分子的集合，每個元胞具有0和1兩個不同狀態(tài)。如果自動機元胞Cell(i，j)在t時刻的狀態(tài)S(t)為0，則表示該元胞為空或已降解；如果為1，則表示該元胞仍存在。

3.5 CA的演化規(guī)則

自動機的演化規(guī)則具有局部特性，每個元胞的演化均依賴其相鄰的元胞狀態(tài)。也就是說，元胞的下一時刻狀態(tài)只與它鄰域的元胞狀態(tài)有關(guān)[15]。自動機按照演化規(guī)則進行并行演化。

根據(jù)茶葉農(nóng)藥分子降解機理和推導(dǎo)出來的微分方程，設(shè)計了自動機的演化規(guī)則為

如果Sij=1，那么

(6)

式中：N(i，j，t)為t時刻自動機元胞的鄰域內(nèi)細胞的個數(shù)；P(i，j，t)為t時刻狀態(tài)為1的元胞在下一刻演化為狀態(tài)0的演化概率；Pt為t時刻狀態(tài)為1的元胞在下一時刻演化為0的概率閾值，在仿真實驗中根據(jù)農(nóng)藥殘留作出必要調(diào)整；C為未降解農(nóng)藥和總體農(nóng)藥比值。

4 仿真實驗

以文獻[1]中的甲氰菊酯為研究對象，利用計算機語言，仿真編寫該元胞自動機的演化程序。程序功能包括：1) 計算每一時刻農(nóng)藥質(zhì)量分?jǐn)?shù)的值；2) 在二維空間上，展現(xiàn)該模型隨著時間而變化的演化圖案；3) 對每個離散時刻進行農(nóng)藥質(zhì)量分?jǐn)?shù)的統(tǒng)計，繪制其隨時間變化曲線。本研究采用的農(nóng)藥殘留降解數(shù)據(jù)引自文獻[1]中所提供的參數(shù)。模型參數(shù)分別采用表1數(shù)據(jù)。

表1 模型參數(shù)Table 1 Model parameter

表1中：初始質(zhì)量分?jǐn)?shù)a表示演化未開始時的農(nóng)藥質(zhì)量分?jǐn)?shù)；降解速率常速k表示甲氰菊酯降解的快慢；概率閾值Pt表示符合降解條件的元胞在下一時刻的降解概率，在仿真實驗中根據(jù)農(nóng)藥殘留作出必要調(diào)整；離散時間間隔T0即演化過程中將對每一個kT0時刻計算農(nóng)藥殘留質(zhì)量分?jǐn)?shù)。將這些參數(shù)賦值給元胞自動機的演化程序并執(zhí)行，分別得到農(nóng)藥微分方程質(zhì)量分?jǐn)?shù)和CA農(nóng)藥質(zhì)量分?jǐn)?shù)，結(jié)果如表2所示。

表2 不同時刻農(nóng)藥質(zhì)量分?jǐn)?shù)Table 2 Pesticide molecular concentration at different times

由表2可知：CA模型可以較好地擬合甲氰菊酯微分方程農(nóng)藥降解過程中的質(zhì)量分?jǐn)?shù)；在10～15 d兩者的質(zhì)量分?jǐn)?shù)相差較大，這說明由于前期農(nóng)藥質(zhì)量分?jǐn)?shù)下降后，模型的一些參數(shù)無法適應(yīng)，因此該模型模擬前期農(nóng)藥降解過程具有巨大優(yōu)勢，隨著農(nóng)藥質(zhì)量分?jǐn)?shù)下降到一定數(shù)值后，模型參數(shù)也需要作出必要調(diào)整。根據(jù)表2分別繪制農(nóng)藥微分方程質(zhì)量分?jǐn)?shù)和CA元胞農(nóng)藥質(zhì)量分?jǐn)?shù)曲線，結(jié)果如圖3所示。

圖3 茶葉中甲氰菊酯殘留量隨時間變化Fig.3 Changes of fenpropathrin residues in tea with time

圖3顯示：元胞自動機演化過程中，甲氰菊酯殘留量隨著時間的變化曲線較好地模擬了動力學(xué)微分方程模型確定的甲氰菊酯殘留量的變化曲線，同時也表現(xiàn)出了一定的隨機性?？梢钥吹皆诔跗陔A段CA元胞農(nóng)藥質(zhì)量分?jǐn)?shù)降解速度較快，考慮到可能是農(nóng)藥分子數(shù)量下降，農(nóng)藥降解的概率變大。隨著農(nóng)藥分子質(zhì)量分?jǐn)?shù)的下降，兩者的降解速度均緩慢下降。

如圖4所示，從元胞自動機演化過程不同時刻元胞狀態(tài)在二維空間上的演化圖案可以看出，演化圖案隨著時間的推移不斷變淺。剛開始農(nóng)藥殘留量達到最大值，隨著時間的推移，農(nóng)藥殘留量不斷下降。農(nóng)藥降解速度分為兩個階段：快速降解階段和緩慢降解階段。剛開始降解時，由于農(nóng)藥質(zhì)量分?jǐn)?shù)高，農(nóng)藥降解的速度較快，隨著農(nóng)藥質(zhì)量分?jǐn)?shù)的下降，農(nóng)藥降解的速度也開始減慢，可以推斷出農(nóng)藥降解速度與農(nóng)藥質(zhì)量分?jǐn)?shù)呈現(xiàn)正相關(guān)關(guān)系?？梢娫詣訖C演化圖案較好地演示了農(nóng)藥殘留降解的整體過程。

圖4 不同時刻的農(nóng)藥分子殘留狀態(tài)Fig.4 Pesticide residues at different times

5 結(jié) 論

在農(nóng)藥殘留降解機理及其動力學(xué)微分方程模型的基礎(chǔ)上，建立了模擬農(nóng)藥殘留降解的元胞自動機模型，并且對該模型進行了特性理論分析和仿真實驗。特性理論分析和仿真實驗結(jié)果均表明：該模型可以復(fù)現(xiàn)動力學(xué)微分方程所描述的農(nóng)藥殘留降解過程，二維圖像則直觀地展現(xiàn)了農(nóng)藥殘留降解的隨機性、不確定性等演化行為。本研究是假定單一外部環(huán)境變量的情況下，從一級動力學(xué)方程角度研究的農(nóng)藥殘留降解的元胞自動機模型。本次模擬討論點在于農(nóng)藥降解的整體過程，未對降解產(chǎn)物、中間體等物質(zhì)以及農(nóng)藥降解過程所發(fā)生的化學(xué)反應(yīng)進行研究。因此本課題下一個研究重點是：1) 利用元胞自動機模型，更具體深入地模擬農(nóng)藥分子殘留降解的各個過程，包括降解產(chǎn)物的生成，中間體以及外部因素對農(nóng)藥殘留降解過程的影響等；2) 為了從分子角度模擬農(nóng)藥殘留降解過程，以降解過程中的光解反應(yīng)和水解反應(yīng)為研究對象，以光解反應(yīng)和水解反應(yīng)方程定義演化規(guī)則，探究不同迭代步數(shù)情況下的農(nóng)藥降解形貌，分別研究光解反應(yīng)和水解反應(yīng)對農(nóng)藥殘留降解的影響等。