考虑风险因素影响的引水隧洞施工方案优选

余 佳，焦 铮，苏 哲，任炳昱，王国浩，肖 尧

(天津大学水利工程仿真与安全国家重点实验室，天津 300072)

引水隧洞施工过程面临着活动时间随机变化、施工机械故障以及塌方、涌水等风险因素。然而，现有的施工方案优选研究缺乏考虑风险因素的影响，这导致无法获得合理的施工方案。因此，在考虑风险因素影响下进行引水隧洞施工方案优选具有重要的意义。考虑到隧洞施工待选方案数量庞大，有必要首先采用多目标优化方法筛选出Pareto方案集，然后综合考虑专家意见从中优选出最终方案。

在进行多目标优化问题[1]求解时，粒子群优化(PSO)算法由于结构简单、收敛速度快、参数设置少等优点，已被广泛应用[2-3]。针对传统PSO算法易陷入早熟收敛等缺点，国内外学者通过改进算法惯性权重和加速因子[4]，与人工鱼群算法[5]、模拟退火算法[6]、蚁群算法[7]等结合来优化粒子的搜索精度和效率。本文提出Levy飞行自适应混沌粒子群(LFACPSO)算法，利用Levy飞行的强跳跃特性来改进粒子的更新方式，提高粒子跳出局部最优的能力；此外，采用混沌算法来初始化粒子，并自适应地调整粒子惯性权重系数以提高粒子的局部和全局搜索能力。

国内外学者对水电地下工程和大坝工程施工方案优选进行了研究，如毕磊等[8]提出了基于数据包络分析的隧洞施工方案评价与优选方法，赵延喜等[9]利用层次分析法和模糊集法构建了TBM施工风险综合评判计算模型，何伟[10]利用灰色系统方案决策理论进行隧洞施工方案优选，张孝远等[11-13]采用模糊多属性决策法、逼近理想解法、数据驱动的D-AHP方法对大坝工程施工方案进行了优选。虽然这些研究考虑了模糊、灰色等不确定性，但是现有研究缺乏从活动时间随机变化、施工机械故障以及塌方、涌水等风险因素入手分析其对优化过程的影响，且在方案优选过程中缺乏考虑决策者的犹豫度。因此，本文提出考虑风险因素影响的引水隧洞施工方案优选方法，建立综合考虑风险因素影响的施工方案优化模型，采用LFACPSO算法求解模型获得Pareto方案集，并采用多属性群决策直觉模糊熵权幂平均(IFEWPA)方法，在考虑专家犹豫度条件下从Pareto方案集中优选出最佳方案。

1 考虑风险因素影响的引水隧洞施工方案优选方法

1.1 考虑风险因素影响的引水隧洞施工方案优化数学模型

建立考虑风险因素影响的引水隧洞施工方案优化数学模型。以考虑风险因素影响的工期和成本均值最小为目标建立目标函数。由于引水隧洞工程可能涉及多个施工标段，因此优化目标中常包含多个工期目标，目标函数g为

(1)

约束条件为

(2)

1.2 考虑风险因素影响的引水隧洞施工仿真方法

引水隧洞施工中面临着活动时间随机波动、施工机械故障以及塌方、涌水等风险的影响。采用网络计划技术(CPM)和循环网络仿真技术(CYCLONE)耦合方法[14]来模拟各类风险综合影响下的隧洞施工过程。其中，活动时间随机波动通过概率分布来表示[15]；施工机械故障风险的影响通过为钻机、自卸汽车、装载机等机械赋予故障间隔时间及故障处理时间参数，并将其考虑到钻孔和出渣活动模拟过程中；塌方、涌水等突发地质风险事件的影响通过采用贝叶斯网络[16]计算风险发生概率并将其嵌入CYCLONE模型来实现。采用贝叶斯网络计算风险事件发生概率的过程详见文献[17]。

1.3 基于LFACPSO算法的引水隧洞施工方案多目标优化方法

1.3.1 LFACPSO算法

1.3.1.1 基于Levy飞行的粒子更新方法

Levy飞行是一种非高斯随机过程，是一种短距离搜索与偶尔长距离行走相间的行走方式，其步长S服从Levy分布。目前常采用Mantegna等[18]提出的算法模拟Levy分布：

(3)

(4)

在优化方案求解过程中，当粒子停留在某一局部位置的次数大于预设值时(本文取10次[20])，则在粒子速度更新公式中引入Levy飞行，以使其跳出局部最优。基于Levy飞行的粒子速度和位置更新公式如下:

υi,t+1=ωL(xi,t)+c1ζ(pi,t-xi,t)+c2η(pt-xi,t)

(5)

xi,t+1=xi,t+υi,t+1

(6)

其中L(xi,t)=0.01Sxi,trandom(size(xi,t))

式中：xi,t——第i个粒子在第t代进化过程中的位置；pi,t——第i个粒子搜索到的历史最优位置；pt——整个粒子群搜索到的历史最优位置；υi,t+1——第i个粒子在t+1代进化过程中的速度；ω——惯性权重；c1、c2——学习因子，取c1=c2=2.0；ζ、η——符合[0,1]均匀分布的随机数；size(xi,t)——xi,t的维度；random(size(xi,t))——在xi,t各维度上生成0～1随机数。

1.3.1.2 粒子混沌初始化

混沌运动具有遍历性、随机性、对初值的敏感性等特点[21]。采用混沌序列来初始化粒子有助于提高种群多样性及粒子搜索的遍历性，克服随机初始化的不足。采用如下的Logistic混沌方程：

Zq+1=dZq+(1-Zq)

(7)

式中：Zq——第q次迭代后的混沌变量取值；d——混沌系统的控制参数，取d=4。初始值Z0∈(0,1)且Z0∉{0.25,0.50,0.75}时，系统处于完全的混沌状态。

1.3.1.3 自适应惯性权重

采用自适应调整惯性权重的策略，使惯性权重随迭代次数的改变而变化。惯性权重ω的调整方式如下：

(8)

式中：ωt——第t次迭代的惯性权重值；ωmax、ωmin——最大、最小惯性权重；tmax——最大迭代次数。

1.3.2 基于LFACPSO算法的施工方案多目标优化流程

基于LFACPSO算法的施工方案多目标优化流程如下：(a) 确定每个工序的待选方案，采用考虑风险因素影响的仿真方法计算各方案下的目标值；(b) 初始化算法参数，包括粒子数量、最大迭代次数、学习因子c1和c2、最大和最小惯性权重等；(c)采用混沌理论初始化粒子，其中一个粒子对应一种施工方案组合；(d) 计算初始粒子的适应度值，并将粒子当前的适应度值设为其初始最优适应度值；(e)选择领导粒子，判断粒子位置未更新的次数是否大于10次，若未大于10次，则采用传统的粒子更新公式[22]来更新速度和位置，否则，采用式(5)(6)来更新粒子的速度和位置；(f)选取非支配解，并根据式(8)更新惯性权重；(g)判断算法是否满足终止条件，若满足，则输出Pareto方案集合，否则，返回步骤(e)继续迭代。

1.4 施工方案多属性群决策IFEWPA方法

施工方案多属性群决策IFEWPA方法首先采用IFPWA算子，在考虑专家犹豫度的条件下进行专家意见融合，获得集成矩阵；然后采用直觉模糊交叉熵来计算属性权重，以优选出最终方案。

1.4.1 基于IFPWA算子的专家意见融合方法

将专家评价用直觉模糊数来表示，形成直觉模糊评价矩阵。第k个专家的直觉模糊矩阵Ak为

(9)

式中：αi,j,k——专家k对方案i中目标j的直觉模糊评价，αi,j,k=(μi,j,k,vi,j,k,γi,j,k)；μi,j,k、vi,j,k、γi,j,k——专家评价的隶属度、非隶属度、犹豫度，μi,j,k,vi,j,k,γi,j,k∈[0,1]；m——待选方案数量；n——优化目标数量。

采用Song等[23]提出的方法计算直觉模糊集之间的相似度，并替代P-A算子中的支持度函数，得到IFPWA算子，该算子的支持度函数Sup计算公式为

(10)

(11)

利用IFPWA算子将p个直觉模糊评价矩阵Ak(k=1,2,…,p)融合为集成评价矩阵A：

(12)

式中：ξi,j,k——各专家评价直觉模糊集之间的相互支持程度，ξ越大则相互支持程度越高。

1.4.2 基于直觉模糊交叉熵的属性权重计算方法

存在直觉模糊集A={〈x,μA(x),vA(x)〉|x∈X}和B={〈x,μB(x),vB(x)〉|x∈X}，其中μA(x)、vA(x)分别为有限非空集X中元素x对集合A的隶属度和非隶属度，μB(x)、vB(x)分别为元素x对集合B的隶属度和非隶属度。记a=1+μA(x)-vA(x)，b=1+μB(x)-vB(x)，c=1-μA(x)+vA(x)，d=1-μB(x)+vB(x)，则A与B的直觉模糊交叉熵为

(13)

由于上述直觉模糊交叉熵不满足对称性，故直觉模糊集的对称交叉熵为

C*(A,B)=C(A,B)+C(B,A)

(14)

在多属性决策过程中，属性Gj的相似度F(Gj)可表示为

(15)

F(Gj)越高，表明属性Gj越不重要，一般赋予较低的权重值。因此，对于m个方案n个属性的多属性决策问题，属性Gj(j=1,2,…,n)的权重wj可由下式计算得到：

(16)

式中：ri,j——集成评价矩阵中第i个方案下第j个目标的直觉模糊数。

确定所有属性权重wj(j=1,2，…,n)后，将m个待选方案的n个指标综合起来，得到全局偏好值Ri(式(17))，根据该值将所有待选方案排序即可获得最优方案。

(17)

式中：μij、vij——集成矩阵A中第i个方案下第j个属性的直觉模糊评价的隶属度与非隶属度。

方案i的得分函数值ti为

(18)

2 工 程 应 用

某引水隧洞C2标段采用钻爆法施工，包括1号隧洞和2号隧洞的西端洞段，1号隧洞长4 700 m，2号隧洞长6 000 m。1号、2号隧洞均分成上部和下部开挖，且两条隧洞之间每间隔500 m布置一条横通道，整个标段共设置8条横通道。采用考虑风险因素影响的施工方案优选方法进行C2标段施工方案优选。

2.1 基于LFACPSO算法的施工方案多目标优化

表1 C2标段可选机械参数

采用考虑风险因素影响的引水隧洞施工仿真方法获得各工序在不同施工机械方案下的工期及成本，然后采用LFACPSO算法进行施工方案多目标求解。对200个粒子进行混沌初始化，设置最大迭代次数为1 000次，算法学习因子c1=c2=2.0，最大惯性权重ωmax=1.0，最小惯性权重ωmin=0.2。通过算法求解获得27个Pareto方案如图1所示，图中圆的半径越大，表示成本越高。

图1 Pareto方案集Fig.1 Pareto scheme set

2.2 基于IFEWPA方法的施工方案多属性群决策

图2 专家集成评价结果Fig.2 Integrated evaluation results of experts

根据获得的各目标权重及集成的专家评价结果，由式(17)计算得到各施工方案的全局偏好值，并根据式(18)获得相应的得分函数值。由图3可知，各方案的得分函数值均不相同，且方案17的得分函数值最大，为综合最优方案。在方案17中，1号隧洞上部配置21台手风钻、10台自卸汽车、2台装载机，2号隧洞上部配置21台手风钻、12台自卸汽车、2台装载机，1号隧洞下部配置6台手风钻、2台自卸汽车、1台装载机，2号隧洞下部配置10台手风钻、3台自卸汽车、1台装载机，横通道配置6台手风钻、4台自卸汽车、1台装载机。

图3 施工方案的全局偏好值及得分函数值Fig.3 Global preference values and score function values of construction scheme

2.3 对比分析与讨论

为了证明本文所提出方法的有效性，以2号隧洞上部800 m开挖段为对象，将优选方案17下考虑风险综合影响的仿真进度与不考虑风险影响的仿真进度、实际进度、计划进度进行对比，结果如图4所示。由图4可知，在800 m段开挖过程中，相比于不考虑风险影响的仿真进度和计划进度，考虑风险综合影响的仿真进度与实际进度更加贴近。在完成800 m开挖时，实际工期为175 d，考虑风险综合影响的平均仿真工期为169 d，不考虑风险影响的平均仿真工期为155 d，计划工期为159 d，考虑风险综合影响的仿真工期更接近实际工期。

图4 2号隧洞上部施工进度对比Fig.4 Comparison of construction schedule of number 2 tunnel upper section

为了证明提出的LFACPSO算法在进行隧洞施工方案多目标优化方面的优越性，选取不带Levy飞行的自适应混沌粒子群(ACPSO)算法、传统的PSO算法和带精英策略的非支配排序遗传(NSGA-Ⅱ)算法分别进行优化模型的求解。各算法的初始种群均设置为200个，迭代次数为1 000次，每种算法分别重复计算10次，并对10次计算的结果进行对比分析，结果如图5所示。

图5 方案优化模型的不同算法求解结果对比Fig.5 Comparison of solution results of different algorithms for scheme optimization model

从算法的优化结果来看，LFACPSO算法所得的平均Pareto优化方案数量最大，为24.6个，而ACPSO算法、PSO算法、NSGA-Ⅱ算法所得的平均方案数量分别为20.6个、21个、20.7个。由此可知，相比其他算法，LFACPSO算法求解方案多目标优化模型的能力更强，所得Pareto方案集更加完备。

从算法的鲁棒性来看，LFACPSO、ACPSO、PSO、NSGA-Ⅱ算法10次计算所得Pareto方案数量的标准差分别为2.61、5.12、5.10、3.95。对比其他算法，LFACPSO算法求解结果更稳定，鲁棒性更高。

3 结 语

本文提出了考虑风险因素影响的引水隧洞施工方案优选方法。该方法首先在综合考虑活动时间随机变化、施工机械故障以及塌方、涌水等风险因素影响下建立施工方案优化模型，克服了现有施工方案优选研究缺乏考虑风险因素影响的不足。其次，针对传统PSO算法局部搜索能力差、易陷入早熟收敛的不足，提出了LFACPSO算法来求解优化模型，该算法采用基于Levy飞行的粒子更新方式来提高粒子跳出局部最优的能力，同时采用混沌算法来初始化粒子，并自适应地调整粒子惯性权重系数以提高算法的全局与局部搜索能力。最后，针对现有施工方案优选研究未能体现专家犹豫度的不足，提出了施工方案多属性群决策的IFEWPA方法，在考虑专家犹豫度的条件下优选出最终方案。相比计划进度和不考虑风险影响的仿真进度，综合考虑风险影响的仿真进度与实际进度更加吻合，证明了本文提出方法的有效性。将LFACPSO算法与ACPSO、PSO和NSGA-Ⅱ算法进行对比，结果表明，LFACPSO算法所得的Pareto方案集更加完备，结果更稳定，证明了该算法的有效性与鲁棒性。