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

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

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

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

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

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

(1)

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

y=ae-kt

(2)

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

(3)

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

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

(4)

式(4)為元胞自動(dòng)機(jī)演化的設(shè)計(jì)基礎(chǔ)。

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

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

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

2.2 元胞自動(dòng)機(jī)模型建模的基本思想

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

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

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

3 元胞自動(dòng)機(jī)設(shè)計(jì)

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

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

(5)

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

3.1 CA的離散時(shí)間

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

3.2 CA的元胞和元胞空間

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

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

3.3 CA的元胞鄰域

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

圖2 元胞自動(dòng)機(jī)鄰域Fig.2 Neighbors of cellular automata

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

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

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

3.5 CA的演化規(guī)則

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

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

如果Sij=1，那么

(6)

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

4 仿真實(shí)驗(yàn)

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

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

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

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

由表2可知：CA模型可以較好地?cái)M合甲氰菊酯微分方程農(nóng)藥降解過(guò)程中的質(zhì)量分?jǐn)?shù)；在10～15 d兩者的質(zhì)量分?jǐn)?shù)相差較大，這說(shuō)明由于前期農(nóng)藥質(zhì)量分?jǐn)?shù)下降后，模型的一些參數(shù)無(wú)法適應(yīng)，因此該模型模擬前期農(nóng)藥降解過(guò)程具有巨大優(yōu)勢(shì)，隨著農(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 茶葉中甲氰菊酯殘留量隨時(shí)間變化Fig.3 Changes of fenpropathrin residues in tea with time

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

如圖4所示，從元胞自動(dòng)機(jī)演化過(guò)程不同時(shí)刻元胞狀態(tài)在二維空間上的演化圖案可以看出，演化圖案隨著時(shí)間的推移不斷變淺。剛開始農(nóng)藥殘留量達(dá)到最大值，隨著時(shí)間的推移，農(nóng)藥殘留量不斷下降。農(nóng)藥降解速度分為兩個(gè)階段：快速降解階段和緩慢降解階段。剛開始降解時(shí)，由于農(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)系。可見元胞自動(dòng)機(jī)演化圖案較好地演示了農(nóng)藥殘留降解的整體過(guò)程。

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

5 結(jié) 論

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