施加松，常芸芬，吴涢晖，李发明

(防化研究院，北京 102205)

放射性烟羽扩散反问题是利用辐射环境监测数据反解放射性烟羽特征问题的总称(本文简称反问题)。由于解决反问题的技术方法可用于信息深度融合、扩散参数拟合、源项反演、核危害趋势预测等方面，因此在辐射防护、辐射环境保护、核与辐射应急等领域，关于这些方法和模型的研究一直是一个热点，尤其是2011年福岛核事故以后，反解核事故源项的技术研究达到了新的高度，大量研究方法涌现[1-7]，如最优插值法(OI)[8-10]、遗传算法(GA)[2-3]、模拟退火法(SA)[5]、人工神经网络法(ANN)[7]、卡尔曼滤波法(KF)及由其发展得到的扩展卡尔曼滤波法(EKF)和集合卡尔曼滤波法(EnKF)[4,6]、逆扩散模拟技术(IDT)[11]、贝叶斯推导算法[12-13]。这些方法各有其特点，在不同的应用领域各有其优势。美国国家大气释放咨询中心(NARAC)在充分考察多种优化算法后认为，启发式算法(包括GA、SA、ANN等)、逆向扩散算法和贝叶斯推导算法在解决源项反演问题方面都具有良好的应用条件。相对于其他算法，遗传算法具有很好的收敛性，在计算精度要求不是很高的情况下，计算时间少、鲁棒性高，其缺点是编程实现较复杂、依赖参数较多，易产生早熟收敛现象。本文拟通过优化的遗传算法来构建放射性烟羽反问题求解模型，并通过模拟数据和风洞数据的反解问题对模型进行验证，从而为扩散参数拟合和源项反演等反问题性质的工程应用提供一种高效而精确的技术方法。

1 放射性烟羽扩散反问题数学建模

(1)

其中，θ为反解参数集。

当然也可对方程(1)进行一些细致的处理，如加权重、取对数等[14]。这种基于最小二乘法思想的构建方法对数据要求较高，对于误差分布满足高斯分布的情况，计算结果一般较好，但对于其他分布或由于监测数据中有不可靠数据而造成存在非常大偏差的情况，计算结果会很不理想。实际上，由于监测数据的误差和质量问题与烟羽模型的计算误差之间的耦合关系，使得整体误差分布很难满足高斯分布，此时采用方程(1)作为目标函数，其求解结果很难保证。为解决该问题，本文通过Steiner[15]提出的最频值理论来构建目标函数：

(2)

其中：n为监测值的个数；ε为分散系数，表征误差分布分散程度；k为无量纲系数，当误差分布服从柯西分布时k=1，当分布未知时k=2，当分布中存在尖峰形态(如杰佛里斯分布)时k=3，本文取k=2。ε可通过方程(3)的迭代计算获得：

(3)

目前，放射性烟羽扩散模型主要有直线高斯烟羽模型、烟团模型、随机扩散模型等，本文采用直线高斯烟羽模型作为扩散模型，其数学表达式如下：

(4)

其中：C为下风向放射性烟羽的浓度；(x0,y0)为释放点位置；H为释放高度；Q为释放率；λ为放射性核素的衰变常量；u为释放高度为H时的风速；(σy,σz)为下风向距离x-x0处横向和垂直方向的扩散参数。

至此，放射性烟羽扩散反问题已完全转换为求解目标函数(2)的最小值的数学问题。

2 遗传算法实现

遗传算法是由美国Holland于1975年首先提出的，该算法是通过模拟达尔文生物进化论的自然选择过程搜索最优解。遗传算法已广泛应用于组合优化、机器学习、信号处理、自适应控制和人工生命等领域，是现代有关智能计算中的关键技术。

遗传算法从代表问题潜在解集的一个种群开始，按照适者生存和优胜劣汰的原理，逐代演化产生出越来越好的近似解。在每一代演化过程中，需根据个体的适应度大小选择个体(选择算子)，并借助于基因进行组合交叉(交叉算子)和变异(变异算子)，产生出代表新的解集的种群。这个过程将导致种群像自然进化一样，后代种群较前代更适应环境，末代种群中的最优个体可作为问题近似最优解。遗传算法其实是一种搜索最优解的思想，它没有给出如何选择、如何交叉、如何遗传，甚至没有给出采用什么形式来编码(形成种群中的个体)，具有非常宽泛的适用性。为使其更好地适用于放射性烟羽扩散的反问题，采用Fortran语言编写适用于方程(4)反解问题的遗传算法。具体如下。

1) 编码

以反解参数作为基因，反解参数集形成染色体(种群中个体)。

2) 适应度选取

将式(2)的倒数作为适应度，则min(F2(θ))的求解问题即变成最大值求解问题。

3) 选择算子

为在充分体现适者生存的进化法则的同时，也兼顾个体广泛性，本文采用了一种优化处理方法，具体如下。

(1) 通过式(5)计算个体在下一代群体中生存的期望数目：

(5)

其中：N为个体数；Fi为第i个个体的适应度；int函数表示取整；Ni为第i个个体进入下一代遗传选择操作的次数，通过此方法可确定∑Ni个个体进入下一代。

4) 交叉算子

为加快收敛性，本文采用精英交叉策略对选择出的下一代个体进行交叉操作。精英交叉策略的基本思想如下。

(1) 首先选择1个精英e，在此处即为适应度最大的个体。

(2) 对于每个个体mi，均随机产生1个大于0小于1的数p，如果满足p

本文对精英交叉概率pc进行改进，提出了一种自适应的精英交叉概率，即每个个体可根据其适应度来调整其概率选取值，具体公式如下：

(6)

此处，Fmax取1。通过方程(6)可看出，对于适应度较大的个体，其交叉的概率变小了，这样做的好处是保护优秀基因的个体，进一步加快收敛的速度。

5) 变异算子

变异算子是维持种群多样性、避免过早收敛到局部最优值(即“早熟”)的核心操作。本文采用当前基因值为中心、求解域为方差的高斯随机抽样的方法来获取基因突变后的个体。其中变异概率同样采用自适应的方法动态调整，其公式如下：

(7)

6) 终止条件

遗传算法终止一般有以下几种方法：(1) 采用设定最大代数终止算法运行；(2) 判定群体的收敛程度终止算法运行；(3) 通过算法离线性能和在线性能判断是否终止算法。本算法采用第1种和第2种的混合方法，设定最大代数为20 000代，或适应度值为1时结束搜索。计算程序流程图示于图1。

图1 计算程序流程图Fig.1 Computational flow chart

3 实验验证

实验验证的目的是验证该算法在放射性烟羽扩散反问题领域应用的可行性和可靠性，为进一步研究提供参考经验。验证采用2种方法：模拟实验验证和风洞实验验证。其中模拟实验验证是通过模拟实验数据来验证该算法解决放射性烟羽扩散反问题的可行性；风洞实验验证是通过风洞实验数据来验证算法的可靠性。

3.1 模拟实验验证

通过方程(4)产生1组实验数据，释放点坐标为(0,0)，监测点在x方向(西东方向)坐标50～4 000 m和y方向(南北方向)坐标-1 000～1 000 m矩形区域内以x方向80 m、y方向40 m的间隔均匀分布。释放核素选取137Cs，其他模拟参数列于表1。

表1 模拟实验参数Table 1 Simulation experiment parameter

当扩散参数的选取满足Briggs拟合曲线[16](标准曲线)时，可模拟计算出监测点的放射性浓度，其渲染图如图2所示，这是一个标准的放射性烟羽扩散的模拟结果(或称理论结果)，本文用该数据反解释放和气象参数来验证算法的可行性。

1) 反解参数

从方程(4)的直线高斯烟羽模型可知，释放参数与气象参数一起反解，会出现重解问题。因此，本文固定风速和风向来反解释放点坐标(x0,y0)、释放高度H、释放率Q和大气稳定度s(s为1～6的整数，s=1时大气强不稳定，即稳定度等级为A，s=2时大气不稳定，即稳定度等级为B，…，s=6时大气稳定，即稳定度等级为F)。因此，反解数为5个，其矢量表达形式为(x0,y0,H,Q,s)。

2) 反解结果

模拟实验中不同模拟监测点的浓度差别很大，最小浓度为2.39×10-29Bq/m3，最大浓度达2.79×106Bq/m3，为快速计算，忽略影响较小的监测点数据，使用大于1 000 Bq/m3的数据进行反解，结果列于表2。

表2中的适应值反映了扩散模型与模拟数据的适应程度，若等于1则为完全适应。从表2可看出，适应值非常接近1，这与本次模拟实验监测数据是由扩散模型计算得到的事实相符，说明该算法在解决放射性烟羽扩散反问题方面是可行的。

图2 模拟实验数据渲染图Fig.2 Contour map of simulation experiment data

表2 模拟实验反解结果Table 2 Reverse solution result of simulation experiment

3.2 风洞实验验证

通过1组风洞实验数据来验证反解算法的可靠性。风洞实验的等效参数列于表3，其中的风速为释放高度处的风速。风洞实验数据的渲染图示于图3，其中xg为下风向距离(本文为由西向东方向)，y为垂直下风向距离。本文对监测数据进行归一化处理，即将浓度除以释放率，并定义为扩散因子，单位为s/m3。

表3 风洞实验等效设计参数Table 3 Design parameter of wind tunnel experiment

1) 反解参数

实验设计的风向是西风(风向角270°)，但由图3可看出，实验烟羽中心线偏离了实验设计方案，因此，在反解参数中要考虑风向(φ)的变化。但对于烟羽模型，不同的风向会出现重解的问题，因此必须单独反解。另外，由于风洞实验的扩散特性未知，所以加入扩散参数σy和σz的拟合公式系数的反解，拟合公式如下：

σy=ayxbyσz=azxbz

(8)

其中：(ay,by)和(az,bz)分别为横向(y方向)和垂直方向(z方向)扩散参数的拟合系数。

图3 风洞实验数据渲染图Fig.3 Contour map of wind tunnel experiment data

所以，本文考虑2层反解：第1层，反解风向和扩散参数；第2层，用反解出的风向和扩散参数来反解释放点坐标(x0,y0)、释放高度H和风速u。因此，第1层的反解参数为5个，其矢量表达形式为(φ,ay,az,by,bz)，其求解条件为：

φ∈0°,360°

ay>0,by∈0,20

az>0,bz∈0,20

其中，by和bz取上限的目的是避免计算搜索时浮点溢出，可根据计算机软硬件环境调整。

第2层的反解参数为4个，其矢量表达式为(x0,y0,H,u)，其求解条件为：

H>0

u>0

2) 反解结果

(1) 风向和扩散参数的反解结果

首先用所有大于0的实验数据(全部数据)来反解，结果列于表4。从适应值看，扩散模型与数据的适应性并不理想，从数据分布图也可看出有些数据不满足高斯烟羽分布，为此去除大于1.0×10-5和小于1.0×10-8的数据(裁剪数据)，重新反解，结果列于表5。

从表5可看出，适应性有了较大提高，风向影响不大，但扩散规律变化较大。为比较其关系，对扩散参数与下风向的距离作图，并与中性条件(D类大气稳定度)的Briggs拟合的大气扩散参数曲线进行比较，结果示于图4。由图4可见，利用全部数据反解出的扩散参数曲线与Briggs的曲线变化趋势较接近，其中横向扩散在近距离(几百米以内)时较接近城市地形规律，远距离时较接近标准地形，而垂直扩散参数在相同的下风向距离时较Briggs曲线大。

表4 基于所有大于0的实验数据的风向和扩散参数反解结果Table 4 Inverse solution result of wind direction and diffusion parameters with all experimental data greater than 0

表5 基于部分数据的风向和扩散参数反解结果Table 5 Inverse solution result of wind direction and diffusion parameters with prune data

图4 横向(a)和垂直方向(b)扩散参数拟合曲线 Fig.4 Transverse (a) and vertical (b) diffusion parameter fitting curve

从以上反解结果看，风向结果受数据的裁剪与否影响不大，但对于扩散参数，数据的多少和数据的质量对结果影响很大。综合起来，全部数据反解的结果较好，本文采用其反解出的风向和大气扩散参数来反解释放点、释放高度和风速。

(2) 释放参数的反解结果

同理，首先用所有大于0的实验数据来反解，再用裁剪后的数据进行反解，结果分别列于表6、7。从反解结果看，二者有一定的差别，其中释放高度的差别较大，通过数据裁剪后，这种差别有一定的改善，这也可从适应值变大看出，由于裁剪后的监测数据与烟羽模型符合程度变高，所以反解的结果也随着改善，这说明反解结果的好坏对烟羽模型和监测数据的符合程度有较高的要求。

表6 基于所有大于0的实验数据的释放参数反解结果Table 6 Inverse solution result of release parameters with all experimental data greater than 0

表7 基于部分数据的扩散参数反解结果Table 7 Inverse solution result of release parameters with prune data

4 总结

本文采用最频值理论来建立目标函数，并与遗传算法相结合，构建了一种放射性烟羽扩散反问题求解模型；采用Fortran语言编写相应的程序，其中为提高遗传算法的搜索效率，对遗传算法中的选择、交叉和变异算子进行了优化；最后利用模拟实验数据和风洞实验数据进行反问题求解，以验证该模型的可行性和可靠性。验证结果表明：1) 该模型应用范围很广，可用于不同数据类型的多参数的反问题求解；2) 在扩散模型与监测数据适应性较高的情况下，该模型精度很高，但当适应性降低时，解的精度也降低，这是由于高斯烟羽模式本身在很多具体情况下会带来很大误差。