[0013] 本发明方法具体是：
[0014] 步骤1：设定河流区域，并测量参数以及确立水质方程。
[0015] 设定溯源实验区域(仿真假定河流区域宽200m，长100m)。
[0016] 根据具体应用场景，测量计算河流的横纵向水流速度以及污染物的扩散系数等参数，确立如下水质方程。水质方程又称之为对流扩散方程。三维的水质运输方程如下：
[0017]
[0018] 式中，C是预测点在t时刻点(x,y,z)处的污染物浓度；t是时间；x,y,z是以污染物排放点为坐标原点的坐标；k为污染物的降解系数；ux,uy,uz分别是河流的纵向水流平均速度、横向水流平均速度、垂向水流平均速度；Dx,Dy,Dz分别是污染物在水体的纵向、横向和垂向扩散系数。当河流深度远远小于河流宽度，且污染物垂向扩散的时间远远短于横向扩散时间尺度时，可以将三维水质运输方程简化为二维水质运输方程如下：
[0019]
[0020] 上式中，M0是污染物的初始点排放量，是以污染物排放点为坐标原点的坐标浓度。当污染源的位置以及初始排放总量M0已知时，就可以正向求解出指定河流区域内任意点处的浓度。我们需要根据某一时刻不同观测点处的浓度监测值或不同时刻的某一地点的浓度监测值，求出污染源位置、排放浓度和时间以及其数量。
[0021] 步骤2：具体溯源条件下污染物浓度的解析解分析。
[0022] 若污染物排放方式为连续排放。污染物连续排放在时间趋近于无穷大时，流域中各点的浓度均已达到平衡，不在变化，故与时间无关。假设初始条件为：当初始污染物投放坐标x＝0时，假定初始浓度C＝M0；当limx→∞C＝0,B为河流宽度，考虑一般有溯源需求的河流宽度都较宽，所以假设岸边的反射次数为一次，得连续排放条件下污染物浓度的解析解如下，扩散示意图如图1所示：
[0023]
[0024] 若污染物排放方式为瞬时排放。瞬时排放比连续排放多了个时间参数，假设初始条件为：当初始污染物投放坐标x＝0时，假定初始浓度C＝M0；当limx→∞C＝0,B为河流宽度，考虑一般有溯源需求的河流宽度都较宽，所以假设岸边的反射次数为一次，得瞬时排放条件下污染物浓度的解析解如下，扩散示意图如图2所示：
[0025]
[0026] 步骤3：计算观测井处的浓度观测序列。
[0027] 在步骤1设定的河流区域每隔一段距离(仿真中假设每隔100m)设定观测井，定点测量该点污染物浓度随时间的浓度时间序列，污染源与各观测井位置如图3所示。
[0028] 步骤4：溯源算法参数初始化。
[0029] 令 为算法的优化目标函数，N是决策变量的数量，xi是每个决策变量(i＝1,2,…,N)，是包含所有xi的决策变量的向量，记做 则优化问题可描述
为：
[0030] 使得 在 i＝1,2,…,N内取得最大或最小值。
[0031] 上式中，xmin,xmax分别是决策变量xi(i＝1,2,…,N)的最小、最大边界值。使用HS算法进行优化求解，此步骤中需要初始化和声记忆库的大小(HMS)、和声记忆库取值概率(HMCR)、音调微调概率(PAR)、音调微调带宽(BW)、每条和声向量的决策变量个数(NVAR)、创作的次数(Tmax)。
[0032] 根据每个决策变量的边界范围生成并初始化大小为HMS*(N+1)的和声矩阵。该矩阵包括和声库中所有的解向量与其对应的目标函数值。初始化生成的和声库(HM)如下所示：
[0033]
[0034] 步骤5：计算目标函数。
[0035]
[0036] 式中，C(xi,yi,Mi,t)是模型在观测井处的计算值， 是在已知源点浓度下的观测井处的观测值，n是观测点的个数，T是观测计算总时间。
[0037] 步骤6：迭代求解
[0038] 1.随机生成一组新的和声
[0039] 通过和声记忆库HM选择、片段调整和在HM内随机选择产生一组新的和声Xnew：
[0040] 首先，随机生成一个0到1的数r1，若r1，小于等于HMCR，则新的决策变量Xnew从和声记忆库中随机选择一组和声生成，然后对这个变量进行微调；若r1大于HMCR，则新的决策变量Xnew利用HM初始化方法根据每个决策变量的边界范围生成。其即兴生成新和声的步骤如下：
[0041]
[0042] 式中，NVAR是和声库每组和声决策变量的个数，LB、UB分别为决策变量的取值上下界。经仿真验证，NVAR取值为10时，模型性能最佳。
[0043] 2.更新和声组
[0044] 若上一步骤的新的和声xnew是由和声记忆库HM中随机生成得到，则需要评估是否需要进行动态调整以获得适应度较高的一组新和声，类似于在音乐中调音过程要改变频率，意味着以这种方式获得一种略有不同的和声。此过程由PAR参数控制，这里提出一种动态调整PAR和BW的更新策略，当随机数小于等于0.5时，随机在当前和声库中选取2组和声进行与新生成的一组和声进行线性交叉组合；当随机数大于0.5时，则选用以往基于自适应的音调微调带宽更新。更新过程如下所示：
[0045]
[0046]
[0047] 其中，单点源连续排放条件下，污染源源点浓度随迭代次数变化如图5所示，算法适应度值与迭代次数变化如图6所示。
[0048] 步骤7：求解污染源源点个数。
[0049] 迭代求解至达到迭代次数或目标函数达到误差允许范围内后，根据实际情况，自适应增加模型中污染源个数进行计算并与原目标函数做出比对，若误差小于原目标函数，则污染源源点个数加一，否则不变。溯源算法求解流程图如图4所示。
[0050] 步骤8：溯源完成。
[0051] 根据计算得出的污染源参数最优解作为模型输出。
