CN117147706B - 一种基于准纵波群速度的复合材料各向异性成像方法 - Google Patents

Abstract

本发明公开了一种基于准纵波群速度的复合材料各向异性成像方法。方法包括利用材料的本构模型确定出准纵波相关弹性常数，并根据准纵波相关弹性常数获得三维准纵波群速度离散分布；然后建立三维准纵波群速度方程，并对该方程进行修正，使该方程满足在各向异性条件下使用，且能够匹配任意铺层内群速度分布；利用修正的三维准纵波群速度方程，针对实际检测复合材料样品铺层顺序，针对检测区域进行离散化并据此方程校正延时，完成对复合材料内部缺陷成像。本发明基于各向异性校正的全聚焦算法缺陷成像方法，简化了各向异性校正过程，有助于实现对复合材料的自动化无损检测与评价。

Description

一种基于准纵波群速度的复合材料各向异性成像方法

技术领域

本发明属于高端装备制造领域，尤其是涉及了一种基于准纵波群速度的复合材料各向异性成像方法。

背景技术

纤维增强复合材料（Fiber-Reinforced Polymer，以下称FRP材料）凭借其优异的力学性能，已经被用于航空航天、海洋、新能源等高端装备领域。但纤维增强复合材料的叠层结构使得其在制造或服役过程中容易产生诸如分层、脱粘、气孔、褶皱等缺陷。因此对纤维增强复合材料内部缺陷的表征与评价具有重要的意义。常见的无损检测技术主要有涡流、磁粉、X射线、电子计算机断层扫描(CT)、超声等检测技术。其中，超声无损检测技术以其安全、设备方便携带、能够实现原位检测、检测深度大、检测分辨率高等优点，成为纤维增强复合材料内部缺陷检测的主要技术之一。在超声无损检测技术中，结合相控阵探头（Phasedarray probe）的全矩阵数据采集（Full matrix capture ,FMC）和全聚焦成像（Totalfocus method, TFM）的无损检测技术和评价方法，被称为黄金法则。通过对相控阵探头阵元不同位置收发信号进行采集，对成像区域内每个点依次进行聚焦得到缺陷图像，极大地提高了成像信噪比和横向分辨率，有利于对内部缺陷大小、位置、形状等进行评价。

各向异性材料如碳纤维增强复合材料（Carbon-FRP,以下称CFRP材料）、玻璃纤维增强复合材料（Glass-FRP，以下称GFRP材料）等，根据各向异性材料本构模型计算，声波在内部传播时会发生波形畸变。因此需要对常规全聚焦算法进行延时校正，才能使缺陷位置正确聚焦，从而提高成像质量。

目前延时校正方法主要集中有两种：一为均质化方法，将纤维增强复合材料当成均匀各向异性介质，使用背板反射法（Backwall reflection method，BRM）得到其某一截面处整体群速度分布；二为非均质化方法，将纤维增强复合材料每一铺层单独考虑，借助费马原理或Dijkstra算法计算不同角度下延时。这两种方法都对精确的群速度分布提出了要求。Rokhlin等已经给出了基于弹性常数的纤维增强复合材料准纵波、准横波、纯横波群速度在对称面上计算公式。但由于材料各向异性影响，群速度相对于相速度方向会产生偏斜，即给定相速度方向计算得到群速度以及群速度方向并不等于初始输入的相速度方向。这对于全聚焦算法校正延时需要精确方向对应群速度造成了困难。

发明内容

为了解决背景技术中存在的问题，本发明的目的是提供一种基于准纵波群速度的复合材料各向异性成像方法。

本发明采用的技术方案为：

步骤S1、首先，利用复合材料的本构坐标系确定出四个准纵波相关弹性常数，并根据准纵波相关弹性常数和波矢量获得三维准纵波群速度离散分布；

所述的步骤S1具体为：

步骤S11、对复合材料的单向板进行一次相控阵全矩阵数据采集，得到一组全矩阵数据；

步骤S12、在复合材料本构坐标系中，将平行于纤维排布方向定义为x轴，垂直于纤维排布方向定义为y轴，垂直于xy平面方向定义为z轴，然后根据复合材料本构坐标系建立材料刚度矩阵；

步骤S13、根据步骤S11仅一次采集获得的全矩阵数据，得到材料刚度矩阵中的四个准纵波相关弹性常数C₁₁、C₃₃、C₅₅、C₁₃；

其中，C₁₁、C₃₃、C₅₅、C₁₃分别为材料刚度矩阵中第一排第一列、第三排第三列、第五排第五列、第一排第三列元素的取值；

步骤S14、利用准纵波相关弹性常数获得三维准纵波群速度离散分布；

步骤S2、然后，利用准纵波相关弹性常数和三维准纵波群速度离散分布建立三维准纵波群速度方程；

步骤S3、接着，设置幅值修正因子和相位修正因子，根据幅值修正因子和相位修正因子对三维准纵波群速度方程进行修正，得到修正后的三维准纵波群速度方程；

步骤S4、针对任意铺层顺序复合材料，使用修正的三维准纵波群速度方程对成像区域进行各向异性校正延时，对复合材料的缺陷区域进行全聚焦成像。

所述步骤S14中，三维准纵波群速度离散分布按照以下公式处理得到：

M=arctan((V_gx ²+V_gy ²)^1/2/(V_gz))

V_gx=[n_x/(ρV_p)][-(C₁₁+C₅₅)ρV_p ²+2C₁₁C₅₅n_x ²+(C₁₁C₃₃-C₁₃ ²-2C₁₃C₅₅)n_z ²]/[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²-2ρV_p ²]

V_gy=[n_y/(ρV_p)][-(C₃₃+C₅₅)ρV_p ²+2C₃₃C₅₅(n_y ²+n_z ²)+(C₁₁C₃₃-C₁₃ ²-2C₁₃C₅₅)n_x ²]/[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²-2ρV_p ²]

V_gz=[n_z/(ρV_p)][-(C₃₃+C₅₅)ρV_p ²+2C₃₃C₅₅(n_y ²+n_z ²)+(C₁₁C₃₃-C₁₃ ²-2C₁₃C₅₅)n_x ²]/[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²-2ρV_p ²]

V_p=(2ρ)^-1[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²+(K)^-1/2]^-1/2

n_x=sinθcosφ；n_y=sinθsinφ；n_z=cosθ

其中，M表示面内群速度角度，arctan( )表示求反正切函数，V_gx、V_gy、V_gz分别表示三维准纵波群速度V_g在x、y、z轴上的分量，n_x、n_y、n_z分别表示三维准纵波相速度方向在x、y、z轴上的分量，θ表示入射面内相速度方向，φ表示入射面的极化角，ρ表示复合材料的密度，K表示相速度矢量四次项参数。

所述步骤S2中的三维准纵波群速度方程表示为：

(ρV_g ²)^-1=1/C₁₁N_x ²+1/C₃₃(N_y ²+N_z ²)+D(N_x ²+N_z ²)(N_x ²+N_y ²)

N_x=sinMcosQ；N_y=sinMsinQ；N_z=cosM

其中，V_g表示三维准纵波群速度，N_x、N_y、N_z分别表示三维准纵波群速度方向在x、y、z轴上的分量，D表示非椭圆系数，M表示面内群速度角度，Q表示铺层角度，ρ表示复合材料的密度。

所述的步骤S3具体为：

步骤S31、将三维准纵波群速度方程拆分为二次项和四次项两个部分：

1/(ρV_g ²)=f(M, Q)+ g(M, Q)

其中，f(M, Q)表示三维准纵波群速度方程的二次项部分，g(M, Q)表示三维准纵波群速度方程的四次项部分；

步骤S32、设置幅值修正因子α₁和相位修正因子α₂，根据幅值修正因子α₁和相位修正因子α₂对初始的三维准纵波群速度方程进行修正，修正后的三维准纵波群速度方程表示为：

V_g(M, Q)= {1/ρ [f(M, Q)+ g’(M, Q)] }^1/2

g’(M, Q) =α₁D{[sin(M±α₂)cosQ]²+[cos(M±α₂)]²}{[sin(M±α₂)cosQ]²+[sin(M±α₂)sinQ]²}

其中，V_g(M, Q)表示修正后的三维准纵波群速度；g’(M, Q)表示三维准纵波群速度方程经修正后的四次项部分。

所述步骤S32中，幅值修正因子α₁和相位修正因子α₂的确定方式如下：

首先，给定一组入射面内相速度方向集合[θ]：

[θ]=[θ₁,θ₂, …,θ _l , …,θ_n]

其中，入射面内相速度方向集合[θ]主要由多个入射面内相速度方向θ组成，θ _l 表示第l个入射面内相速度方向，下标n表示入射面内相速度方向的总个数；

然后，按照步骤S1中的公式，根据入射面内相速度方向集合[θ]获得一组三维准纵波群速度集合[V_g]以及一组对应的面内群速度角度集合[M]：

接着，将面内群速度角度集合[M]代入到修正前的三维准纵波群速度方程的二次项部分和四次项部分中，得到一组二次项取值集合[F]和一组四次项取值集合[G]：

最后，按照以下公式处理得到幅值修正因子α₁和相位修正因子α₂：

α₁=max([G])/max(|1/(ρ[V_g]²)-[F]|)

α₂=θ_i-θ_j

其中，max( )表示取最大值函数，||表示取绝对值函数，θ_i表示[G]中最大值对应的入射面内相速度方向，θ_j表示|1/(ρ[V_g]²)- [F]|取最大值时所对应的入射面内相速度方向。

所述的步骤S4具体为：

S41、首先对复合材料的成像区域进行离散以在成像区域中形成若干个像素点，然后利用修正后的三维准纵波群速度方程中的三维准纵波群速度，对成像区域中各个像素点的整条声传播路径进行延时校正，得到相控阵探头中由第i阵元发射、第j阵元接收、成像区域对应像素点(x,y)的整条声传播路径延时t_ij(x,y)：

S42、将全矩阵数据中的所有A扫信号进行叠加，得到校正后的图像：

Im(x,y)=∑_i=1∑_j=1[H(Am(t_ij(x,y)))]

其中，Im(x,y)为成像区域中像素点(x,y)的像素值，H( )表示对数据进行希尔伯特变换，Am( )表示A扫信号的时间与幅值函数。

本发明中，准纵波群速度方程为准纵波群速度（V_g(QP), 以下称V_g）、三维准纵波群速度方向（N_x, N_y, N_z）以及准纵波相关弹性常数(C₁₁, C₃₃, C₅₅, C₁₃)相关方程。在使用本发明方法对复合材料内部缺陷进行成像时，需要得到对应准纵波群速度传播方向与准纵波群速度大小的值。

利用四个准纵波相关弹性常数：C₁₁, C₃₃, C₅₅, C₁₃可得到三维准纵波群速度离散分布。针对某一铺层复合材料，将不加修正的三维准纵波群速度方程与离散群速度分布做差值。不加修正的三维准纵波群速度方程中提取四次项成分，将四次项与差值比较计算得到幅值修正因子α₁与相位修正因子α₂。利用修正因子对三维准纵波群速度方程进行修正，得到针对某一铺层复合材料修正后的三维准纵波群速度方程。重复上述过程，便可以得到任意铺层复合材料内部三维准纵波群速度方程。结合实际检测复合材料层合板铺层顺序，便可以使用该方程对样品进行基于各向异性校正的全聚焦算法缺陷成像。

为了避免利用公式计算只能得到离散群速度分布的影响，本发明在使用基于各向异性校正的全聚焦算法对任意铺层顺序复合材料层合板进行内部缺陷成像时，能够直接得到成像所需要确定群速度角度所对应群速度的值。无需再进行拟合或精细化等操作，便能够得到任意群速度方向上群速度的值，实现全聚焦算法针对任意铺层顺序复合材料层合板的各向异性校正。

本发明能够根据准纵波相关弹性常数得到任意铺层三维准纵波群速度方程，无需求解离散分布并进行拟合或精细化等操作，契合基于三维准纵波群速度的复合材料各向异性成像方法，简化了各向异性校正过程，有助于实现对复合材料的自动化无损检测与评价。

本发明的有益效果为：

1、本发明能够通过四个准纵波相关弹性常数计算三维准纵波群速度方程的修正因子。在基于修正因子对三维准纵波群速度方程四次项进行修正的基础上，得到任意铺层复合材料内部三维准纵波群速度与群速度方向相关方程。实现全聚焦算法针对任意铺层顺序复合材料层合板的各向异性校正。

2、本发明以准纵波相关弹性常数为基础，在修正的基础上完善了三维准纵波群速度方程，实现了任意铺层复合材料内部三维准纵波群速度分布特性的直接计算获取，与非均质各向异性校正全聚焦算法完美契合，大大简化了常规全聚焦算法的各向异性校正过程，有助于实现对复合材料的自动化无损检测与评价。

3、本发明方法计算过程简单，可靠性好，实时性高，能够获得对不同材料任意铺层的纤维增强复合材料三维准纵波群速度分布曲线方程，契合基于各向异性校正的全聚焦算法缺陷成像方法，具有良好的前景。

附图说明

图1是本发明方法实施例的工作流程图；

图2是针对橄榄木材料未修正的三维准纵波群速度方程得到声速分布与计算得到离散分布对比示意图；

图3是针对CFRP材料未修正的三维准纵波群速度方程得到声速分布与计算得到离散分布对比示意图；

图4是本发明方法实施例的理论离散值与公式计算值、四次项值、和二次项值结果对比示意图；

图5是对复合材料层合板实施全聚焦算法校正时对成像区域离散示意图；

图6是本发明方法实施例中对样品一侧钻孔缺陷使用修正与未修正的三维准纵波群速度方程成像FWHM结果对比示意图；

图7是本发明方法实施例中对样品二侧钻孔缺陷使用修正与未修正的三维准纵波群速度方程成像FWHM结果对比示意图。

具体实施方式

下面结合具体实施案例对本发明进行详细说明，以下实施案例将有助于本领域的技术人员进一步理解本发明，但不以任何形式限制本发明。

本发明的实施例及其实施过程、情况如下：

本实施例所使用的复合材料样品由碳纤维增强聚合物基预浸料制备而成，单向板由全0度层预浸料制成。复合材料单向板本构模型为横观各向同性，具有三个弹性对称面，其中纤维垂直面为各向同性面。

如图1所示，基于准纵波群速度的复合材料各向异性成像方法具体步骤如下：

步骤S1、首先，利用复合材料的本构坐标系确定出四个准纵波相关弹性常数，并根据准纵波相关弹性常数和波矢量获得三维准纵波群速度离散分布；

步骤S1具体为：

步骤S11、使用超声相控阵采集设备对复合材料的单向板进行一次相控阵全矩阵数据采集，得到一组全矩阵数据；

具体实施中，使用超声相控阵采集设备的相控阵探头对复合材料的内部缺陷进行超声全矩阵数据采集和全聚焦成像，超声全矩阵数据采集指的是使用相控阵中某一阵元发射，所有阵元接收，依次循环发射所有阵元的过程，将超声全矩阵数据采集所得到的数据矩阵称为全矩阵数据。

步骤S12、在复合材料本构坐标系中，将平行于纤维排布方向定义为x轴，垂直于纤维排布方向定义为y轴，垂直于xy平面方向定义为z轴，声波在xz对称面入射时，因为各向异性影响，存在准纵波、准横波和纯横波三个波阵面，接着根据复合材料本构坐标系建立材料刚度矩阵；

步骤S13、根据步骤S11仅一次采集获得的全矩阵数据，提取不同传播角度下A扫信号的准纵波与准横波的回波渡越时间，然后根据回波渡越时间得到材料刚度矩阵中的四个准纵波相关弹性常数C₁₁、C₃₃、C₅₅、C₁₃；

其中，C₁₁、C₃₃、C₅₅、C₁₃分别为材料刚度矩阵中第一排第一列、第三排第三列、第五排第五列、第一排第三列元素的取值；

步骤S14、利用准纵波相关弹性常数得到三维准纵波相速度V_p：

V_p=(2ρ)^-1[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²+(K)^-1/2]^-1/2

K=[(C₃₃-C₅₅)(n_y ²+n_z ²)-(C₁₁-C₅₅)n_x ²]²+4(C₁₃+C₅₅)²(n_y ²+n_z ²)n_x ²

n_x=sinθcosφ

n_y=sinθsinφ

n_z=cosθ

其中，n_x、n_y、n_z分别表示三维准纵波相速度方向在x、y、z轴上的分量，θ表示入射面内相速度方向，φ表示入射面的极化角，ρ表示复合材料的密度，K表示相速度矢量四次项参数。

步骤S15、根据三维准纵波相速度和波矢量获得三维准纵波群速度离散分布；

步骤S15中，三维准纵波群速度离散分布按照以下公式处理得到：

M=arctan((V_gx ²+V_gy ²)^1/2/(V_gz))

V_gx=[n_x/(ρV_p)][-(C₁₁+C₅₅)ρV_p ²+2C₁₁C₅₅n_x ²+(C₁₁C₃₃-C₁₃ ²-2C₁₃C₅₅)n_z ²]/[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²-2ρV_p ²]

V_gy=[n_y/(ρV_p)][-(C₃₃+C₅₅)ρV_p ²+2C₃₃C₅₅(n_y ²+n_z ²)+(C₁₁C₃₃-C₁₃ ²-2C₁₃C₅₅)n_x ²]/[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²-2ρV_p ²]

V_gz=[n_z/(ρV_p)][-(C₃₃+C₅₅)ρV_p ²+2C₃₃C₅₅(n_y ²+n_z ²)+(C₁₁C₃₃-C₁₃ ²-2C₁₃C₅₅)n_x ²]/[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²-2ρV_p ²]

其中，M表示面内群速度角度，arctan( )表示求反正切函数，V_gx、V_gy、V_gz分别表示三维准纵波群速度V_g在x、y、z轴上的分量；

在一个确定平面内，给定一个入射面内相速度方向θ，便可以求解得到对应群速度大小与群速度方向。但由于各向异性影响，群速度与相速度方向并不相等，因此并不能直接得到想要的对应群速度方向的群速度值。

步骤S2、然后，利用准纵波相关弹性常数和三维准纵波群速度离散分布建立三维准纵波群速度方程；

步骤S2中的三维准纵波群速度方程为：

(ρV_g ²)^-1=1/C₁₁N_x ²+1/C₃₃(N_y ²+N_z ²)+D(N_x ²+N_z ²)(N_x ²+N_y ²)

D=[(C₁₁+C₃₃)-2(C₁₃+2C₅₅)]/C₁₁C₃₃

N_x=sinMcosQ

N_y=sinMsinQ

N_z=cosM

其中，V_g表示三维准纵波群速度，N_x、N_y、N_z分别表示三维准纵波群速度方向在x、y、z轴上的分量，D表示非椭圆系数，M表示面内群速度角度，Q表示铺层角度，ρ表示复合材料的密度。

如图2所示，其中横坐标V_g[100]表示三维准纵波群速度V_g在x轴方向分量大小，纵坐标V_g[001]表示V_g在z轴方向分量大小，所提出的三维准纵波群速度方程在对弱各向异性材料如橄榄木等进行准纵波群速度表征时，具有较好的贴合效果。但在图3中，其中横坐标V_g[100]表示V_g在x轴方向分量大小，纵坐标V_g[001]表示V_g在z轴方向分量大小，针对CFRP等强各向异性材料表征时，在某些角度范围存在一定误差，因此该方程需要进行修正。

步骤S3、接着，设置幅值修正因子和相位修正因子，根据幅值修正因子和相位修正因子对三维准纵波群速度方程进行修正，得到修正后的三维准纵波群速度方程；

具体地，先根据步骤S14中的公式确定出一组三维准纵波相速度，然后将该组三维准纵波相速度代入步骤S15的公式中，可得到一组三维准纵波群速度和对应的一组面内群速度角度，将面内群速度角度带入三维准纵波群速度方程中可得到一组方程值，将三维准纵波群速度方程中二次项的值与离散分布计算值做差与四次项进行比较，计算得到幅值修正因子和相位修正因子，最终完成对三维准纵波群速度方程的修正；

具体实施中，步骤S3具体包括以下步骤：

步骤S31、将三维准纵波群速度方程拆分为二次项和四次项两个部分：

1/(ρV_g ²)=f(M, Q)+ g(M, Q)

f(M, Q)=1/C₁₁N_x ²+1/C₃₃(N_y ²+N_z ²)

g(M, Q)=D(N_z ²N_x ²+N_y ²N_x ²)

其中，f(M, Q)表示三维准纵波群速度方程的二次项部分，g(M, Q)表示三维准纵波群速度方程的四次项部分，V_g表示三维准纵波群速度，N_x、N_y、N_z分别表示三维准纵波群速度方向在x、y、z轴上的分量，D表示非椭圆系数，ρ表示复合材料的密度；

如图4所示，f(M, Q)提供方程主体部分，g(M, Q)提供方程在小范围内波动部分。因此需要对g(M, Q)进行修正，使g(M, Q)满足强各向异性介质。

步骤S32、设置幅值修正因子α₁和相位修正因子α₂，根据幅值修正因子α₁和相位修正因子α₂对初始的三维准纵波群速度方程进行修正，修正后的三维准纵波群速度方程表示为：

V_g(M, Q)= {1/ρ [f(M, Q)+ g’(M, Q)] }^1/2

f(M, Q)=1/C₁₁N_x ²+1/C₃₃(N_y ²+ N_z ²)

g’(M, Q) =α₁D{[sin(M±α₂)cosQ]²+[cos(M±α₂)]²}{[sin(M±α₂)cosQ]²+[sin(M±α₂)sinQ]²}

其中，V_g(M, Q)表示修正后的三维准纵波群速度；g’(M, Q)表示三维准纵波群速度方程经修正后的四次项部分，M表示面内群速度角度，Q表示铺层角度；

具体实施中，根据幅值修正因子α₁对D进行幅值修正，根据相位修正因子α₂对M进行相位修正。

步骤S32中，幅值修正因子α₁和相位修正因子α₂的确定方式如下：

首先，给定一组入射面内相速度方向集合[θ]：

[θ]=[θ₁,θ₂, …,θ _l , …,θ_n]

其中，θ _l 表示第l个入射面内相速度方向，下标n表示入射面内相速度方向的总个数；

然后，按照步骤S1中的公式根据入射面内相速度方向集合[θ]获得一组三维准纵波群速度集合[V_g]以及一组对应的面内群速度角度集合[M]：

[V_g]=[V_g1, V_g2, …, V_gl , …,V_gn]

[M]=[M₁, M₂, …, M _l , …, M_n]

其中，V_gl 表示第l个三维准纵波群速度，M _l 表示第l个面内群速度角度；

接着，将面内群速度角度集合[M]代入到三维准纵波群速度方程的二次项和四次项中，得到一组二次项取值集合[F]和一组四次项取值集合[G]：

[F]=[ f₁, f₂, …, f _l , …,f_n]

[G]=[ g₁, g₂, …, g _l , …,g_n]

其中，f _l 表示第l个二次项部分的取值，g _l 表示第l个四次项部分的取值；

最后，按照以下公式处理得到幅值修正因子α₁和相位修正因子α₂：

α₁=max([G])/max(|1/(ρ[V_g]²)-[F]|)

α₂=θ_i-θ_j

其中，max( )表示取一组数据中的最大值函数，||表示取绝对值函数，θ_i表示[G]中最大值对应的入射面内相速度方向，θ_j表示|1/(ρ[V_g]²)- [F]|取最大值时所对应的入射面内相速度方向。

步骤S4、针对任意铺层顺序复合材料层合板，使用修正的三维准纵波群速度方程对成像区域进行各向异性校正延时，对复合材料的缺陷区域进行全聚焦成像。

步骤S4具体为：

S41、首先对复合材料的成像区域进行离散以在成像区域中形成若干个像素点，然后利用修正后的三维准纵波群速度方程，对成像区域中各个像素点的整条声传播路径进行延时校正：

t_ij(x,y)= ∑_k=1[w_k/ V_g(M, Q)]

w_k=(d²+l_k ²)^1/2

M=arctan(d/ l_k)

其中，t_ij(x,y)表示相控阵探头中由第i阵元发射、第j阵元接收、成像区域对应像素点(x,y)的整条声传播路径延时，∑表示累加求和函数，w_k表示在复合材料的第k层内声传播路径长度，k=1, 2, …, N，N表示复合材料实际层的总层数，V_g(M, Q)表示修正后的三维准纵波群速度，d表示纵向离散层间距，l_k表示复合材料的第k层内声传播路径长度的横向分量，具体实施中，离散层与材料实际层重合；

像素点的整条声传播路径是指发射阵元到像素点以及像素点到接收阵元之间的声波传播路径。

S42、将全矩阵数据中的所有A扫信号进行叠加，得到校正后的图像：

Im(x,y)=∑_i=1∑_j=1[H(Am(t_ij(x,y)))]

其中，Im(x,y)为感兴趣成像区域中像素点(x,y)的像素值，H( )表示对数据进行希尔伯特变换，Am( )表示A扫信号的时间与幅值函数。

针对全矩阵数据中任意发射接收阵元组合的A扫信号，针对区域中某一点聚焦，考虑离散区域每一层间的折射，利用步骤S3得到修正后的三维准纵波群速度方程校正每一段延时，并对整条声传播路径上延时进行叠加。

根据待检测复合材料样品的不同铺层（例如Q=30°, 45°, 60°, 90°等），均可利用步骤S3得到的修正后三维准纵波群速度方程。利用该方程便可以得到对应铺层的完整准纵波群速度分布，据此完成全聚焦算法的校正成像。在全聚焦成像中，需要对复合材料的成像区域进行离散，如图5所示，d表示纵向离散层间距，M表示面内群速度角度，l₁表示复合材料的第1层内声传播路径长度的横向分量，i表示发射阵元的序数，j均表示接收阵元的序数，N表示复合材料实际层的总层数，其中在纵向上离散间隔依据厚度进行设置，为了保证考虑层间折射影响，在层间必须分布离散点；横向离散间隔则根据精度进行设置。

本发明计算过程简单可靠，效率高，仅利用四个准纵波相关弹性常数，便可以得到三维准纵波群速度方程。避免了直接利用公式计算只能得到离散群速度分布造成影响，无需再进行拟合或精细化等操作，便能够得到任意群速度方向上群速度的值。

且该算法实现了任意铺层复合材料内部三维准纵波群速度分布特性的直接计算获取，与非均质各向异性校正全聚焦算法完美契合，大大简化了常规全聚焦算法的各向异性校正过程，有助于实现对复合材料的自动化无损检测与评价。

下面通过本发明方法对具有不同铺层顺序的CFRP复合材料层合板进行内部缺陷检测，以验证本申请方法准确性与可靠性：

实验使用CFRP层合板，由碳纤维/环氧树脂预浸料（USN15000/6511）通过热压工艺制成。样品一总厚度为20mm，铺层顺序为全单向板，在材料中间距离底部10mm处有一直径为3mm大小侧钻孔，用以验证修正的三维准纵波群速度方程准确性。样品二总厚度为12mm，铺层顺序为[-45/0/45/90]₂₀，在距离材料底部分别为3mm的位置各有一个直径为2mm的侧钻孔，用以验证使用修正的三维准纵波群速度方程对全聚焦算法进行各向异性校正的准确性。样品一与样品二的密度均为1590kg/m3.

实验使用型号为5L64-38.4相控阵探头，中心频率为5MHz，共包含64阵元，阵元间距为0.6mm。信号采集设备为64/64 OEM-PA便携式相控阵采集仪。

首先利用相控阵探头在复合材料单向板平行于纤维排布方向采集全矩阵数据，根据不同角度A扫信号回波渡越时间，得到所需的四个准纵波相关弹性常数C₁₁、C₃₃、C₅₅、C₁₃。

然后，利用四个准纵波相关弹性常数C₁₁、C₃₃、C₅₅、C₁₃，给定一组入射面内相速度方向集合[θ]值计算得到某一面内三维准纵波群速度离散分布，即一组三维准纵波群速度与群速度方向对应角度值。

接着，得到对应面内的三维准纵波群速度方程；并对三维准纵波群速度方程进行修正，使该方程能够满足在强各向异性条件下使用。根据步骤S1的公式，将所得到的一组三维准纵波群速度代入三维准纵波群速度方程中，得到一组由三维准纵波群速度方程计算得到的值。将三维准纵波群速度方程二次项与三维准纵波群速度离散分布的差值与三维准纵波群速度方程四次项进行比较，计算得到幅值修正因子与相位修正因子。

将修正的三维准纵波群速度方程代入全聚焦算法，通过将成像区域离散，完成各向异性校正，最终实现对复合材料缺陷的成像结果。

为了进一步证明本发明方法修正的三维准纵波群速度方程以及对缺陷成像的准确性可靠性，首先针对样品一将本发明方法与未修正的三维准纵波群速度方程进行成像的结果对比，其次针对样品二将本发明方向基于各向异性校正的全聚焦算法与不加校正的全聚焦算法（各向同性全聚焦算法）成像结果进行对比。使用半高全宽（full width at halfmaxima，FWHM）进行成像质量比较。

本发明方法对于CFRP样品一的侧钻孔缺陷使用修正与未修正的三维准纵波群速度方程进行成像对比，FWHM结果如图6所示。本发明方法基于修正的三维准纵波群速度方程对全聚焦算法进行校正成像，FWHM结果明显优于使用未修正的三维准纵波群速度方程，同时因为准确的校正，与各向同性TFM相比分辨率也得到了提升。

本发明方法对于具有不同铺层顺序的CFRP样品二的侧钻孔缺陷同样使用修正与未修正的三维准纵波群速度方程进行成像对比，FWHM结果如图7所示。本发明方法基于修正的三维准纵波群速度方程对全聚焦算法进行校正成像，FWHM结果表明首先本发明方法在缺陷成像的强度上高于使用各向同性TFM和使用未修正的方程，其次在半高全宽上本申请方法同样优于其余两种方法，证明本发明方法对复合材料缺陷进行全聚焦校正成像具有良好的图像信噪比和分辨率。

Claims (2)

1.一种基于准纵波群速度的复合材料各向异性成像方法，其特征在于，包括以下步骤：

步骤S1、首先，利用复合材料的本构坐标系确定出四个准纵波相关弹性常数，并根据准纵波相关弹性常数获得三维准纵波群速度离散分布；

所述的步骤S1具体为：

步骤S11、对复合材料进行一次数据采集，得到全矩阵数据，然后根据复合材料本构坐标系建立材料刚度矩阵；

步骤S12、根据步骤S11仅一次采集获得的全矩阵数据，得到材料刚度矩阵中的四个准纵波相关弹性常数C₁₁、C₃₃、C₅₅、C₁₃；

其中，C₁₁、C₃₃、C₅₅、C₁₃分别为材料刚度矩阵中第一排第一列、第三排第三列、第五排第五列、第一排第三列元素的取值；

步骤S13、利用准纵波相关弹性常数和入射面内相速度方向θ得到三维准纵波相速度V_p：

V_p=(2ρ)^-1[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²+(K)^-1/2]^-1/2

K=[(C₃₃-C₅₅)(n_y ²+n_z ²)-(C₁₁-C₅₅)n_x ²]²+4(C₁₃+C₅₅)²(n_y ²+n_z ²)n_x ²

n_x=sinθcosφ

n_y=sinθsinφ

n_z=cosθ

其中，n_x、n_y、n_z分别表示三维准纵波相速度方向在x、y、z轴上的分量，θ表示入射面内相速度方向，φ表示入射面的极化角，ρ表示复合材料的密度，K表示相速度矢量四次项参数；

步骤S14、根据三维准纵波相速度V_p获得三维准纵波群速度离散分布；

所述步骤S14中，三维准纵波群速度离散分布按照以下公式处理得到：

M=arctan((V_gx ²+V_gy ²)^1/2/(V_gz))

V_gx=[n_x/(ρV_p)][-(C₁₁+C₅₅)ρV_p ²+2C₁₁C₅₅n_x ²+(C₁₁C₃₃-C₁₃ ²-2C₁₃C₅₅)n_z ²]/[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²-2ρV_p ²]

V_gy=[n_y/(ρV_p)][-(C₃₃+C₅₅)ρV_p ²+2C₃₃C₅₅(n_y ²+n_z ²)+(C₁₁C₃₃-C₁₃ ²-2C₁₃C₅₅)n_x ²]/[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²-2ρV_p ²]

V_gz=[n_z/(ρV_p)][-(C₃₃+C₅₅)ρV_p ²+2C₃₃C₅₅(n_y ²+n_z ²)+(C₁₁C₃₃-C₁₃ ²-2C₁₃C₅₅)n_x ²]/[(C₃₃+C₅₅)(n_y ²+n_z ²)+(C₁₁+C₅₅)n_x ²-2ρV_p ²]

其中，M表示面内群速度角度，arctan( )表示求反正切函数，V_gx、V_gy、V_gz分别表示三维准纵波群速度V_g在x、y、z轴上的分量；

步骤S2、然后，利用三维准纵波群速度离散分布建立三维准纵波群速度方程；

所述步骤S2中的三维准纵波群速度方程为：

(ρV_g ²)^-1=1/C₁₁N_x ²+1/C₃₃(N_y ²+N_z ²)+D(N_x ²+N_z ²)(N_x ²+N_y ²)

D=[(C₁₁+C₃₃)-2(C₁₃+2C₅₅)]/C₁₁C₃₃

N_x=sinMcosQ

N_y=sinMsinQ

N_z=cosM

其中，V_g表示三维准纵波群速度，N_x、N_y、N_z分别表示三维准纵波群速度方向在x、y、z轴上的分量，D表示非椭圆系数，M表示面内群速度角度，Q表示铺层角度，ρ表示复合材料的密度；

步骤S3、接着，设置幅值修正因子和相位修正因子，根据幅值修正因子和相位修正因子对三维准纵波群速度方程进行修正，得到修正后的三维准纵波群速度方程；

所述的步骤S3具体为：

步骤S31、将三维准纵波群速度方程拆分为二次项和四次项两个部分：

1/(ρV_g ²)=f(M, Q)+ g(M, Q)

f(M, Q)=1/C₁₁N_x ²+1/C₃₃(N_y ²+N_z ²)

g(M, Q)=D(N_z ²N_x ²+N_y ²N_x ²)

其中，f(M, Q)表示三维准纵波群速度方程的二次项部分，g(M, Q)表示三维准纵波群速度方程的四次项部分；

步骤S32、设置幅值修正因子α₁和相位修正因子α₂，根据幅值修正因子α₁和相位修正因子α₂对初始的三维准纵波群速度方程进行修正，修正后的三维准纵波群速度方程表示为：

V_g(M, Q)= {1/ρ [f(M, Q)+ g’(M, Q)] }^1/2

f(M, Q)=1/C₁₁N_x ²+1/C₃₃(N_y ²+N_z ²)

g’(M,Q)=α₁D{[sin(M±α₂)cosQ]²+[cos(M±α₂)]²}{[sin(M±α₂)cosQ]²+[sin(M±α₂)sinQ]²}

其中，V_g(M, Q)表示修正后的三维准纵波群速度；g’(M, Q)表示三维准纵波群速度方程经修正后的四次项部分，M表示面内群速度角度，Q表示铺层角度 ；

所述步骤S32中，幅值修正因子α₁和相位修正因子α₂的确定方式如下：

首先，给定一组入射面内相速度方向集合[θ]；

所述的入射面内相速度方向集合[θ]主要由多个入射面内相速度方向θ组成；

然后，根据入射面内相速度方向集合[θ]获得一组三维准纵波群速度集合[V_g]以及一组对应的面内群速度角度集合[M]；

接着，将面内群速度角度集合[M]代入到三维准纵波群速度方程的二次项部分和四次项部分中，得到一组二次项取值集合[F]和一组四次项取值集合[G]；

最后，按照以下公式处理得到幅值修正因子α₁和相位修正因子α₂：

α₁=max([G])/max(|1/(ρ[V_g]²)-[F]|)

α₂=θ_i-θ_j

其中，max( )表示取最大值函数，||表示取绝对值函数，θ_i表示[G]中最大值对应的入射面内相速度方向，θ_j表示|1/(ρ[V_g]²)- [F]|取最大值时所对应的入射面内相速度方向；

步骤S4、针对任意铺层顺序复合材料，使用修正的三维准纵波群速度方程对成像区域进行各向异性校正延时，对复合材料的缺陷区域进行全聚焦成像。

2.根据权利要求1所述的一种基于准纵波群速度的复合材料各向异性成像方法，其特征在于：

所述的步骤S4具体为：

S41、首先利用修正后的三维准纵波群速度方程中的三维准纵波群速度，对复合材料的成像区域中各个像素点的整条声传播路径进行延时校正，得到相控阵探头中由第i阵元发射、第j阵元接收、成像区域对应像素点(x,y)的整条声传播路径延时t_ij (x,y)；

S42、将全矩阵数据中的所有A扫信号进行叠加，得到校正后的图像。