混凝土被广泛应用于工程建设，其内部缺陷的准确预测与预报，对于确保结构物的正常工作和避免建筑物的灾难性损伤具有重大意义。目前，混凝土缺陷检测技术主要有超声脉冲法^[1]、雷达扫描法^[2]、红外热谱法^[3]等。这些方法应用广泛，但均存在一定的局限性。如混凝土对超声的强吸收问题影响了超声脉冲技术在混凝土检测中的应用。雷达法检测仪器多采用国外进口，成本较高，且无法对原始数据进行分析，该技术在国内尚处于起步阶段。红外热谱法应用较多，但其主要应用于定性判断和表面质量检测，难以实现对缺陷定量及目标内部的可视化检测。

近年来，电阻率成像技术(ERT, electrical resistance tomography)因其成本低、无损伤及功能成像等优点^[4]被广泛用于工程检测领域^[5]。ERT技术是一种新型的图像重建技术。ERT技术是通过向被测对象注入电流，测量边界电压(或电势)，经过反演计算，重构被测对象内部的电导率，进而显示物体内部状态的一种成像技术。但是，ERT是一个病态反问题^[6]，存在着非线性、欠定性、严重不适定性和计算量大等难点，成像的精度和稳定性很难保证。为提高图像重建质量，文献[7]提出阈值修正法，可去除部分图像伪影，改善成像质量，但选取阈值大小时，主观因素影响较大，难以寻找到合理阈值。董峰等^[8]在阈值修正法基础上提出一种快速自适应稀疏重建算法，研究了噪声对图像重建效果的影响，获得了较好效果。肖理庆^[9]提出一种改进牛顿拉夫逊绝对图像重建算法，采用灵敏度矩阵自动更新的策略，结合图像结果自动选取阈值修正电阻率，可有效降低图像的相对误差。Cheney等^[10]提出牛顿单步残差法(NOSER, Newton’s one-step error reconstructor)，设定均匀分布的电导率为初值，用一次牛顿迭代后的结果进行成像，可实现快速成像，但其成像伪影较大，电阻率与真实值相差较大，仅能定性分析，无法实现定量分析。

将ERT技术应用于混凝土检测领域时，由于混凝土自身电导率较小，测量受外界因素影响较大，且混凝土基体与缺陷间的电导率为突变而不是渐变关系等原因，用于医学检测的传统ERT算法很难将缺陷与基体间的锐利边界重构出来，限制了ERT技术在混凝土无损检测中的应用。

采用ERT技术进行混凝土缺陷检测时，被测对象有两个特点：基体与缺陷有明显边界(即电导率发生突变而不是渐变)；混凝土的电导率可事先测量。针对这些特点和上述算法的问题，文中提出一种优化算法——约束牛顿迭代算法(NCRM, Newton’s constrained reconstruction method)，依据被测对象的这些信息作为先验条件，结合上一次迭代结果来调整初始向量，逐次递进反算初值，通过进一步迭代来得到较为准确的内部电导率分布，实现混凝土缺陷无损检测目的。

1 ERT工作原理

ERT技术实质上是低频电流场的反问题计算，主要包括正问题和反问题的研究。正问题为已知场域电阻率分布，求解边界电极处电压值，计算过程稳定，有唯一解；而其反问题是已知边界电压，通过反运算推测场域内各单元的电阻率分布，一般情况下能获得的边界电压值有限，而单元数远多于电压数量，计算呈现明显的病态性和欠定性，数学上的解决方案多采用Newton-Raphson迭代和正则化方法来降低病态性，用最小二乘法获得近似解。

从物体表面注入低频恒定电流，忽略介电常数，将场域内的电流场认为是准静态场，则电阻率σ与电势分布函数φ满足Laplace方程^[11]，即

$ \mathit{\Omega }:\nabla \cdot (\sigma \nabla \varphi ) = 0, $

(1)

$ {\mathit{\Gamma }_1}:\mathit{\boldsymbol{\varphi }} + {z_l}\sigma \frac{{\partial \varphi }}{{\partial \mathit{\boldsymbol{\bar n}}}} = {U_l}, $

(2)

$ {\mathit{\Gamma }_2}:\int \sigma \frac{{\partial \varphi }}{{\partial \mathit{\boldsymbol{\bar{n}}}}}d\mathit{\Gamma } = - {I_l}, $

(3)

式中：φ为场域Ω内的电势分布函数；z_l代表第l个测量电极的接触阻抗；Ι_l是电极处的电流密度；n为外法向单位矢量；σ为场域内的电导率分布函数。

正问题是根据边界条件与场域内电导率分布σ求解电势分布φ。采用有限元法求解时，可写成矩阵形式

$ \mathit{\boldsymbol{S\varphi }} = \mathit{\boldsymbol{F}}, $

(4)

式中：S为系数矩阵，是由单元系数矩阵S_ij^e组装而成，里面包含单元的电导率，此处假设每个单元的电导率为常数，S_ij^e可由单元形函数推导得到；F为电流列向量，电流注入流出的位置分别为I和-I，其余项均为0。

ERT的反问题是由边界测量电压计算被测物内部的电导率分布，一般采用最小二乘法搜索出一合适的电导率分布σ，使在此电导率分布的计算电压U_i(σ)与测量电压V_i足够接近，即(V_i-U_i(σ))(i=1, 2, …, N)达到最小，即：

$ E(\sigma)=\sum\limits_{i=1}^{N}| | V_{i}-U_{i}(\sigma)| |^{2}=\sum\limits_{i=1}^{N} \sum\limits_{j=1}^{N}\left(V_{i j}-U_{i j}(\sigma)\right)^{2}, $

(5)

式中：E(σ)是电导率分布σ的函数；v_ij是第i次激励第j次测量得到的测量电压；U_ij(σ)是在电导率分布为σ时由ERT正运算得到第i次激励第j次计算电压；N是电极数量。

寻找电导率分布σ时，使E(σ)达到最小，可令其偏导数为零，即

$ {F_m}\left( \sigma \right) = \frac{{\partial E\left( \sigma \right)}}{{\partial {\sigma _m}}} = - 2\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {\left( {{V_{ij}} - {U_{ij}}\left( \sigma \right)} \right)} } \frac{{\partial {U_{ij}}\left( \sigma \right)}}{{\partial {\sigma _m}}} = 0,\left( {m = 1,2, \cdots ,M} \right)。$

(6)

根据牛顿迭代法，有

$ \sigma^{(k+1)}=\sigma^{(k)}-\left[\boldsymbol{J}\left(\sigma^{(k)}\right)\right]^{-1} \cdot F\left(\sigma^{(k)}\right),(k=0,1, \cdots), $

(7)

式中：k为迭代次数，一直到满足精度要求为止；J(σ^k)为雅克比矩阵，计算公式为

$ {J_{n,m}}\left( {{\sigma ^{\left( k \right)}}} \right) = \frac{{\partial {F_N}\left( {{\sigma ^{\left( k \right)}}} \right)}}{{\partial {\sigma _m}}},\left( {n,m = 1,2, \cdots ,M} \right)。$

(8)

以上是ERT正问题和反问题的数学表示，反问题求解的关键是尽量找到较好的、与真实电导率分布相近的初值σ^(k)。

2 NCRM算法图像重构原理

当采用牛顿类算法计算时，如果初始电导率分布接近真实电导率分布，那么迭代次数就会减少很多。NOSER方法^[12]用整个区域的平均电导率作为迭代的初值，仅通过一次迭代来实现ERT成像，结果虽不精确，但能从重建电图像中得到场域内目标的位置和大致尺寸，可实现快速成像。采用的NCRM算法在此基础上，利用已知的先验条件来多次调整初值，以获得更为准确的结果，由于被测物体内电导率分布未知，初值仍是假设场域内电阻率均匀分布，计算出第一次的初值X₁，即X₁=c(1, 1, …, 1)^T, 设l⁽⁰⁾=(1, 1, …, 1)^T，向量l⁽⁰⁾为所有元素等于1的单位向量，计第k次迭代l^(k)。c的大小利用最小二乘法确定，代入E(σ)中求得

$ E\left( {c{\mathit{\boldsymbol{l}}^{(k)}}} \right) = \sum\limits_{i = 1}^N | |{V_i} - {U_i}\left( {c{\mathit{\boldsymbol{l}}^{(k)}}} \right)|{|^2} = \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {{{\left( {{V_{ij}} - {U_{ij}}\left( {c{\mathit{\boldsymbol{l}}^{(k)}}} \right)} \right)}^2}} } , $

(9)

式中:c为初始电导率分布中基体的电导率; l^(k)为第k次迭代时对单位向量l⁽⁰⁾的修正向量；V_i为相邻电极上的测量电压值；U_i(cl^(k))为利用有限元方法计算的相邻电极电压值; c为常数，可转换为U_ij(cl^(k))=cU_ij(l^(k))，求得c值为

$ c = \sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {\left( {{V_{ij}}{U_{ij}}\left( {{\mathit{\boldsymbol{l}}^{(k)}}} \right)} \right)} } /\sum\limits_{i = 1}^N {\sum\limits_{j = 1}^N {{{\left( {{U_{ij}}\left( {{\mathit{\boldsymbol{l}}^{(k)}}} \right)} \right)}^2}} } ,{(i,j = 1,2, \cdots ,N)}。$

(10)

利用式(10)求解系数c，将X₁=c(1, 1, …, 1)^T作为初值带入方程(7)求得近似电导率分布Y₁=(σ₁⁽¹⁾, σ₂⁽¹⁾, …, σ_m⁽¹⁾)^T，同时也获得第一次重构图像，由结果可以确定最有可能是缺陷的位置(单元)，但伪影较多，原因是均匀分布的电导率与实际电导率分布相差较大，没有很好地满足牛顿法对初值的要求。从实验结果看，缺陷的识别范围也明显大于实际尺寸，缺陷的中心位置可由此判定，试验证明中心位置与实际缺陷中心位置也基本一致。以该近似电导率分布Y₁为基础，NCRM算法将混凝土基体及缺陷的电导率(σ_c、σ_f)作为先验信息，依据第一次重构图像确定的缺陷中心位置，将混凝土电导率σ_c和缺陷电导率σ_f作为先验条件代入算法中，调整分布向量l⁽⁰⁾为l⁽¹⁾，即由第一次的均匀分布l⁽⁰⁾=(1, 1, …, 1, …, 1)^T调整为l^(k)=(1, 1, …, $ \frac{{{\sigma }_{f}}}{{{\sigma }_{c}}} $, …, 1)^T，l^(k)中$ \frac{{{\sigma }_{f}}}{{{\sigma }_{c}}} $为缺陷电导率与基体电导率的比值，在向量l^(k)中的位置(序号信息)就是前一次分析结果中的缺陷中心位置对应的单元编号。将引入先验信息的向量l^(k)代入公式(10)中，式中v_ij为测量电压值，由试验获得；U_ij(l^(k))为假设电导率分布为(1, 1, …, $ \frac{{{\sigma }_{f}}}{{{\sigma }_{c}}} $, …, 1)^T时的边界电压值，可通过正运算分析得到，由公式(10)可得到优化后的被测对象电导率分布X_k=σ(k)_cl^(k)，显然，X_k考虑了第一次成像的结果，也引入了先验信息，相比均匀分布X₁更接近于真实电导率分布，将此电导率分布作为初值代入公式(6)，计算后可得到更接近于真实值的电导率分布Y_k，按照此方法对第k次初始电导率分布X_k施加约束，以此类推，至所有缺陷都被定位后停止迭代。

需要指出的是，如果测量电压不含噪音，是理想的数据源时，由于有严格的数学推演，采用文中方法可基本消除伪影，获得数学上最优, 并与实际相符的反运算结果。该算法原理也适用于多目标成像，如试件中含有多种缺陷(如孔洞、裂缝或其他异物等)，σ_f即该缺陷的测量电导率值，将$ \frac{{{\sigma }_{f}}}{{{\sigma }_{c}}} $代入要调整的先验向量l^(k)，即可获得满意的成像结果。此外，作为先验信息的混凝土和缺陷的电导率值并不要求精确，只需取其近似的测量值，能够反映出两者之间的差异程度即可获得较好的成像结果，因为在该算法中，其电导率值被用来调整迭代过程中的初始值，即X，且每次调整均需满足公式(9)以及迭代公式(7)，所以，调整后的初值均比均匀分布的预估初值更加接近真实解。

为对比施加约束前后重构图像性能指标，文中拟采用相对误差及相关系数进行量化评价。研究中混凝土基体电导率以实测值作为先验条件，所以，在图像指标评定中，试件真实电导率可假设为缺陷位置电导率为0 ms/cm，其余部分为均匀基体，电导率为实测值。因实验中试件采用的是细集料，均匀性较高，所以假设的真实电导率分布具有参考价值。

相对误差：

$ R_{\mathrm{E}}=\frac{\left\|x^{*}-x\right\|}{\left\|x^{*}\right\|_{2}}, $

(11)

相关系数：

$ C_{\mathrm{C}}=\frac{\sum\nolimits_{i}^{N}\left(x_{i}-\overline{x}\right)\left(x_{i}^{*}-\overline{x}^{*}\right)}{\left\|x_{i}-\overline{x}\right\|_{2}\left\|x_{i}^{*}-\overline{x}^{*}\right\|_{2}}, $

(12)

式中：x为电导率分布测试结果；x^*为与测试结果对应的真实电导率分布(假设值)；“-”符号表示平均值。相对误差越小、相关系数越高则成像质量越好。

3 实验材料与方法

实际混凝土检测中，由于模型误差、接触阻抗、测量误差是不可避免的，输入的测量电压一定含有噪音数据。为验证NCRM约束算法的可行性和实用性，文中采用预设孔洞的圆盘形水泥试件模拟混凝土缺陷。在重构算法迭代过程中，利用缺陷与混凝土电导率的比值作为先验条件，该缺陷为孔洞，其电导率趋近于0，计算时直接采用σ_f=0作为先验条件干预反运算的迭代过程，即可获得较为理想的重构图像。

3.1 试样制作

水泥砂浆试件的材料为P42.5普通硅酸盐水泥和标准砂(粒径范围0.25~0.65 mm, 密度2.59 g/cm³)，其中，水灰比为0.5。实验中共有2组试件，如图 1所示。第1组试件为圆盘形试件，直径10 cm，高2 cm，分别预设1个和2个直径为9 mm的圆孔，用于模拟混凝土缺陷的电阻成像检测实验。第2组试件为方块形试件，尺寸为5 cm×5 cm×2 cm，用于检测混凝土电导率，为电阻成像检测试验提供先验信息。2组试件采用同种配比同时浇筑，在相同条件下养护。

3.2 成像硬件系统简介

电阻层析成像系统主要由电极组件、电流源及数据采集系统、数据处理系统组成。文中采用自行搭建的电阻成像系统，如图 2所示。硬件部分包括：①计算机(包含数据处理和成像算法程序)；②可编程控制转换开关；③电流源(美国吉时利Model 6221，能够为各种元件和材料的交流特征分析提供2 pA~100 mA的交流电流源)；④电压采集仪(安捷伦34 970 A)；⑤电极夹具。软件部分是基于Matlab环境下自编译的NCRM约束算法程序。工作原理为注入激励电流，用可编程控制转换开关实现相邻激励模式^[13]，计算机通过接口将指令传达给控制开关和电压采集仪，控制开关执行指令从16个电极中选择激励电极对，同时电压采集仪测量其他电极对电压并将电压测量值返回计算机作为图像重建的原始数据。

因混凝土表面粗糙且坚硬，为保证电压测量精度，在电极接触方面尤其重要。考虑混凝土表面特性，自制电极夹具可有效减小电极之间的接触阻抗差。每个电极端部用导电泡沫(导电布包裹泡沫)，导向杆采用激光钻孔精确定位(见图 2)，通过弹簧施加压力，使导电泡沫与混凝土试件表面接触良好并受力均匀。

4 实验与结果分析 4.1 试件电阻检测试验

NCRM约束算法利用了缺陷电导率及混凝土电导率范围的先验信息，试验中的缺陷为孔洞，室内环境下空气电导率接近零，在计算中缺陷电导率选取为σ_f=0 ms/cm。混凝土电导率采用LCR精密电桥(Applent AT2818)测得，如图 3所示。

电阻层析成像中交流电频率需尽可能低，以避免由于连接到测量单元的电缆过长而产生的电极电感效应，另外，频率必须足够高，以避免电极极化效应。对于电阻成像来说，当相位角最小时，电阻抗虚部最小，对应的交流电频率最为合适。将试件从标准养护箱中取出置于实验室12 h，实验室环境为23 ℃，相对湿度70%，此时，试件中含水率较高，电阻较干燥状态小，有利于验证NCRM约束算法的可行性。利用LCR仪器检测试件电阻抗时采用1~300 kHz的交流电压，检测结果如表 1所示，根据相位角的大小，采用激励电流频率为50 kHz。实验测得电阻抗为4.082 kΩ，相位角8.48°，电阻为4.037 kΩ，根据试件尺寸计算其电导率为σ_c=0.123 ms/cm。

4.2 图像重建结果与分析

与电阻抗测量试件一致，试件从标准养护箱中拿出来后置于室内环境下12 h后开始试验，室内湿度70%，温度23 ℃，实验中电极布置方式如图 4(a)所示。重构成像在Matlab平台上进行，有限元模型如图 4所示。单元数量256个，节点数145，16个电极均匀分布在试件边缘。

实验采用直径为10 cm的混凝土试件，设置了1个直径为9 mm的空心孔来模拟单个缺陷(见图 4(a))，设置2个直径为9 mm的空心孔来模拟多个缺陷(见图 4(b))，将前面测得的混凝土试件电导率(0.123 ms/cm)作为先验信息，得到缺陷与混凝土机体的电导率比值$ \frac{{{\sigma }_{f}}}{{{\sigma }_{c}}} $。实验采用电流相邻激励模式，电流源提供10 mA，50 kHz的正弦交流电源，通过转换开关在16电极上进行切换，每次切换过程中通过数据采集模块得到16个电压数据，整个数据采集可获得16×16=256个电压数据，扣除无效电压数据(即在每次电流注入电极位置采集的3个电压含有较大的接触阻抗，误差较大，视为无效电压不进行运算)，由于对称性，有效电压值为16×13/2=104个，其中，最大电压值为8.05 V，最小电压值为0.06 V。将测量结果代入初值计算式(9)，得到第1个试探值，σ_c⁽¹⁾=0.116 ms/cm，与测量得到的电导率0.123 ms/cm相差较小，第1次迭代假设试件的电阻率全范围均匀分布，第1次迭代的初值采用实际测量值0.123 ms/cm，当然也可以采用约束公式的计算值。通过第1次计算后的结果Y₁可以获得重建图像(见图 5(a)和图 5(b)所示)，图像中颜色最深的位置与试件圆孔位置大致相同，按照NCRM算法，选取Y₁结果的最小值以及所对应的单位编号，将已知的约束条件将计算结果Y₁代入公式(10)进行重新计算，这里的已知约束条件包括：混凝土的电导率、孔洞的电导率(假设为0)以及孔洞和混凝土的锐利边界关系。得到X₂(1, 1, 1, 0，…，1, 1, 1)，这里“0”所在的位置即为Y₁结果最小值所在位置(即对应单元编号)。将X₂代入程序计算得到Y₂，如图 5(c)和图 5(d)所示。根据Y₂的成像结果可以看出伪影明显减少，缺陷的位置也更加明确，计算结果Y₂中仍然有比基体数值明显大的，但考虑到实验试件中缺陷实际大小与划分的有限元模型单元不能完全匹配重合，所以，实验中只进行了2步迭代。若需进一步约束迭代计算，此时再设置X₃的列向量(1, 1, 1, 0，…, 0，…, 1, 1, 1)，即第1次试探得到的单元位置保持不变，第2次找到的单位位置赋值为0，得到X₃，计算结果Y₃。

图 5为NCRM算法进行试样缺陷检测的2次计算结果Y₁和Y₂的成像图和电导率散点图。图 5(a)和图 5(b)为Y₁的结果，图 5(c)和图 5(d)为Y₂的结果。当采用均匀电导率分布为初值计算时，Y₁成像图有较多的伪影，但可以准确定位缺陷的中心位置，由对应的图 5(b)针状图可以看出，反运算后电导率的分布呈阶梯状逐渐变化，虽可以定位缺陷位置，但无法分辨缺陷与基体的边界，也就无法判断缺陷的大小程度、分布范围；当对初值X₁进行约束一次后，计算结果Y₂成像图明显减少了伪影，显示的缺陷位置集中分布在实际孔洞处，对应针状图(图 5(d))可以看出，基体与缺陷的电导率呈大致的二值性，梯度渐变的趋势减弱，反映在成像图(见图 5(c))上即出现了更符合真实情况的较为锐利的边界过渡。对于单目标缺陷检测，仅迭代2次即可获得较为满意的结果。

图 6为NCRM算法进行试样缺陷检测的4次计算结果Y₁，Y₂，Y₃及Y₄的成像图和电导率散点图。图 6(a)、图 6(c)、图 6(e)、图 6(g)为电导率重构图像，图 6(b)、图 6(d)、图 6(f)、图 6(h)为电导率二维散点图。当采用均匀电导率分布为初值计算时，Y₁成像图有较多的伪影，如图 6(a)所示，但可以准确定位缺陷的中心位置，由对应的图 6(a)针状图可以看出，逆运算后电导率的分布呈阶梯状逐渐变化，虽可以定位缺陷位置，但无法分辨缺陷与基体的边界，也就无法判断缺陷的大小程度、分布范围；对初值X依次进行3次约束后，可以看到，计算结果Y₂，Y₃及Y₄成像图的伪影成递减状态，显示的缺陷位置集中分布在实际孔洞处，虽然重构得到的缺陷中心与实际有一定的偏差，但是相对于不加约束前的重构图像，成像质量有了明显的改善。同样，从对应图 6(d)、图 6(f)、图 6(h)可以看出，基体与缺陷的电导率呈大致的二值性，梯度渐变的趋势减弱，反映在成像图上即出现了更符合真实情况的较为锐利的边界过渡。对于2个目标缺陷检测，经过3次迭代可获得较为满意的结果。

4.3 图像评价指标

按照式(11)和式(12)分别计算2组实验每次迭代计算后的相对误差和相关系数，结果如表 2和表 3所示。

由表 2和表 3可以看出，施加约束后，电导率分布重构图像的相对误差明显减小，单孔缺陷施加1次约束就得到了较好的成像结果。两孔缺陷在施加3次约束后，相对误差的减小较为显著，图像相关系数也明显提高，说明NCRM约束方法可以用于多目标缺陷的检测。通过实验, 验证了文中提出的NCRM约束算法相对传统算法能够更好地进行混凝土试件缺陷检测。计算结果仍含有部分伪影，主要原因有：

1) 原始测量数据无法完全消除噪音，因为测量误差及接触阻抗在各电极间的差异客观存在；

2) 数值建模采用三角形单元，单元划分是任意的，不能事先考虑缺陷的位置和范围，并且有限元分析中假设同一单元内的电导率是常数，不可避免地出现某些单元同时跨越了实际被测对象的基体和缺陷，而这些单元的电导率应该同时包含两者信息，这些系统误差也会影响成像精度；

3) 逆运算并不是单元越密效果越好，对于16电极系统，独立有效的测量数据为104个，通过这些独立测量数据反演模型中各个单元的电导率本身就是一个欠定问题，单元数越大，问题欠定性越强，结果越不可控。

因此，NCRM算法增加的约束条件即为降低反问题病态性，相当于增加了一些有效数据给逆运算系统。

5 结论

对于混凝土检测来说，由于ERT传统算法偏重于解决稳定性问题，不能很好地反映机体与缺陷之间锐利的边界，使得重建图像与真实情况有较大差别。NCRM算法利用先验信息对ERT反问题进行优化，重建结果表明该方法是有效的，可以降低伪影，提高目标位置的定位精度，迭代次数少，计算速度快。

该算法增加了约束条件来干预反演过程，需要事先获知被测对象及缺陷的电导率范围，且只适用于缺陷与基体间有较明显的电导率突变现象的情况。因此，文中NCRM算法用于混凝土孔洞、裂缝的ERT检测效果较好。对于混凝土渗水检测，因其内部电导率在干湿界面处是渐变的，与约束条件不符，故不能用此算法求解。