随着互联网智能控制技术和计算机仿真技术的快速发展，我国特高拱坝建设正由“数字大坝”^[1]向“智能大坝”^[2-3]进行转变，这带来了温控防裂技术和管理水平的提高^[4]。在温度监测方面，分布式光纤应用更为广泛^[5-6]。然而要实时了解特高拱坝温度性态，仅靠光纤监测数据还是不够的，还需通过数值仿真对温度未来变化趋势进行预测。在特高拱坝温度分布状态预报方面，关键在于确定边界条件和热学参数^[7]，目前施工期热学参数多选用室内试验值^[8]。实际上，混凝土施工期温度变化是自身水化热、通水冷却、环境气温、流水养护和表面保温等共同作用的结果，热学参数采用室内试验值是不准确的。因此，在特高拱坝温度仿真计算过程中，有必要基于大坝温度观测值和群体智能优化模型，反演获得真实浇筑环境下的混凝土热学参数^[9-10]。

在混凝土热学参数智能识别方面，王振红等^[11]通过现场试验和遗传算法进行反演；黄达海等^[12]结合室外试验结果，反演获得混凝土导温系数；苏怀智等^[13]考虑水管冷却效果，基于复合形法进行反演；强晟等^[14]结合现场试验，反演获得了混凝土表面散热系数；肖照阳等^[15]基于现场资料和BP神经网络进行反演；周建兵等^[16]考虑一期冷却期间多挡通水，基于单纯形法进行反演；Liu等^[17]基于现场实验，将改进的粒子群算法用于Dittus-Boelter方程参数识别。上述研究通常是将混凝土视为单一均质材料，实际上混凝土温度是多种材料共同作用的结果。对于不同浇筑仓来说，强度等级不同，混凝土水化热特性不同。对于同一浇筑仓来说，通常使用多种级配类型的混凝土，级配不同，胶凝材料用料、组分和水灰比等参数存在差异，混凝土水化热亦不同。因此有必要对不同强度混凝土和同一强度不同级配混凝土进行反演，确定待反演热学参数，这其中亦需要考虑长浇筑历时带来的环境气温差异、冷却通水动态变化。

粒子群算法^[18-20](particle swarm optimization, PSO)通过种群粒子之间“探索”和“开发”来求得最优解。由于粒子同一趋向于最优解，局部搜寻能力较差。为了避免算法陷入局部最优以及进行高效搜索，学者们对惯性权值进行了优化研究，基本思想就是惯性权值非固定值，以曲线形式进行递减变化。Shi等^[21]给出了惯性权值线性递减策略；陈贵敏等^[22]提出了3种非线性惯性权值递减策略，发现大多数情况下凹函数权值递减策略最优；张俊等^[23]给出了自适应指数惯性权值。

为了提高算法的精度和收敛速度，本文将遗传算法中的交叉和变异操作引入粒子群优化算法，交叉算子在粒子群中随机选择两个粒子进行速度和位置交换组合，父串的优秀特征可以遗传给子串，该操作可以使搜索范围扩大，全局搜索能力增强。变异算子是以事先设定的概率判断粒子是否进行变异，若变异对粒子位置作变动，从而使粒子具有局部搜索能力，维持粒子群多样性。最终建立了混合粒子群优化模型(hybrid particle swarm optimization, HPSO)，凹函数权值递减策略代替线性递减惯性权值^[24]，放弃传统PSO通过追踪粒子群极值来更新粒子位置的思想。相反，通过与个体极值和全局极值交叉并在粒子群中变异来寻找最优解，从而提高算法的全局和局部搜索能力。李红亚等^[25]对遗传算法和粒子群算法两者融合方式进行了分析，给出了并行式、串行式(先遗传算法后粒子群算法、先粒子群算法后遗传算法)、嵌入式(遗传算法嵌入、粒子群算法嵌入)3种混合方法，并指出混合粒子群算法克服了传统粒子群算法易局部最优的缺点，增强了遗传算法收敛精度。目前，HPSO已应用于混凝土坝黏弹性参数反演求解^[26]和配电系统分布式电源优化选址^[27]。Gao等^[28]研究表明，与标准遗传算法和模拟退火算法相比，HPSO更有效。

本文基于分布式光纤测温数据，将混合粒子群优化算法应用于白鹤滩拱坝低温季节C40和C35混凝土、C40不同级配混凝土热学参数反演，考虑冷却水换挡以及沿水管水温变化影响，选用实际通水冷却方案以及环境气温，从而实时跟踪混凝土发热规律。

1 混合粒子群算法 1.1 标准粒子群算法

作为一种群体随机优化方法，粒子群算法主要用到两个概念：一是全局搜索，指的是粒子脱离原先路线向新轨迹搜索的能力；二是局部搜索，指的是粒子在原先搜索轨迹上更精细搜索的能力。

假定在D维搜索空间有n个粒子，第i个粒子位置为X_i=(X_i1, X_i2, …, X_iD)^T，速度为V_i=(v_i1, v_i2, …, v_iD)^T，个体极值为P_i=(p_i1, p_i2, …, p_iD)^T，粒子群全局极值为P_g(p_g1, p_g2, …, p_gD)^T。

粒子速度和位置可通过式(1)和(2)更新为

$ v_{i d}^{k+1}=\omega v_{i d}^{k}+c_{1} r_{1}\left(p_{i d}^{k}-X_{i d}^{k}\right)+c_{2} r_{2}\left(p_{g d}^{k}-X_{i d}^{k}\right), $

(1)

$ \begin{array}{c} X_{i d}^{k+1}=X_{i d}^{k}+v_{i d}^{k+1} ,\\ d=1,2, \cdots, D, \quad i=1,2, \cdots, n . \end{array} $

(2)

式中：ω为惯性权重；c₁和c₂为学习因子，其值非负；r₁和r₂为[0, 1]内的随机数。式(1)中ωv_id^k表示粒子保持当前速度下的全局搜索能力，$ {c_1}{r_1}\left( {p_{id}^k - X_{id}^k} \right)和 {c_2}{r_2}\left( {p_{gd}^k - X_{id}^k} \right) $则表示粒子局部搜索能力。

1.2 凹函数惯性权重系数

惯性权重ω为粒子保持运动惯性的能力，ω值大则全局搜索能力较强，ω值小则局部搜索能力较强。在算法迭代过程中，有必要依据不同时期搜索情况对ω值进行调节。本文选用的凹函数递减惯性权值为

$ \omega(k)=\omega_{\max }-\left(\omega_{\max }-\omega_{\min }\right)\left(\frac{2 k}{T_{\max }}-\left(\frac{k}{T_{\max }}\right)^{2}\right). $

(3)

式中：k为当前迭代次数，T_max为允许最大进化次数。在大多数应用中^[29]，ω_max=0.9，ω_min=0.4。

图 1给出了凹函数递减惯性权值变化曲线，可以看出，在进化前期，凹函数递减权值比线性递减权值减少更快；在进化后期，相比线性递减权值，凹函数递减权值减少放慢，更有利于粒子跳出局部最优而获取全局最优。

1.3 粒子交叉及变异策略

基于交叉及变异的粒子群算法，是指对个体极值交叉并在粒子群中变异产生新的种群来寻找最优解。该方法相比于传统的粒子群算法，摒弃了通过追踪粒子群极值来更新粒子位置的思想，提高了粒子群多样性和局部寻优能力，算法流程如图 2所示。

该方法在每次进化过程中，根据交叉率p₁选择给定数量的粒子放入交叉池中，池中的粒子随机两两杂交，其速度和位置更新如下：

$ v_{1}^{\prime} =\frac{v_{1}+v_{2}}{\left\|v_{1}+v_{2}\right\|}\left\|v_{1}\right\| , $

(4)

$ v_{2}^{\prime} =\frac{v_{1}+v_{2}}{\left\|v_{1}+v_{2}\right\|}\left\|v_{2}\right\| , $

(5)

$ X_{1}^{\prime} =p_{1} X_{1}+\left(1-p_{1}\right) X_{2}, $

(6)

$ X_{2}^{\prime} =p_{1} X_{2}+\left(1-p_{1}\right) X_{1} . $

(7)

依据概率p₂采用Gauss算子，在粒子群中选取一定数量的粒子进行变异。

$ X_{i}^{j}=\frac{X_{i}^{j, \max }+X_{i}^{j, \min }}{2}+\frac{X_{i}^{j, \max }-X_{i}^{j, \min }}{6}\left[\sum\limits_{n=1}^{12} R_{n}-6\right]. $

(8)

式中：i为粒子编号，j为粒子第j维位置，X_i^{j, max}和$ X_i^{j, \max }和 X_i^{j, \min } $分别为第i个粒子j维位置最大值与最小值，R_n∈(0, 1)服从均匀分布。因此，$ \left[ {\sum\limits_{n - 1}^{12} {{R_n} - 6} } \right] $为均匀分布的随机数，用于替代满足正态分布(均值为$ \frac{{X_i^{j, \max } + X_i^{j, \min }}}{2} $，方差为$ \frac{{X_i^{j, \max } - X_i^{j, \min }}}{6} $的随机数。

1.4 混合粒子群算法性能测试

选取10维Girewank、Rastrigrin、Rosenbrock和Sphere函数进行寻优性能测试，试验100次，粒子数取20，进化100代，计算结果如表 1所示。

很明显，对于这4个测试函数，HPSO陷入次优解次数更少，更有利于跳出局部最优解寻得全局最优解。对于收敛精度和收敛速度，HPSO所得结果均优于PSO。

2 混凝土热学参数优化反演 2.1 初期有热源水管冷却问题

混凝土绝热温升通常采用双曲线式进行表达，然而不易进行计算。为了计算方便并提高计算精度，朱伯芳^[7]提出了组合指数式，表示如下：

$ \theta(\tau)=\theta_{0} s\left(1-\mathrm{e}^{-m_{1} \tau}\right)+\theta_{0}(1-s)\left(1-\mathrm{e}^{-m_{2} \tau}\right) . $

(9)

绝热温升双曲线式表示如下：

$ \theta(\tau)=\theta_{0} \frac{\tau}{n+\tau}. $

(10)

式(9)和式(10)中参数之间关系表示如下：

$ s=0.6, m_{1}=\frac{1.45}{n}, m_{2}=0.1 m_{1}. $

(11)

在一期通水冷却阶段，往往换挡多次，此时混凝土的平均温度为^[16]

$ T(t)=T_{w i}+\left(T_{i}-T_{w i}\right) \phi_{i}(t)+\theta_{0} \psi_{i}(t), $

(12)

$ \phi_{i}(t)=\mathrm{e}^{-p_{i} t}, $

(13)

$ \begin{array}{l} \psi_{i}(t)=\frac{s m_{1}}{m_{1}-p_{i}}\left(\mathrm{e}^{-p_{i} t-m_{1} t_{i}}-\mathrm{e}^{-m_{1}\left(t+t_{i}\right)}\right)+ \\ \ \ \ \ \ \ \frac{(1-s) m_{2}}{m_{2}-p_{i}}\left(\mathrm{e}^{-p_{i} t-m_{2} t_{i}}-\mathrm{e}^{-m_{2}\left(t+t_{i}\right)}\right). \end{array} $

(14)

式中：T_wi和T_i分别为第i挡时水温和混凝土温度，ϕ_i(t)、ψ_i(t)和p_i分别为第i挡时水冷函数、水冷温升函数和水冷参数，t_i为流量或水温改变时刻。当流量或水温改变时t必须从0开始。

水冷参数公式为

$ p=g k a / D^{2}, $

(15)

$ k=2.09-1.35 \xi+0.320 \xi^{2}, $

(16)

$ \xi=\frac{\lambda L}{c_{w} \rho_{w} q_{w}}. $

(17)

对于聚乙烯水管

$ g=\frac{a \ln 100}{\ln (b / c)+\left(\lambda / \lambda_{1}\right) \ln \left(c / r_{0}\right)}, $

(18)

$ b=0.583\ 6 \sqrt{S_{1} S_{2}}. $

(19)

式中：a和λ分别为混凝土的导温系数和导热系数，L为冷却水管长度，c_w、ρ_w和q_w分别为冷却水比热容、密度和流量，D、b和c分别为等效均质混凝土圆柱体直径、外半径和内半径，r₀和λ₁分别为水管的内半径和导热系数，S₁和S₂分别为冷却水管水平间距和铅直间距。

2.2 冷却水沿途水温计算

文[30]给出了冷却水沿途温升的3种推算公式，本文选取迭代算法。如图 3所示，假定在i和i+1截面之间，沿水管外缘C₀的积分$ \int_{{C_0}} {\frac{{\partial T}}{{\partial r}}} {\rm{d}}s $沿水管长度呈线性变化，于是ΔL区间内水温的增量为

$ \Delta T_{\mathrm{w} i}=\frac{\lambda \Delta L}{2 c_{\mathrm{w}} \rho_{\mathrm{w}} q_{\mathrm{w}}}\left[\left(\int_{C_{0}} \frac{\partial T}{\partial r} \mathrm{~d} s\right)_{i}+\left(\int_{C_{0}} \frac{\partial T}{\partial r} \mathrm{~d} s\right)_{i+1}\right], $

(20)

$ T_{\mathrm{w}i}=T_{\mathrm{wl}}+\sum\limits_{j=1}^{i-1} \Delta T_{\mathrm{w}{j}} \quad i=2,3,4, \cdots. $

(21)

重复迭代计算，直至各截面的水温趋于稳定，迭代结束控制指标为

$ \mathop {\max }\limits_i \left| {\frac{{T_{{\rm{w}}i}^k - T_{{\rm{w}}i}^{k + 1}}}{{T_{{\rm{w}}i}^{k + 1}}}} \right| \le \varepsilon . $

(22)

2.3 优化模型求解方法

混凝土热学参数主要包括密度ρ、比热容c、导热系数λ、导温系数a、绝热温升θ₀和温升规律n，其中室内试验获得的密度值、比热容值、导热系数值和导温系数值，其受环境影响变化很小。而绝热温升θ₀和温升规律n作为关键控制因素，直接影响最高温度和温度变化过程，有必要进行反演分析计算。混凝土浇筑后低温季节需表面覆盖保温被，直接影响等效表面放热系数值β。因此，本文选取θ₀、n、β作为反演参数。

目标函数为

$ F(\boldsymbol{X})=\sum\limits_{j=1}^{s} \sum\limits_{i=1}^{t}\left(\frac{T_{i j}^{\mathrm{C}}-T_{i j}^{\mathrm{M}}}{T_{i j}^{\mathrm{M}}}\right)^{2} \rightarrow \min. $

(23)

式中：X=[x₁, x₂, x₃]^T=[θ₀, n, β]^T，s为测点数，t为时长，T_ij^C和T_ij^M分别为温度计算值和监测值。

3 工程算例分析

我国西南地区在建的白鹤滩大坝^[31-32]，混凝土标号有C40、C35和C30。为了跟踪低温季节C40和C35混凝土、C40不同级配混凝土发热规律，考虑环境气温变化以及实际通水冷却方案，对混凝土热学参数进行优化反演。

3.1 C40和C35混凝土反演

如图 4所示，建立19#坝段013~017的5个浇筑仓有限元模型，每个浇筑仓厚3 m，其中013~015混凝土标号为C40，016~017混凝土标号为C35。混凝土基本热学参数为：C40混凝土密度ρ为2 601 kg/m³，比热容c为0.851 kJ/(kg·℃)，导温系数a为0.070 8 m²/d；C35混凝土ρ为2 553 kg/m³，比热容c为0.873 kJ/(kg·℃)，导温系数a为0.073 m²/d。通过六面体八节点单元划分了50 000个单元，54 621个节点。

浇筑仓分6个坯层，每个坯层厚0.5 m，在第3坯层底面以上0.15 m处布置分布式光纤，埋设位置如图 5所示。取光纤平均温度作为实测值，19#013浇筑仓光纤温度不参与反演。第1和第4坯层均埋设1层冷却水管，水管的内外半径分别是14和16 mm。19#坝段013~015浇筑仓水管水平间距1 m，铅直间距1.5 m。19#坝段016~017浇筑仓水管水平间距和铅直间距均为1.5 m。

依据浇筑进度取实际仓间间歇时间，选用智能通水系统来选取冷却通水时间、水温和流量，环境气温取0.5 d平均气温。时间步长取为0.5 d，总时间为58 d。选用30个粒子，参数反演结果如表 2所示。

图 6给出了通过混合粒子群算法和传统粒子群算法计算得到的适应度随迭代次数变化曲线，随着迭代次数的增加，收敛速度明显加快，混合粒子群算法更优。

C40和C35混凝土绝热温升和温升规律收敛过程如图 7所示，混凝土热学参数反演效果较好，20次迭代后目标函数值已稳定，收敛速度平稳。19#坝段014~017浇筑仓温度计算结果与实测结果对比如图 8所示，变化趋势一致，反演获得结果能很好反映真实浇筑环境下的混凝土热学参数。

3.2 C40不同级配混凝土反演

建立19#坝段012~014的3个浇筑仓有限元模型如图 9所示，每个浇筑仓厚3 m，共有28 800个单元，31 899个节点，浇筑仓有二级配、三级配以及四级配混凝土。二级配混凝土密度ρ为2 532 kg/m³，比热容c为0.885 kJ/(kg·℃)，导温系数a为0.075 6 m²/d；四级配混凝土密度ρ为2 674 kg/m³，比热容c为0.863 kJ/(kg·℃)，导温系数a为0.072 48 m²/d。光纤埋设方案如图 10所示，选取点B作为二级配混凝土中心点温度，顺河向光纤平均值作为四级配混凝土中心点温度。利用迭代算法获得各个不同级配混凝土进口水温，对二级配混凝土绝热温升和温升规律、四级配混凝土绝热温升和温升规律、等效表面放热系数这5个参数进行反演，选取30个粒子，计算结果如表 3所示。并据此进行19#坝段013~014浇筑仓二级配及四级配混凝土温度计算，实测结果和计算结果如图 11和图 12所示。对于温度偏差，二级配比四级配混凝土更大，主要是二级配混凝土典型监测点并未完全位于二级配混凝土中心点部位，冷却通水流量不平稳。

4 结论

(1) 特高拱坝热学参数反演实际上是多参数组合的大空间搜索问题，本文建立了混合粒子群模型，将凹函数权值递减系数引入到粒子群算法中，并与遗传算法交叉变异策略结合。通过算法性能测试，混合粒子群模型更能逼近全局最优解。

(2) 考虑多挡通水以及冷却水沿途水温变化影响，依据已埋仓光纤测温数据，将混合粒子群模型用于白鹤滩拱坝低温季节不同材料热学参数识别。工程算例表明，混合粒子群模型可以提高反演的参数精度。