土遗址是古代人类活动留下的以土为主要材料的历史建筑遗迹,我国的地面遗存以西北干旱地区分布较为集中[1-3].然而,在长期自然营力和人为影响下,夯土建筑因受拉、受剪而开裂,裂缝延伸渐至贯穿,呈崩塌、失稳等破坏形式.在文物“修旧如旧”的保护理念指导下,目前普遍采用“锚固”这种隐蔽的力学稳定性控制技术处理土遗址稳定问题[4].
传统岩土锚固设计主要考虑安全性与经济性,参数取值通常较为保守,而土遗址兼具有文物和夯土建筑双重属性,要求在保证遗址安全的前提下尽可能避免由加固导致的遗址二次伤害.因此,科学的锚固参数优化方法在土遗址加固中显得尤为重要.
当前国内外对锚固参数优化的相关研究主要集中于传统岩土体锚固领域,而对于土遗址锚固参数优化方面的研究比较匮乏.谌文武等[5]以加固后交河故城崖体稳定性满足要求为约束条件,以加固性价比为目标函数,运用快速Lagrange差分分析和遗传算法,求得了锚固设计参数的最优解.姜谙男等[6]基于FLAC 3D三维数值模拟技术,应用神经网络、遗传算法建立了进化神经网络演化有限差分方法,给出了确定约束条件的优化方法与步骤.尹顺德等[7]运用神经网络建立了设计参数与稳定性及造价间的非线性映射关系,以造价为目标函数,利用并行遗传算法在全局范围搜索出最优锚固设计参数.
可以看出,现有锚固参数优化方法主要以遗传算法、神经网络和有限差分法为基础,计算量巨大且难以收敛,不利于大规模工程应用.此外,实际工程中遇到的问题主要为多重目标响应优化,此类优化问题较为复杂,而单一目标响应仅能对最简单的优化问题进行建模[8-9].因此,如何实现多目标需求下土遗址锚固参数设计的定量化与简易化是需要进一步研究的课题.
本文引入统计中的期望函数法,对土遗址锚固系统进行了单一响应与多重响应优化,目的是在满足锚固力需求的前提下,尽量减小钻孔对遗址基体的破坏(本文在锚杆直径为定值的前提下采用注浆量大小间接表示钻孔对遗址的伤害程度).试验设计为全因子试验,取两个关键锚固参数(锚固长度与锚孔直径)各5种水平,利用响应面法构建模型,而后采用期望函数法寻找满足目标响应需求的最优参数水平组合,从而获得最大锚固力与最小遗址伤害的良好平衡.
因子设计是刻画多变量过程的有效工具,其可以区分重要因素与次要因素,并识别它们之间可能存在的相互影响.不同的因子设计规划由所有因素在不同水平下的全部组合构成.
RSM是一种建模分析工程问题的数学与统计的综合技术,它有助于解决工程中的建模、分析、拓展、改进与优化等问题,也是构造优化模型的有效工具[10].RSM由以下三部分组成:1)用试验方法初步探索输入因素或过程是否有改进空间;2)用经验统计模型去构建过程变量与响应之间的适当近似关系;3)用优化方法寻求过程变量的最优水平(值)组合以获得满意的响应输出值[11].RSM同样可以帮助量化单个或多个已测得的响应和关键输入因素之间的关系.
RSM的第一步是应用全因子设计(full factorial design)或部分因子设计(fractional factorial design)方法构建试验,如Box-Behnken试验设计或中心复合设计等.在许多工程领域中,某一输出变量y与一组可控变量x1,x2,…,xn之间存在某种关系,而在多数RSM问题中,独立变量(影响因素)与响应之间也存在某种不确定的近似关系.因此,RSM的下一步就是为某一响应与一组独立变量之间的真实函数关系寻找一个适当的近似方法,通常采用高阶多项式(如二阶多项式)模型进行分析.
例如,当影响因素数量为n时,响应y的二阶线性回归模型形式如下:
该模型包含了多种回归系数,其中,β1,β2,…,βn为主效应系数;β11,β22,…,βnn为二次主效应系数;β12,β13,…,βn-1,n 为双因素相互作用效应系数.
Derringer和Suich[12]在1980年应用同步优化技术对多重响应进行了优化分析.该过程介绍了期望函数法的概念:利用目标函数D(x)(称为复合期望函数),并将预期响应转化为无量纲量d*(独立期望指标),期望区间为0~1(0为最不期望,1为最佳期望).使D取得最大值时的参数设定即为最优参数组合.
①采用文献[12]提出的公式计算与各目标响应相对应的无量纲量d*(单一目标响应合意性).期望函数依据响应特征有3种形式:望目特性、望小特性和望大特性.本小节依据工程实际需求(最大锚固力P与最小注浆量V),采用望大特性来确定最大锚固力P的值,采用望小特性来确定注浆量V的值(i和k由实际工程中望大、望小响应的数量决定,且0≤i≤n,0 ≤ k≤ n,i+k=n,本文取 i=1,k=1) .
图1 技术路线
Fig.1 The technical route
(a)望大特性,即yi的值越大越好.当yi超过某一特定值时认为最佳,期望值为1;当yi小于某一特定值时认为不可接受,期望值为0;当yi介于前述两者之间时,期望值即为计算值.望大特性的期望函数形式可写为
式中T为响应yi的第i个目标值,L为该响应可接受的最小值,wi为权重.当wi=1时,函数为线性;当wi>1时,则获得yi的目标响应更为重要;当wi<1时,为获得响应yi的目标值所分配的权重较小.
(b)望小特性,即yk的值越小越好.当yk小于某一特定值时认为最佳,期望值为1;当yk超过某一特定值时认为不可接受,期望值为0;当yk介于前述两者之间时,期望值即为计算值.望小特性的期望函数形式可写为
式中T为响应yk的第k个目标值,U为该响应可接受的最大值.
②计算多目标响应的复合合意性D.D的值可由全部n个目标响应的d*值组合计算得到,即
式中wn为响应n的重要性系数,1≤n≤i+k.wn的值越大,说明第n个响应相对于其他响应更为重要.下一步分析的目的就是要寻找一种最优的响应设定以获得最大的D值.
③确定最优参数水平组合.D值越大则过程质量越好.因此,在已获得D值的基础上,可对每个可控参数的影响及最优水平进行估计.
本文首先设计全因子试验,采用多元线性回归建立模型,利用响应面方法研究了控制因素与单一目标响应的关系,而后将统计中的期望函数法引入多重目标响应优化中,获得了多重目标响应需求下的锚固参数可行域范围.技术路线如图1所示.
课题组前期对遗址土、楠竹锚杆、改性泥浆等材料的物理力学性能进行了详细测试(试验方案如图2所示)[3],并通过现场拉拔试验给出了楠竹-改性泥浆界面的黏结-滑移模型(如图 3 所示),其中主要模型参数为 τ1=1.89 MPa,τ2=0.44 MPa,s1=3.34 mm,s2=10.73 mm,s3=22.56 mm[13].此外,通过对不同影响因素水平下锚固系统的极限锚固力、破坏模式、滑移量等的分析,明确了锚固长度、锚固剂厚度、锚固剂强度、锚杆表面粗糙度及布设角度等均会影响锚固系统的锚固力,其中锚固长度对锚固性能的影响最为显著,其次为锚杆表面粗糙度和锚固剂厚度的影响[13-14].
图2 锚固系统拉拔试验示意图
Fig.2 Schematic diagram for the pull-out test of the anchorage system
图3 楠竹-夯土锚固界面黏结-滑移模型
Fig.3 The bond-slip model for the bamboo-rammed earth interface
本文锚固参数优化的目的是为了获得最大锚固力与最小遗址伤害的良好平衡.各影响因素中对锚固力影响最大的因素为锚固长度,而对遗址伤害程度影响最大的是锚孔直径.因此,本文选取锚固长度Ld与锚孔直径d作为控制因素,应用52次全因子试验来分析锚固长度Ld与锚孔直径d对锚固力P与注浆量V(在锚杆直径d0一定时,V值间接表示了钻孔对遗址的伤害程度)的影响.可控因素分别为锚孔直径d与锚杆直径d0之比d/d0和锚固长度Ld与锚杆直径d0之比Ld/d0.为减小误差,锚杆均选用外直径约35 mm的楠竹制成,试样所采用的锚固剂均按土、水泥、粉煤灰的质量比为85∶10∶5配制,水灰比0.31,在此基础上加入浓度5%的硅溶胶制成,且每种因素水平组合均设2个试样重复试验.各控制因素及其取值水平如表1所示.每个试样测得的最大锚固力P及相关参数列于表2,由于每种组合均进行了2个试样重复试验,因此每组参数均对应2个P的试验值,其中V=Ldπ(d)/4.
注1 d,d0与前述无量纲量及D有本质区别.
表1 控制因素与因素水平
Table 1 Control factors and factor levels
control factor factor level 1 2 3 4 5 Ld/mm 800 1 200 1 500 2 000 3 000 d/mm 75 80 85 100 110 Ld/d0 22.86 34.29 42.86 57.14 85.71 d/d0 2.14 2.28 2.43 2.86 3.14 anchor diameter d0/mm 35(mean)
表2 楠竹锚杆现场拉拔试验结果[13-14]
Table 2 Bamboo bolt in-situ pull-out test results[13-14]
Ld/d0 d/d022.86 34.29 42.86 57.14 85.71 P/kN V/mm3 P/kN V/mm3 P/kN V/mm3 P/kN V/mm3 P/kN V/mm3 2.14 9.88 13.22 2.76E+6 15.36 13.28 4.14E+6 16.14 14.63 5.18E+6 19.11 21.09 6.91E+6 23.54 20.03 1.04E+7 2.28 14.61 17.93 3.25E+6 22.19 19.79 4.87E+6 27.38 31.33 6.09E+6 29.89 33.14 8.12E+6 32.98 34.27 1.22E+7 2.43 19.39 19.39 3.77E+6 28.36 25.22 5.65E+6 35.72 34.47 7.07E+6 39.82 38.44 9.42E+6 40.50 38.63 1.41E+7 2.86 22.08 19.87 5.51E+6 29.94 30.17 8.27E+6 34.38 36.99 1.03E+7 40.87 38.12 1.38E+7 37.17 39.65 2.07E+7 3.14 18.98 19.41 6.83E+6 29.84 28.10 1.02E+7 33.13 35.19 1.28E+7 36.81 29.83 1.71E+7 34.21 38.74 2.56E+7
将试验中参数d/d0和Ld/d0在各水平下对应的P值和V值应用统计方法进行回归分析,结果以回归方程的形式表达.本文选取二次式对回归模型的显著性及失拟项进行了测试.在未对响应进行任何变换的基础上,拟合结果显示两个响应的二次式模型在统计上是显著的(拟合优度R2>80%,且否定原假设的概率P<0.05),因此可以用其进行进一步分析.模型中的某些影响较小的项可以被忽略,从而获得简化的模型.失拟项对两个响应的影响均不显著,说明失拟项与纯误差的相关性较小.以下为依据试验数据获得的响应P和V的二元二次多项式线性回归模型:
为了更清楚地了解各过程变量对响应P和V的影响规律,根据模型公式(式(5)和(6))分别绘制了各单一响应的预测响应面图和等值线图(图4和图5).
图4为用过程参数表示的响应P的响应面图.可以看出,响应P随d/d0的增大呈先增后减的变化趋势,而随Ld/d0的增加P值呈非线性增长.
为分析响应P对控制参数的敏感程度,绘制了响应P的等值线图(图5).据图可看出,Ld/d0和d/d0对P值的影响均较为显著,响应P的最优区域位于图形顶部偏右,此区域内最大拉拔荷载P的值将超过42 kN.
图4 P的响应面图
Fig.4 The response surface diagram of P
图5 P的等值线图
Fig.5 The contour plot of P
图6 和图7分别为响应V的响应面图和等值线图,两图均显示V的值将随d/d0和Ld/d0的增大而增大,同时当d/d0取值较大时,V随Ld/d0的增幅也更大.响应V随d/d0和Ld/d0的总体变化趋势可依据响应面图形状与等值线图得到,V的最优区域位于图7的左下角,即当d/d0和Ld/d0均较小时响应V取得最优值,此区域内V值将小于5×106 mm3.
上述分析明确了各响应随不同过程参数的变化趋势,确定了最优响应对应的大致区域,但单一响应分析中各过程参数间的相对重要性依然不明确,响应的最优值及响应的最优参数组合仍需进一步确定.图8为根据试验结果获得的P的最大响应值与相应的最优因素组合.响应P 的最大化搜索范围为 9.88~40.87 kN .
由图8可以看出,响应P随d/d0和Ld/d0均呈先增大后减小的变化趋势,P的预测最大值为43.04 kN,该值超过了目标值40.87 kN,因此d*值为1,相应的过程参数d/d0和Ld/d0的最优取值组合为 2.766 3 和 73.013.
响应V的最小值可通过统计分析获得.依据试验结果可知,响应V的最小化搜索范围为2.76×106~2.56×107 mm3.最小响应值与相应的最优因素组合如图 9 所示.
图6 V的响应面图
Fig.6 The response surface diagram of V
图7 V的等值线图
Fig.7 The contour plot of V
图8 响应P的最优参数组合
Fig.8 The optimal parameter combination of response P
由图9可以看出,响应V随d/d0和Ld/d0的增大而增大,变化趋势近似线性,V的预测最小值为 2.875×106 mm3,该值略大于目标值 2.76×106 mm3,因此 d* 的最大值为 0.994 95,为在d/d0和Ld/d0均取最小值时获得.
综上所述,可以发现,响应P的最优值在d/d0和Ld/d0均较大时获得,而响应V的最优值在两个参数取值均最小时获得,二者所对应的最优参数取值间存在冲突.因此,需要寻找一组相对均衡的参数组合使两个响应的值均能满足需求.
图9 响应V的最优参数组合
Fig.9 The optimal parameter combination of response V
本节将采用多重响应优化方法来克服单一响应优化中出现的参数取值冲突问题.在此方法中,需要为所有响应赋予期望权重.由于在土遗址文物保护工程中,锚固后遗址的安全稳定与遗址本体的原样保持同样重要,因此足够的锚固力与遗址本体“最小干预”两个目标同等重要,故本文采用同量加权,即取w=1.对于所有响应的联合影响,期望取决于不同的参数值.本节采用基于RSM的数学模型,分析过程参数对D值的影响.下式为根据试验数据采用二元二次多项式线性回归获得的复合合意性D的回归模型:
图10为用过程参数表示的复合合意性D的响应曲面图,可以看出D随过程参数呈非线性变化:随d/d0和Ld/d0的增大,D值均呈先增大后减小的变化趋势,曲线均呈“凸”状.图11为复合合意性D的敏感度等值线图.D的最优区域位于图形中部,该区域内D值将大于0.8,即d/d0,Ld/d0过大或过小时均不能获得满意的D值.
将所有单一响应的最优过程参数组合均绘制于同一等值线图中,得到多重响应的重叠等值线图,该图有助于设计者依据不同区域的响应值和目标需求可视化地对参数取值进行取舍.对每一种响应,应先根据该响应的实际需求设定可接受的响应值区间(高和低的等值线),然后绘制重叠等值线图.图12所示为响应P和V的重叠等值线图,图中实线表示下限,虚线表示上限.可行域(图12中白色区域)是两个响应上下限内的重叠区域,该区域内响应P和V的可接受值均位于各自的等值线之间,可行域内任意点对应的过程参数水平组合均能满足前期设定的响应需求.例如 Ld/d0=53 和 d/d0=2.5;Ld/d0=42 和 d/d0=2.7.
上述分析明确了复合合意性D随不同过程参数的变化趋势,确定了最优响应的大致区域,但为了更准确地确定各过程参数水平的最优组合,仍需对多重响应各参数间的相对重要性进行分析.假定设计者期望获得P的目标值为40.87 kN,尽管任何大于40.87 kN的P值均认为最佳;期望获得V的目标值为2.76×106 mm3,尽管所有小于2.76×106 mm3的V值均认为最佳.
图10 D的响应面图
Fig.10 The response surface diagram of D
图11 D的等值线图
Fig.11 The contour plot of D
图12 响应P和V的重叠等值线图
Fig.12 The overlaid contour plot of P and V
表3 各过程参数的取值区间
Table 3 The interval value of each process parameter
parameter acceptable lower limit L acceptable upper limit U weight w d/d0 2.14 3.14 1 Ld/d0 22.86 85.71 1
表4 各响应的期望目标及取值区间
Table 4 The interval value and desired goal of each response
response goal acceptable lower limit L minimun acceptable upper limit U maximum weight w P/kN 9.88 40.87 1 V/mm3 2.76E+6 2.56E+7 1
表3和表4中列出了各过程参数与目标响应的取值区间,其中过程参数取值区间由表1设定的因素水平确定,目标响应取值区间由锚杆现场拉拔试验结果(见表2)确定.本文采用这些设定,利用统计方法对该优化问题进行求解.
采用统计软件(SPSS、MINITAB 等)分析可得,当 d/d0=2.604 6,Ld/d0=48.888 8 时复合合意性最佳.将其分别代入式(5)和(6)即可得预测响应P和V的最优解,其值分别为37.780 3 kN和9.5×106 mm3,此时无量纲量d*的计算公式可以化为
由于 9.88 kN≤ P max≤40.87 kN,d*1=(37.780 3-9.88)/(40.87-9.88)=0.900 3.
由于 2.76×106 mm3≤ V min≤2.56×107 mm3,d*2=(2.56×107-9.5×106)/(2.56×107-2.76×106)=0.704 9.
注意到本文需考虑两种响应同等重要,因此式(8)和(9)中权重w1,w2均取为1,则这两个期望函数均为线性.复合合意性D可由式(4)计算获得,有D=(0.900 3×0.704 9)1/2=0.796 6.
图13 响应P和响应V的最优参数组合
Fig.13 Optimal parameter combination of response P and V
图13 为通过上述方法确定的响应P和V的优化设置与优化结果.该图显示了不同因素水平组合对各目标响应产生的影响趋势.图中每一列对应一种因素,首行对应复合合意性D,其余各行均对应一个响应变量(本文中分别为P和V).每列顶部数字分别为试验设计中的当前因素设定水平及其上下限值;每个响应行的左侧为当前因素组合下该响应的预测值和d*值.当前设定下的响应D值位于图13左上角.
图13中首行的每个单元格均显示了响应D随某一过程参数的变化趋势(其余参数固定不变).图中竖实线为当前参数设定值,横虚线表示当前响应值.优化结果显示:最优因素水平组合下(d/d0=2.604 6,Ld/d0=48.888 8) 响应P 对应的d*1 值为0.900 3,响应V对应的d*2 值为0.704 9,两个响应的复合合意性D值为0.796 6.可见,为获得锚固力最大化目标与注浆量最小化目标的良好平衡,两类响应的取值均作出了一定让步,因此其各自所对应的无量纲量d*值均小于1,据此计算获得的复合合意性D值约为0.8,可认为此锚固参数组合(锚孔直径为91 mm,锚固长度为1.7 m)能够较好地满足各个响应的需求.
在实际工程中,可采用类似方法,通过调整优化图中的不同参数水平来满足响应的不同需求.如在尽可能保证最优性能的前提下寻求最小造价;或根据需要的参数组合来预测响应值.
本文将期望函数法引入土遗址锚固参数优化中,分析了以锚固力P最大化和遗址伤害最小化(以注浆量V反应)为目标时锚固参数的组合优化方法,主要结论如下:
1)采用响应面模型可以明确某一独立锚固参数(锚固长度、锚孔直径)对特定目标响应(锚固力最大、遗址伤害最小)的影响规律;采用期望函数法能够在满足可接受目标响应需求的条件下获得锚固参数的最优组合.
2)基于期望函数法的单一响应优化结果显示:响应P的最优值在d/d0和Ld/d0均较大时获得,而响应V的最优值在两者取值均最小时获得,两种响应所对应的最优参数取值间存在冲突.
3)基于期望函数法的多重响应优化结果显示:当d/d0和Ld/d0取值均适中时复合合意性D值最大,该方法同时能够方便地确定多重目标响应需求下的锚固参数可行域范围,方便设计者根据实际工程条件对锚固参数进行可视化取值.
[1] 孙满利,王旭东,李最雄.土遗址保护初论[M].北京:科学出版社,2010.(SUN Manli,WANG Xudong,LI Zuixiong.Initial Discussion of Earthen Sites[M].Beijing:Science Press,2010.(in Chinese))
[2] FODDE E,WATANABE K,FUJII Y.Preservation of earthen sites in remote areas:the buddhist monastery of Ajina Tepa,Tajikistan[J].Conservation and Management of Archaeological Sites,2013,9(4):194-218.
[3] ZHAO D,LU W,WANG Y L,et al.Experimental studies on earthen architecture sites consolidated with BS materials in arid regions[J].Advances in Materials Science and Engineering,2016(2):1-13.DOI:10.1155/2016/6836315.
[4] 张景科,郭青林,李最雄,等.土遗址锚固机理初探[M].兰州:兰州大学出版社,2014.(ZHANG Jingke,GUO Qinglin,LI Zuixiong,et al.Preliminary Study on Anchorage Mechanism of Rammed Earth Sites[M].Lanzhou:Lanzhou University Press,2014.(in Chinese))
[5] 谌文武,张宇翔,和法国,等.基于FLAC和遗传算法的斜坡加固方案优化方法[J].中南大学学报(自然科学版),2011,42(11):3507-3514.(CHEN Wenwu,ZHANG Yuxiang,HE Faguo,et al.Optimization method for slope reinforcement design based on FLAC and genetic algorithms[J].Journal of Central South University(Science and Technology),2011,42(11):3507-3514.(in Chinese))
[6] 姜谙男,冯夏庭,刘建,等.基于三维数值模拟的地下大型洞室锚固参数智能优化[J].岩石力学与工程学报,2004,23(10):1700-1705.(JIANG Annan,FENG Xiating,LIU Jian,et al.Intelligent optimization of anchoring parameters for large underground houses based on three dimensional numerical simulation[J].Chinese Journal of Rock Mechanics and Engineering,2004,23(10):1700-1705.(in Chinese))
[7] 尹顺德,冯夏庭,张友良,等.滑坡加固方案优化的并行进化神经网络方法研究[J].岩石力学与工程学报,2004,23(16):2698-2702.(YIN Shunde,FENG Xiating,ZHANG Youliang,et al.Study on optimum design of landslide stabilization by parallel evolutionary neural network method[J].Chinese Journal of Rock Mechanics and Engineering,2004,23(16):2698-2702.(in Chinese))
[8] GHAZY M F.Effect of using mortar interface and overlays on masonry behavior by using Taguchi method[J].ACI Material Journal,2012,109(5):509-516.
[9] AKALIN O.Eco-cement optimization using statistical mixture design method[J].ACI Materials Journal,2014,111(4):391-398.
[10] NARAYANAN N,RAMAMURTHY K.Prediction relations for compressive strength of altered concrete[J].ACI Materials Journal,2000,97(3):367-373.
[11] MONTGOMERY D C.Design and Analysis of Experiments[M].New York:John Wiley &Sons,2011.
[12] DERRINGER G,SUICH R.Simultaneous optimization of several response variables[J].Journal of Quality Technology,1980,12(4):214-219.
[13] LU W,ZHAO D,MAO X F,et al.Experimental study on bond-slip behavior of bamboo boltmodified slurry interface under pull-out load[J].Advances in Civil Engineering,2018(5):1-23.DOI:10.1155/2018/6960285.
[14] 毛筱霏.高昌故城遗址保护与锚杆加固应用技术研究[D].博士学位论文.西安:西安建筑科技大学,2009.(MAO Xiaofei.The protection of Gaochang site and research on consolidated technology of bolt[D].PhD Thesis.Xi’an:Xi’an University of Architecture and Technology.(in Chinese))
A Composite Optimization Method for Anchorage Parameters of Rammed Earth Sites Based on Desirability Functions
芦苇,赵冬,李东波,毛筱霏.基于期望函数的土遗址锚固参数组合优化方法[J].应用数学和力学,2019,40(6):650-662.
LU Wei,ZHAO Dong,LI Dongbo,MAO Xiaofei.A composite optimization method for anchorage parameters of rammed earth sites based on desirability functions[J].Applied Mathematics and Mechanics,2019,40(6):650-662.