基于元胞自动机的茶叶农药残留降解演化过程

陈洪立，屠贵林，乔 欣，彭继宇

(1.浙江理工大学 机械与自动控制学院，浙江 杭州 310018；2.浙江工业大学 特种装备制造与先进加工技术教育部/浙江省重点实验室，浙江 杭州 310023)

茶是我国最重要的饮品，茶叶的质量直接影响其产品出口。目前，随着消费者群体对身体健康的重视，其对产品质量的要求越来越高，西方发达国家茶叶农残的进口标准日趋严格，茶叶的出口阻力加大。农药残留直接影响到茶叶整体质量安全[1-2]，因此茶叶农药残留问题是当下需要迫切解决的问题之一。目前，茶叶农药残留检测的方法普遍存在费时、费力和检测精度不高等问题。通过检索国内外文献，发现从分子角度解析茶叶农药降解规律的相关研究较少，合理的分子数值模拟方法结果可以作为实验验证的依据，弥补实验工作的不足[3]。

近年来许多茶叶农药残留的研究是基于偏微分方程的模型，较为典型的有指数负增长函数模型、多项式回归模型灰色预测GM(1，1)模型等[4]。元胞自动机是一种离散的演化系统，它由许多小的组建构成，具有模拟复杂系统的能力。于乃功等[5]在形态学结构动力学模型的基础上，建立了青霉素发酵的元胞自动机模型。因此，利用元胞自动机模拟农药残留降解过程有着巨大优势。笔者采用仿真实验与数值模拟结果比较的研究方法，探索元胞自动机算法模拟农药残留降解过程的可行性。本研究以茶叶农药残留的机理和推导出的动力学微分方程为基础，建立农药残留降解的元胞自动机模型，并根据设计的规则编程设计，进行仿真实验，使得演化结果可视化。

1 农残降解的动力学微分方程

农药在植物上的衰减速度与农药的质量分数相关，质量分数越高，衰减速度越快，可用一级反动力学公式来表示，该公式是不考虑外在条件的理想化偏微分方程模型[4]，即

(1)

式中：y为t时刻的农药质量分数；t为施药后时间；k为农药降解速率常数；a为农药的初始质量分数。解微分方程(1)可得

y=ae-kt

(2)

当t取足够小，则有

(3)

式中：t=kT0，T0为离散时间间隔。由式(2～3)推导出方程的离散形式，即

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

(4)

式(4)为元胞自动机演化的设计基础。

2 元胞自动机理论在农药降解中的应用

2.1 农药残留降解中的元胞

农药降解的元胞自动机模型中，在每个时刻kT，元胞都具有不同的状态：存在或已降解。在该模型中，农药分子的残留降解过程其实就是元胞自动机的各个元胞在不同状态值之间的转变过程[6-7]。演化规则作用于每个元胞节点，它是一个状态转移函数，下一时刻元胞的状态仅与当前自身状态和邻居状态有关[8]，若当前元胞存在，则扫描其周围邻居的状态，如果周围邻居状态满足该元胞降解的条件，那么该元胞在下一时刻有一定概率降解。

2.2 元胞自动机模型建模的基本思想

农药残留降解过程与施药量、风雨和时间等多种因素有关，其消解过程同放射性物质衰变一样，符合一级动力学方程。为了更好地模拟农药残留降解过程，利用数值计算的方法对降解场景进行仿真[9-10]，对该模型作出以下假设：1) 每个元胞代表农药分子的集合；2) 根据建模需要，只考虑农药降解，不考虑中间产物的生成；3) 忽略外部因素的影响。在此假设条件下，分别对元胞、元胞空间、状态集合、元胞邻域和演化规则进行定义。最后利用计算机语言直观地模拟该降解过程。

2.3 元胞自动机的应用

根据农药残留降解过程的动力学特性，建立了相关的元胞自动机模型。采用二维结构作为农药降解的空间，演化规则根据农药降解机理以及一级动力学方程确定。模型中的每个元胞都代表特定数量的农药分子集合，它们具有不同的状态，为了研究农药残留降解速率与其质量分数的关系，采用邻居较多的邻域类型。

3 元胞自动机设计

本次模拟采用二维元胞自动机模型，模拟茶叶农药降解的二维元胞自动机可以定义为一个4 元组[7]，即

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

(5)

式中：t为离散时间；Cells为CA的基本元素，即元胞；Cellspace为元胞空间；Neighborhoods为CA的邻域，即演化规则的定义域；Rules为CA的演化规则。

3.1 CA的离散时间

对农药降解的整个过程进行时间离散化，对于每一个离散时刻t=kT0，元胞自动机模型都会统计当前农药的质量分数，随后根据所有离散时间对应的质量分数，绘制农药降解过程的农残曲线。

3.2 CA的元胞和元胞空间

元胞是元胞自动机的基本组成要素，每个元胞代表一群农药分子的集合，其结构简单，只有0和1两个状态。建立二维元胞空间，如图1所示。元胞总数N为55×55，在农药降解的元胞自动机模型中，农药的质量分数=初始质量分数×(农药元胞个数/元胞总数)。

图1 二维元胞空间Fig.1 Two dimensional cell space

3.3 CA的元胞邻域

在二维元胞自动机模型中，元胞自动机的邻居类型主要有Moore型、扩展Moore型或者标准的Von.Neumann型等[11]，如图2所示。

图2 元胞自动机邻域Fig.2 Neighbors of cellular automata

元胞邻域是其演化规则的演化作用范围，也就是元胞演化规则的定义域[12]。由于该模型需要研究农药分子与其质量分数之间的关系，因此这里采用邻居较多的Moore型邻域，即每个农药分子有8 个邻居。

3.4 CA的元胞状态

为了用数学模型再现农药分子降解动态，用数值计算的方法实现对整个复杂反应过程的模拟[13-14]，并将最后结果显示出来，该过程需要将不同时刻的元胞定义成不同的状态。自动机中每一个元胞代表一群农药分子的集合，每个元胞具有0和1两个不同状态。如果自动机元胞Cell(i，j)在t时刻的状态S(t)为0，则表示该元胞为空或已降解；如果为1，则表示该元胞仍存在。

3.5 CA的演化规则

自动机的演化规则具有局部特性，每个元胞的演化均依赖其相邻的元胞状态。也就是说，元胞的下一时刻状态只与它邻域的元胞状态有关[15]。自动机按照演化规则进行并行演化。

根据茶叶农药分子降解机理和推导出来的微分方程，设计了自动机的演化规则为

如果Sij=1，那么

(6)

式中：N(i，j，t)为t时刻自动机元胞的邻域内细胞的个数；P(i，j，t)为t时刻状态为1的元胞在下一刻演化为状态0的演化概率；Pt为t时刻状态为1的元胞在下一时刻演化为0的概率阈值，在仿真实验中根据农药残留作出必要调整；C为未降解农药和总体农药比值。

4 仿真实验

以文献[1]中的甲氰菊酯为研究对象，利用计算机语言，仿真编写该元胞自动机的演化程序。程序功能包括：1) 计算每一时刻农药质量分数的值；2) 在二维空间上，展现该模型随着时间而变化的演化图案；3) 对每个离散时刻进行农药质量分数的统计，绘制其随时间变化曲线。本研究采用的农药残留降解数据引自文献[1]中所提供的参数。模型参数分别采用表1数据。

表1 模型参数Table 1 Model parameter

表1中：初始质量分数a表示演化未开始时的农药质量分数；降解速率常速k表示甲氰菊酯降解的快慢；概率阈值Pt表示符合降解条件的元胞在下一时刻的降解概率，在仿真实验中根据农药残留作出必要调整；离散时间间隔T0即演化过程中将对每一个kT0时刻计算农药残留质量分数。将这些参数赋值给元胞自动机的演化程序并执行，分别得到农药微分方程质量分数和CA农药质量分数，结果如表2所示。

表2 不同时刻农药质量分数Table 2 Pesticide molecular concentration at different times

由表2可知：CA模型可以较好地拟合甲氰菊酯微分方程农药降解过程中的质量分数；在10～15 d两者的质量分数相差较大，这说明由于前期农药质量分数下降后，模型的一些参数无法适应，因此该模型模拟前期农药降解过程具有巨大优势，随着农药质量分数下降到一定数值后，模型参数也需要作出必要调整。根据表2分别绘制农药微分方程质量分数和CA元胞农药质量分数曲线，结果如图3所示。

图3 茶叶中甲氰菊酯残留量随时间变化Fig.3 Changes of fenpropathrin residues in tea with time

图3显示：元胞自动机演化过程中，甲氰菊酯残留量随着时间的变化曲线较好地模拟了动力学微分方程模型确定的甲氰菊酯残留量的变化曲线，同时也表现出了一定的随机性。可以看到在初期阶段CA元胞农药质量分数降解速度较快，考虑到可能是农药分子数量下降，农药降解的概率变大。随着农药分子质量分数的下降，两者的降解速度均缓慢下降。

如图4所示，从元胞自动机演化过程不同时刻元胞状态在二维空间上的演化图案可以看出，演化图案随着时间的推移不断变浅。刚开始农药残留量达到最大值，随着时间的推移，农药残留量不断下降。农药降解速度分为两个阶段：快速降解阶段和缓慢降解阶段。刚开始降解时，由于农药质量分数高，农药降解的速度较快，随着农药质量分数的下降，农药降解的速度也开始减慢，可以推断出农药降解速度与农药质量分数呈现正相关关系。可见元胞自动机演化图案较好地演示了农药残留降解的整体过程。

图4 不同时刻的农药分子残留状态Fig.4 Pesticide residues at different times

5 结 论

在农药残留降解机理及其动力学微分方程模型的基础上，建立了模拟农药残留降解的元胞自动机模型，并且对该模型进行了特性理论分析和仿真实验。特性理论分析和仿真实验结果均表明：该模型可以复现动力学微分方程所描述的农药残留降解过程，二维图像则直观地展现了农药残留降解的随机性、不确定性等演化行为。本研究是假定单一外部环境变量的情况下，从一级动力学方程角度研究的农药残留降解的元胞自动机模型。本次模拟讨论点在于农药降解的整体过程，未对降解产物、中间体等物质以及农药降解过程所发生的化学反应进行研究。因此本课题下一个研究重点是：1) 利用元胞自动机模型，更具体深入地模拟农药分子残留降解的各个过程，包括降解产物的生成，中间体以及外部因素对农药残留降解过程的影响等；2) 为了从分子角度模拟农药残留降解过程，以降解过程中的光解反应和水解反应为研究对象，以光解反应和水解反应方程定义演化规则，探究不同迭代步数情况下的农药降解形貌，分别研究光解反应和水解反应对农药残留降解的影响等。