传统的复合材料损伤检测方法如超声、涡流和红外成像等在检测时间和应用场景上有一定的局限性[4-7]。根据碳纤维自传感特点及结构损伤电学敏感特性提出的电学阻抗检测方法，具有非侵入、无辐射、响应快等优点，在碳纤维复合材料无损检测领域已逐步受到国内外****的广泛关注[8]。而基于电阻抗检测原理的电阻抗层析成像(Electrical Impedance Tomography, EIT)方法[9]，近年来也被应用于碳纤维增强复合材料无损检测。
在电阻抗层析成像过程中，电阻抗场域内边界电压数据远小于网格刨分的数量，这使得EIT逆问题表现出严重的欠定性和病态性[10]。目前，一般采用正则化的方法来提高EIT逆问题的稳定性和抗噪能力。Clay和Ferree[11]通过奇异值分解正则化并对不同深度的灵敏度矩阵进行归一化处理来研究成像逆问题。Gonzalez等[12]将全变差正则化运用于EIT逆问题中，提高了图像中目标边界的分辨率。李星等[13]将对角权重正则化应用于EIT逆问题，有效提高了EIT的抗噪能力和重建图像的质量。范文茹等[14]采用稀疏正则化的方法，减少了EIT逆问题的重建图像的伪影，提高了EIT逆问题的重建图像损伤边缘分辨率。
本文利用CFRP层压板发生结构损伤导致其电导率值单向变化的特点，结合逆问题求解的非负特性，研究基于修正残差范数最陡下降法(Modified Residual Norm Steepest Descent, MRNSD)，以改善重建图像质量。同时，为进一步提高成像精度、抑制噪声以及解决迭代算法的半收敛现象，结合灵敏度映射的预处理方法和软闭值方法，提出改进的MRNSD算法，并通过仿真和实验对该算法进行验证。
1 CFRP层压板特性 CFRP是由碳纤维和树脂聚合物组成的复合材料，而树脂聚合物是绝缘材料，碳纤维是导电材料，且碳纤维的纤维方向的电导率较高，纤维方向的垂直方向和层压板的厚度方向的电导率较低，因此CFRP的导电性具有各向异性的特点[15]。
通常将CFRP简化为单层均质连续各向异性材料，定义电导率张量为
 | (1) |
式中: σ11为纤维方向电导率；σ22为纤维方向的垂直方向电导率；σ33为层压板厚度方向的电导率。在本文中使用COMSOL软件构建八层CFRP层压板正交([0°/90°]4)模型，其中每层碳纤维(长10 cm、宽10 cm、厚0.04 cm)为各向异性材料。相邻层为不同纤维方向，在COMSOL软件中将0°纤维电导率设置为
 |
| 图 1 正交型CFRP层压板模型 Fig. 1 Orthogonal CFRP laminate model |
| 图选项 |
2 电阻抗层析成像正问题 EIT图像重建是通过边界电压V和输入电流I得到电导率分布σ′。通过建模和有限元离散化就得到了EIT问题的确定性观测模型为
 | (2) |
式中：U(σ′; I)为σ′和I到V映射的正演模型；R(σ′)为σ′到电阻的映射模型。V=R(σ′)I模型与σ′呈非线性关系，与I呈线性关系。在电导率变化很小的情况下，可以考虑通过线性化方程组的方法来精确地求解EIT反问题。
 | (3) |
式中：δσ∈Rn×1为电导率张量模值变化量的分布矩阵，n为重建图像中的像素数；δUv∈Rm×1为材料损伤前后电压测量差值，m为测量值个数；σ0为材料初始电导率分布；J∈Rm×n为Jacobian矩阵。
根据四端口网络的Geselowitz灵敏度定理，灵敏度图的快速计算方法如下：
 | (4) |
式中：u(Id)和u(If)分别为第d次和第f次驱动模式的电势分布；Jdfe为第d次和第f次驱动的灵敏度矩阵；Vdf为边界测量电压；σe为电导率分布矩阵；Ωe为被测区域；dV为对体积的积分。
3 电阻抗层析成像逆问题 3.1 MSRND算法 根据式(3)，EIT采用差分成像方法。当CFRP层压板发生冲击、分层、裂纹等损伤时，即会出现碳纤维断裂，进而导致材料电导率减小。因此，差值成像中的δσ可认为是单向变化量。根据先验信息，保证δU是损伤前与损伤后的电压差，则式(3)中δσ具有非负特性。利用该特点，结合MSRND算法对CFRP层压板的EIT图像重构进行深入研究。
常用的求解EIT逆问题算法，比如共轭梯度算法等无约束Krylov子空间方法，但常用方法没能考虑到解的非负性[16]，导致成像精度较低。Nagy和Strakos[17]提出了一种非负约束极小化算法，设极小化函数为
 | (5) |
式中：σreg为目标电导率分布。
强制非负约束性条件为δσ≥0。令δσ=ez，通过参数化公式(5)可得
 | (6) |
式中：G(z)为变换后的最小化函数，z为z域。
根据链式求导法则，有
 | (7) |
令gradzG(z)=0, 可以得到该约束极小化问题在δσ空间上的非线性规划最佳解的必要条件(Karush-Kuhn-Tucker, KKT)，得到以下迭代格式：
 | (8) |
式中：pk=FkJT(δUv-Jδσk)；αk为搜索参数。通过计算极小化残差范数
 | (9) |
进一步假设sk=JTrk和pk=Fksk，则可进一步推出
 | (10) |
因此，第k+1次迭代的残差向量为
 | (11) |
满足
算法1??MRNSD算法。
输入: Jacobian矩阵J, 迭代初始向量δσ0和迭代终止条件。
输出:重建电导率δσ。
初始化:δσ=δσ0; s=JT(δUv－Jδσ); F=diag(δσ); γ=sTFs。
如果k=1, 2, …，则迭代:
p=-Fs
u=Jp
δσ=δσ+αp
F=diag(δσ)
y=JTu
s=s－αy
γ=sTFs
如满足终止条件，则迭代停止，否则继续迭代。
其中：s为中间变量；α为搜索参数。
3.2 预处理改进 预处理经常被用来加速迭代算法的收敛以及提高计算精度。在线性化算法中，由于EIT的非线性和不适定性，灵敏度矩阵直接影响重建图像的质量。在灵敏度矩阵的计算中采用一些先验信息，可以使场域灵敏度更加均匀，改善求解的不适定性。灵敏度的归一化方法在电容层析成像(Electrical Capacitance Tomography, ECT)中得到了广泛的应用[18]。在此基础上，提出了一种改进的基于灵敏度映射的预处理器。
灵敏度矩阵J在测量区域边界附近的灵敏度最大，这导致了物理项的不定性，导致重建图像的误差更大。因此，使用平滑预条件校正方案，即
 | (12) |
令C=diag(ω1, ω2, …, ωk)且ωk≠0，来保证矩阵C是非奇异，而ωk的有效求法是通过如下求和操作，从而将每个像素的总灵敏度标准化，
 | (13) |
式(5)中极小化问题的预处理优化则可描述为：在每次迭代中使灵敏度矩阵J右乘非奇异矩阵C，边界测量电压变化值向量δσ左乘非奇异矩阵的逆矩阵C-1。然后将迭代法应用于预处理系统，就可以得到算法2，即预处理修正残差范数最陡下降法(Pretreatment Modified Residual Norm Steepest Descent，PMRNSD)。
算法2??PMRNSD算法。
输入: Jacobian矩阵J, 迭代初始向量δσ0和迭代终止条件。
输出:重建电导率δσ。
初始化:δσ=δσ0; s=JT(δUv－Jδσ); F=diag(δσ); γ=sTFs。
如果k=1, 2, …，则迭代:
p=-Fs
u=C-1Jp
δσ=δσ+αp
F=diag(δσ)
y=JTC-Tu
s=s－αy
γ=sTFs
如满足终止条件，则迭代停止，否则继续迭代。
3.3 软闭值改进 使用共轭梯度等无约束迭代正则化算法求解EIT反问题会出现半收敛现象，MRNSD算法相对较为稳定，但因噪声的影响，仍然会存在半收敛现象。即重建图像的质量随着迭代的进行收敛到某个最优值后，近似解逐渐趋近于一个无意义的被噪声高度污染的向量，图像的质量随之快速下降。如果没有选择一个合适的迭代数则可能求解出无意义的解。MRNSD算法出现半收敛现象的主要原因是边界测量电压中的噪声会随着迭代而扩散，改进算法通过在每步迭代中添加去噪声算法来抑制噪声的扩散。
软闭值方法是一种常用的去噪算法。对任意的向量x∈RN，对于给定μ值，定义软闭值算子Sμ为
 | (14) |
在每步迭代中，将软闭值方法作用于残差向量中，得到以下迭代过程:
 |
虽然rk+1与MRNSD算法中的δUv－δσ不同，但仍可以通过极小化
 | (15) |
另外，根据式(15)中的迭代格式可以得到sk=JTrk和pk=Fksk，则有
 | (16) |
改进算法通过在每步迭代中结合正则化方法与软闭值方法来抑制噪声的扩散，进而克服传统算法的半收敛现象，记为算法3，即快速预处理修正残差范数最陡下降法(Fast Pretreatment Modified Residual Norm Steepest Descent，FPMRNSD)。
算法3??FPMRNSD算法。
输入:Jocobian矩阵J, 迭代初始向量δσ0和迭代终止条件。
输出:重建电导率δσ。
初始化:δσ=δσ0; s=JT(δU－Jδσ); F=diag(δσ); γ=sTFs。
如果k=1, 2, …，则迭代:
p=-Fs
u=C-1Jp
δσold=δσ
δσ=δσold+αp
F=diag(δσ)
r=Sμ(r－αp)
y=JTC-Tu
s=s-αy
γ=sTFs
如
其中：σold为上一次迭代的电导率矩阵；τ为给定参数。
4 实验及对比 为验证算法的有效性，分别引入相关系数和相对误差2个指标对重建图像进行评价。相关系数re定义为
 | (17) |
式中：δσ*为电导率的真实分布向量，
 | (18) |
利用相对误差对不同算法的收敛性能以及重建图像的优劣进行评估。
4.1 仿真实验对比 构建8层各向异性CFRP层压板模型(见图 1)，采用16电极EIT结构。根据CFRP常见损伤类型，构建4种损伤模型, 如图 2所示。在板材表面中心区域设置一个高0.16 cm，半径0.5 cm的圆锥体为冲击损伤模型1；在板材表面一条对角线上设置2个高0.16 cm，半径0.5 cm的圆锥体为冲击损伤模型2；在板材内部中心区域设置一个半径为0.5 cm，高为0.12 cm的圆柱体为分层损伤模型；在板材表面中心区域设置一个4 cm×0.2 cm×0.08 cm的长方体形状为裂纹损伤模型。其中4种模型损伤位置的电导率改变量均为50%。选取被测材料厚度方向中心层XY截面计算灵敏度矩阵J，根据边界电压值的变化，利用CGLS、Tikhonov、MRNSD、PMRNSD和FPMRNSD算法进行图像重建，结果如图 2所示。
 |
| Experimental results of CFRP laminates damage simulation 图 2 CFRP层压板损伤仿真实验结果 Fig. 2 Experimental results of CFRP laminates damage simulation |
| 图选项 |
图 2中的小圆圈表示模型的损伤位置轮廓。从对比图中可以看出，4种算法均可有效地显示出损伤区域，但相比Tikhonov和CGLS算法，MRNSD算法的重建图像中的伪影明显减少，损伤图像轮廓更准确，边缘分辨率更高，有效地提高了重建图像质量。表 1给出了图 2中使用不同算法重建图像的相关系数。可以看到，非负约束算法有效提高了重建图像与真实分布的相关系数，提高了图像质量。
表 1 不同算法的相关系数值 Table 1 Value of correlation coefficients for different algorithms
| 损伤模型 | re | ||||
| CGLS | Tikhonov | MRNSD | PMRNSD | FPMRNSD | |
| 冲击损伤1 | 0.375 2 | 0.517 8 | 0.884 7 | 0.891 6 | 0.882 6 |
| 冲击损伤2 | 0.132 9 | 0.278 6 | 0.275 6 | 0.215 3 | 0.220 3 |
| 分层损伤 | 0.328 3 | 0.436 7 | 0.907 5 | 0.903 1 | 0.910 4 |
| 裂纹损伤 | 0.192 1 | 0.300 8 | 0.503 7 | 0.521 2 | 0.563 7 |
表选项
4.2 抗噪声对比 边界测量电压中的噪声会直接影响成像质量。为验证改进算法抗噪声的有效性，在成像仿真实验中对测量数据施加高斯噪声，即
 | (19) |
式中：L为噪声水平；randn为高斯分布的随机矩阵。L=1%高斯随机噪声是实际测量中的典型噪声级，将其添加到边界测量电压中，通过不同算法得到的重建图像如图 3所示。为了进一步验证5种算法对噪声的鲁棒性，将测量电压中的噪声水平增加到5%，得到重建图像如图 4所示。在1%和5%噪声影响下，不同算法在迭代过程中随迭代次数增加而变化的相对误差和残差范数如图 5和图 6所示。
 |
| 图 3 在1%噪声下CFRP层压板损伤仿真实验结果 Fig. 3 Experimental results of CFRP laminates damage simulation under 1% noise |
| 图选项 |
 |
| 图 4 在5%噪声下CFRP层压板损伤仿真实验结果 Fig. 4 Experimental results of CFRP laminates damage simulation under 5% noise |
| 图选项 |
 |
| 图 5 在1%噪声下CFRP层压板损伤仿真实验数据 Fig. 5 Experimental data of CFRP laminates damage simulation under 1% noise |
| 图选项 |
 |
| 图 6 在5%噪声下CFRP层压板损伤仿真实验数据 Fig. 6 Experimental data of CFRP laminates damage simulation under 5% noise |
| 图选项 |
如图 3、图 4所示，在测量数据中混入了1%噪声后，传统算法CGLS和Tikhonov的重建图像比较算法MRNSD、PMRNSD和FPMRNSD算法更为平滑、模糊，成像目标的轮廓和区域均变大，其抗噪声能力明显低于本文中的改进算法。在5%噪声水平下，传统算法的重建图像伪影明显增多，图像更加模糊。而本文中的改进算法的重建图像更加稳定，目标位置更加准确，在高噪声下具有更好的质量。
重建算法CGLS、MRNSD、PMRNSD和FPMRNSD的相对误差随迭代次数增加的收敛曲线如图 5(a)和图 6(a)所示。CGLS算法的半收敛现象明显，MRNSD算法相对比较稳定，但依然存在半收敛现象，而FPMRNSD算法则克服了CGLS和MRNSD算法的半收敛现象。此外，残差范数是
4.3 CFRP层压板实验 搭建EIT损伤检测实验平台(见图 7)对CFRP层压板进行损伤检测。板材长10 cm，宽10 cm，厚0.3 cm，在板材四周均匀放置16个电极，将电流注入和电压测量的多路复用器与电极相连，实现对电极的循环激励和循环测量。实验中，使用精密电流源(KEITHLEY 6221)在相邻电极施加100 mA电流，使用多功能开关测量单元(KEYSIGHT 34980A)和电枢矩阵开关(KEYSIGHT 34932)进行其他电极电位的采集，基于LabVIEW使用系统软件根据量协议进行开关量的切换。通过接线端子模块(KEYSIGHT 34932T)将电位数据导入成像程序得到最终的成像结果。
 |
| 图 7 EIT碳纤维损伤检测实验平台 Fig. 7 EIT experimental platform for carbon fiber damage detection |
| 图选项 |
损伤的CFRP层压板及不同算法的损伤重建图像如图 8所示。改进算法的重建图像效果更好，损伤轮廓边缘分辨率有明显提高。图 9为CFRP层压板损伤检测实验中不同算法的相关系数对比。可以发现，改进算法明显提高了重建图像的相关系数，提高了重建图像质量。
 |
| 图 8 CFRP层压板损伤检测实验结果 Fig. 8 Experimental results of damage detection for CFRP laminates |
| 图选项 |
 |
| 图 9 CFRP层压板损伤检测实验相关系数对比 Fig. 9 Correlation coefficient comparison of damage detection experiment of CFRP laminates |
| 图选项 |
5 结论 本文通过仿真和实验研究，分析和验证了MRNSD及其改进算法对提高CFRP材料的EIT检测重建图像质量的有效性，得到如下结论：
1) MRNSD算法强制执行非负性约束使逆问题产生更精确的近似解，有效地减少了重建图像的伪影，提高了图像质量。
2) 基于PMRNSD算法，使得灵敏度更加均匀，进一步提高了求解精确度。
3) 基于FPMRNSD算法，在不同噪声强度的影响下均实现了最佳迭代次数的自动更新，有效地克服了迭代中的半收敛现象。改进算法的残差范数值和相对误差值均小于其他算法，表明改进算法有效地抑制了噪声。
参考文献
| [1] | 杨乃宾. 新一代大型客机复合材料结构[J]. 航空学报, 2008, 29(3): 596-603. YANG N B. Composite structures for new generation large commercial jet[J]. Acta Aeronautica et Astronautica Sinica, 2008, 29(3): 596-603. DOI:10.3321/j.issn:1000-6893.2008.03.010 (in Chinese) |
| [2] | 邹达懿, 王鹏飞. 复合材料平尾有限元建模方法研究[J]. 国外电子测量技术, 2012, 31(7): 24-27. ZOU D Y, WANG P F. Investigation of finite element method for composite horizontal tail plane[J]. Foreign Electronic Measurement Technology, 2012, 31(7): 24-27. DOI:10.3969/j.issn.1002-8978.2012.07.007 (in Chinese) |
| [3] | 范文茹, 雷建, 董玉珊, 等. 基于四电极法的CFRP结构损伤检测研究[J]. 仪器仪表学报, 2017, 38(4): 961-968. FAN W R, LEI J, DONG Y S, et al. Damage detection of CFRP laminate structure based on four-probe method[J]. Chinese Journal of Scientific Instrument, 2017, 38(4): 961-968. DOI:10.3969/j.issn.0254-3087.2017.04.022 (in Chinese) |
| [4] | 任吉林, 曾亮, 张丽攀, 等. 碳纤维复合材料涂层厚度涡流法测量的研究[J]. 仪器仪表学报, 2011, 32(12): 2662-2668. REN J L, ZENG L, ZHANG L P, et al. Researh on thickness measureing of carbon fiber composite's coating with eddy current testing[J]. Chinese Journal of Scientific Instrument, 2011, 32(12): 2662-2668. (in Chinese) |
| [5] | 程军.碳纤维复合材料的电磁涡流无损检测技术的研究[D].南京: 南京航空航天大学, 2015: 15-20. CHENG J.Nondestructive testing of carbon fiber reinforced polymer composites using eddy current method[D].Nanjing: Nanjing University of Aeronautics and Astronautics, 2015: 15-20(in Chinese). |
| [6] | 王昵辰, 杨瑞珍, 何赟泽, 等. 多模红外热成像检测碳纤维布加固混凝土粘结缺陷[J]. 仪器仪表学报, 2018, 39(3): 37-44. WANG Z C, YANG R Z, HE Y Z, et al. Detection of interface bonding defects in carbon fiber sheet reinforced concrete using multi-mode infrared thermography[J]. Chinese Journal of Scientific Instrument, 2018, 39(3): 37-44. (in Chinese) |
| [7] | 王珍珍, 任鹏, 程鸿伟, 等. 新型加固用智能碳纤维板及感知性能试验[J]. 中国测试, 2016, 42(3): 113-117. WANG Z Z, REN P, CHENG H W, et al. New smart carbon fiber reinforced polymer plate for strengthening and its sensing performance test[J]. China Measurement & Test, 2016, 42(3): 113-117. (in Chinese) |
| [8] | BALTOPOULOS A, POLYDORIDES N, PAMBAGUIAN L, et al. Damage identification in carbon fiber reinforced polymer plates using electrical resistance tomography mapping[J]. Journal of Composite Materials, 2013, 47(26): 3285-3301. DOI:10.1177/0021998312464079 |
| [9] | FAN W R, WANG H X, XUE Q. Modified sparse regularization for electrical impedance tomography[J]. Review of Scientific Instrument, 2016, 87: 034702. DOI:10.1063/1.4943207 |
| [10] | VAVOULIOTIS A, PAIPETIS A, KOSTOPULOS V. On the fatigue life prediction of CFRP laminates using the electrical resistance change method[J]. Composites Science and Technology, 2011, 71(5): 630-642. DOI:10.1016/j.compscitech.2011.01.003 |
| [11] | CLAY M T, FERREE T C. Weighted regularization in electrical impedance tomography with applications to acute cerebral stroke[J]. IEEE Transactions on Medical Imaging, 2002, 21(6): 629-637. DOI:10.1109/TMI.2002.800572 |
| [12] | GONZALEZ G, KOLEHMAINEN V, SEPPANEN A. Isotropic and anisotropic total variation regularization in electrical impedance tomography[J]. Computers & Mathematics with Applications, 2017, 74(3): 564-576. |
| [13] | 李星, 杨帆, 余晓, 等. 基于自诊断正则化的电阻抗成像逆问题研究[J]. 生物医学工程学杂志, 2018, 35(3): 460-467. LI X, YANG F, YU X, et al. Study on the inverse problem of electrical impedance tomography based on self-diagnosis regularization[J]. Journal of Biomedical Engineering, 2018, 35(3): 460-467. (in Chinese) |
| [14] | 范文茹, 王勃, 李靓瑶, 等. 基于电阻抗层析成像的CFRP结构损伤检测[J]. 北京航空航天大学学报, 2019, 45(11): 2177-2183. FAN W R, WANG B, LI J Y, et al. Damage detection of CFRP structure based on electrical impedance tomography[J]. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(11): 2177-2183. (in Chinese) |
| [15] | 孙磊.碳纤维增强树脂基复合材料涡流无损检测有限元分析[D].厦门: 厦门大学, 2014: 14-21. SUN L.Finite element analysis of eddy current testing of carbon fiber reinforced polymers[D].Xiamen: Xiamen University, 2014: 14-21(in Chinese). |
| [16] | KAUFMAN L. Maximum likelihood, least squares, and penalized least squares for PET[J]. IEEE Transactions on Medical Imaging, 1993, 12(2): 200-214. DOI:10.1109/42.232249 |
| [17] | NAGY J G, STRAKOS Z. Enforcing nonnegativity in image reconstruction algorithms[J]. Mathematical Modeling Estimation and Imaging, 2000, 4121: 182-190. DOI:10.1117/12.402439 |
| [18] | 郭志恒, 律德才, 邵富群. 基于差分灵敏度模型的电容层析成像图像重建方法[J]. 中国电机工程学报, 2012, 32(23): 75-82. GUO Z H, LV D C, SHAO F Q. Image reconstruction method for electrical capacitance tomography based on the difference sensitivity model[J]. Proceedings of the CSEE, 2012, 32(23): 75-82. (in Chinese) |
