CN113376256A - 基于变步长网格模型小波包分量波形的厚度遍历反演法 - Google Patents

基于变步长网格模型小波包分量波形的厚度遍历反演法 Download PDF

Info

Publication number
CN113376256A
CN113376256A CN202110629436.7A CN202110629436A CN113376256A CN 113376256 A CN113376256 A CN 113376256A CN 202110629436 A CN202110629436 A CN 202110629436A CN 113376256 A CN113376256 A CN 113376256A
Authority
CN
China
Prior art keywords
grid
thickness
wavelet packet
layer structure
coordinates
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110629436.7A
Other languages
English (en)
Other versions
CN113376256B (zh
Inventor
周世圆
胡怡
孙晓莹
刘逯航
郑翀
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beijing Institute of Technology BIT
Original Assignee
Beijing Institute of Technology BIT
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beijing Institute of Technology BIT filed Critical Beijing Institute of Technology BIT
Priority to CN202110629436.7A priority Critical patent/CN113376256B/zh
Publication of CN113376256A publication Critical patent/CN113376256A/zh
Application granted granted Critical
Publication of CN113376256B publication Critical patent/CN113376256B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/04—Analysing solids
    • G01N29/07—Analysing solids by measuring propagation velocity or propagation time of acoustic waves
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
    • G01N29/4409—Processing the detected response signal, e.g. electronic circuits specially adapted therefor by comparison
    • G01N29/4418—Processing the detected response signal, e.g. electronic circuits specially adapted therefor by comparison with a model, e.g. best-fit, regression analysis
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
    • G01N29/4409—Processing the detected response signal, e.g. electronic circuits specially adapted therefor by comparison
    • G01N29/4427—Processing the detected response signal, e.g. electronic circuits specially adapted therefor by comparison with stored values, e.g. threshold values
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N29/00—Investigating or analysing materials by the use of ultrasonic, sonic or infrasonic waves; Visualisation of the interior of objects by transmitting ultrasonic or sonic waves through the object
    • G01N29/44—Processing the detected response signal, e.g. electronic circuits specially adapted therefor
    • G01N29/4472—Mathematical theories or simulation
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00—Indexing codes associated with group G01N29/00
    • G01N2291/01—Indexing codes associated with the measuring variable
    • G01N2291/011—Velocity or travel time
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00—Indexing codes associated with group G01N29/00
    • G01N2291/02—Indexing codes associated with the analysed material
    • G01N2291/023—Solids
    • G01N2291/0231—Composite or layered materials
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01N—INVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N2291/00—Indexing codes associated with group G01N29/00
    • G01N2291/02—Indexing codes associated with the analysed material
    • G01N2291/028—Material parameters
    • G01N2291/02854—Length, thickness

Landscapes

  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Pathology (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Analytical Chemistry (AREA)
  • Biochemistry (AREA)
  • Health & Medical Sciences (AREA)
  • Immunology (AREA)
  • Signal Processing (AREA)
  • Engineering & Computer Science (AREA)
  • Algebra (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Mathematical Physics (AREA)
  • Pure & Applied Mathematics (AREA)
  • Acoustics & Sound (AREA)
  • Investigating Or Analyzing Materials By The Use Of Ultrasonic Waves (AREA)

Abstract

本发明提供一种基于变步长网格模型小波包分量波形的厚度遍历反演法，针对多层粘接结构不同层大差别区域物理参数，采用变步长网格模型，在保证反演精度的前提下大幅降低了模型规模；通过小波包分解方法对实测全波形和正演观测空间信号时域全波形进行处理，采用M层分解小波包分量波形即可将时域全波形数据量缩减到原始信号的1/2M，在指数级减少波形数据量的同时滤除了高频噪声信号；利用中间层厚度区间极小的特点，提出厚度遍历反演方法；通过上述三项措施，实现了中间薄层的高效高精度反演测量。

Description

基于变步长网格模型小波包分量波形的厚度遍历反演法
技术领域
本发明属于测控技术领域，尤其涉及一种基于变步长网格模型小波包分量波形的厚度遍历反演法。
背景技术
多层结构材料因为其具有良好的机械性能和热力学性能，被广泛应用于航空航天，特种设备，工业通用技术设备，军用武器设备等众多领域。多层结构粘接质量好坏直接影响到结构整体性能，对整体结构的使用性、稳定性以及使用寿命都有重要影响。其中，中间粘接层的厚度异常超差(过厚或者过薄)会导致使用过程中整个结构的失效，造成重大事故和不可估量的经济损失。
由于多层结构粘接层厚度较薄，且多层结构各层之间声阻抗差异较大，在粘接层-下层结构交界面产生的二界面反射回波较弱，无法得到明显的界面回波信息，因此通过常规使用的脉冲回波法无法测量粘接薄层的厚度。
综上所述，亟待研究新型厚度测量方法，解决上述多层结构中间粘接薄层厚度测量难题。地球物理勘探领域复杂地质结构的全波形反演方法能够有效的重构地下物理特性参数，可通过采用全波形反演方法测量多层结构粘接层的厚度，解决其厚度异常超差检测难题。
由于多层结构粘接层厚度较薄，对正演模型的网格剖分不宜采用稀疏网格，而对多层结构各层模型均采用细密网格划分将导致整个网格系统网格量和反演计算量大幅增加，反演效率较低。此外，全波形反演方法使用的实测全波形和正演观测空间信号时域全波形的数据量较大，对实测全波形和正演观测空间信号时域全波形的大量计算同样降低了反演效率。因此在保证反演精度的同时降低全波形反演涉及的计算量对于提高反演效率具有重大意义。
发明内容
为解决上述问题，本发明提供一种基于变步长网格模型小波包分量波形的厚度遍历反演法，在已知三层结构各层纵波声速、横波声速和密度，及上、下层的厚度的情况下，实现了中间薄层的厚度测量；通过建立变步长网格模型与利用小波包分量波形，指数级的减少了反演计算量，有效提高了反演计算效能。
一种基于变步长网格模型小波包分量波形的厚度遍历反演法，待测粘接层位于上层结构和下层结构之间，形成三层结构，待测粘接层的厚度遍历反演法包括以下步骤：
S1：将激励信号从三层结构的上表面垂直入射，得到三层结构的实测回波信号dobs；
S2：构建三层结构的三维网格模型，设定三层结构各层的纵波声速、横波声速、密度以及厚度，同时设定三维网格模型各网格节点处的速度场和应力场的初始值，其中，待测粘接层的网格步长小于上层结构和下层结构的网格步长，且待测粘接层的厚度为可能达到的厚度极大值，网格节点包括网格的顶点以及网格边界的中点；
S3：将激励模拟信号从三维网格模型的上表面垂直入射，得到三层结构的正演回波信号dcal，并从正演回波信号dcal提取上层结构与待测粘接层之间的界面的观测回波信号S1m；
S4：将观测回波信号S1m与实测回波信号dobs分别进行三次以上的小波包分解，对应得到小波包分量观测波形与小波包分量实测波形；
S5：将小波包分量观测波形与小波包分量实测波形之间残差的L1范数作为目标函数值，判断所述目标函数值是否小于设定阈值或迭代次数是否达到设定上限，若满足其中一个条件，则待测粘接层当前设定的厚度为所求的厚度值，若均不满足，则进入步骤S6；
S6：按照设定步长减小三维网格模型中待测粘接层的厚度，并重新划分网格和设定三维网格模型各网格节点处的速度场和应力场，实现三维网格模型的更新，再采用更新后的三维网格模型重新执行步骤S3～S5，直到得到待测粘接层的厚度值。
进一步地，以三层结构厚度方向为笛卡尔坐标系x方向，以三层结构长度方向为笛卡尔坐标系z方向，所述速度场包括x方向的速度分量与z方向的速度分量，所述应力场包括切应力、x方向的正应力以及z方向的正应力，步骤S6中所述三维网格模型各网格节点处的速度场和应力场的设定方法为：
Figure BDA0003102082360000031
Figure BDA0003102082360000032
Figure BDA0003102082360000033
Figure BDA0003102082360000034
Figure BDA0003102082360000035
其中，i为网格顶点在x方向上的编号，j为网格顶点在z方向的编号，
Figure BDA0003102082360000036
为坐标为(i,j)的网格节点处x方向的空间差分算子，
Figure BDA0003102082360000037
为坐标为(i,j)的网格节点处z方向的空间差分算子，
Figure BDA0003102082360000038
为坐标为(i,j)的网格节点处的时间差分算子，ρi,j+1/2为坐标为(i,j+1/2)的网格节点处的密度，
Figure BDA0003102082360000041
为坐标为(i,j+1/2)的网格节点处x方向的速度分量，
Figure BDA0003102082360000042
为坐标为(i,j+1/2)的网格节点处x方向的正应力，
Figure BDA0003102082360000043
为坐标为(i,j+1/2)的网格节点处的切应力，ρi+1/2,j为坐标为(i+1/2,j)的网格节点处的密度，
Figure BDA0003102082360000044
为坐标为(i+1/2,j)的网格节点处z方向的速度分量，
Figure BDA0003102082360000045
为坐标为(i+1/2,j)的网格节点处z方向的正应力，
Figure BDA0003102082360000046
为坐标为(i+1/2,j)的网格节点处的切应力，
Figure BDA0003102082360000047
为坐标为(i+1/2,j+1/2)的网格节点处x方向的正应力，
Figure BDA0003102082360000048
为坐标为(i+1/2,j+1/2)的网格节点处x方向的速度分量，λi+1/2,j+1/2为坐标为(i+1/2,j+1/2)的网格节点处的一阶拉梅常数，μi+1/2,j+1/2为坐标为(i+1/2,j+1/2)的网格节点处的二阶拉梅常数，
Figure BDA0003102082360000049
为坐标为(i+1/2,j+1/2)的网格节点处z方向的速度分量，
Figure BDA00031020823600000410
为坐标为(i+1/2,j+1/2)的网格节点处z方向的正应力，
Figure BDA00031020823600000411
为坐标为(i,j)的网格节点处的切应力，μi,j为坐标为(i,j)的网格节点处的二阶拉梅常数，
Figure BDA00031020823600000412
为坐标为(i,j)的网格节点处x方向的速度分量，
Figure BDA00031020823600000413
为坐标为(i,j)的网格节点处z方向的速度分量。
进一步地，空间差分算子
Figure BDA00031020823600000414
空间差分算子
Figure BDA00031020823600000415
以及时间差分算子
Figure BDA00031020823600000416
的获取方法为：
将激励模拟信号在各网格顶点处的波动进行泰勒展开，得到各网格顶点处的波场值，再根据波场值获取各网格顶点处的空间差分算子
Figure BDA00031020823600000417
空间差分算子
Figure BDA00031020823600000418
以及时间差分算子
Figure BDA00031020823600000419
进一步地，上层结构与待测粘接层之间的界面为上交界面，则上交界面处的网格顶点的空间差分算子
Figure BDA00031020823600000420
和
Figure BDA00031020823600000421
的计算方法为：
依次在x方向上取四个网格顶点执行空间差分算子获取操作，且每次选出的四个网格顶点所在的三个网格单元中，至少有一个、至多有两个属于待测粘接层；其中，所述空间差分算子获取操作为：
将四个网格顶点沿上层结构至待测粘接层的方向分别记为(m-1,j)、(m,j)、(m+1,j)以及(m+2,j)，顶点(m-1,j)和顶点(m,j)所在的网格记为第一网格，顶点(m,j)和顶点(m+1,j)所在的网格记为第二网格，顶点(m+1,j)和顶点(m+2,j)所在的网格记为第三网格，则有：
Figure BDA0003102082360000051
Figure BDA0003102082360000052
Figure BDA0003102082360000053
Figure BDA0003102082360000054
Figure BDA0003102082360000055
其中，
Figure BDA0003102082360000056
为nΔt时刻顶点(m-1,j)处的波场值，
Figure BDA0003102082360000057
为nΔt时刻顶点(m,j)处的波场值，
Figure BDA0003102082360000058
为nΔt时刻顶点(m+2,j)处的波场值，
Figure BDA0003102082360000059
为nΔt时刻顶点(m-1,j)处的波场值，Δt为设定的时间步长，n为设定的步长倍数，Δx为x方向的设定步长，η1～η4为待定系数，c1为第一网格与第二网格的长度比值，c2为第三网格与第二网格的长度比值；
同时，所述空间差分算子
Figure BDA00031020823600000510
的计算公式为：
Figure BDA00031020823600000511
其中，Δz为z方向的设定步长；
上交界面处的网格顶点的时间差分算子
Figure BDA00031020823600000512
的计算方法为：
Figure BDA00031020823600000513
其中，
Figure BDA00031020823600000514
为
Figure BDA00031020823600000515
时刻网格顶点(i,j)处的波场值，
Figure BDA00031020823600000516
为
Figure BDA00031020823600000517
时刻网格顶点(i,j)处的波场值。
进一步地，步骤S4中将观测回波信号S1m与实测回波信号dobs分别进行三次以上的小波包分解时，每次只对观测回波信号S1m与实测回波信号dobs的低频部分进行小波包分解。
进一步地，步骤S6中按照设定步长减小三维网格模型中待测粘接层的厚度时，所述设定步长为待测粘接层划分的网格长度。
有益效果：
1、本发明提供一种基于变步长网格模型小波包分量波形的厚度遍历反演法，针对多层粘接结构不同层大差别区域物理参数，采用变步长网格模型，在保证反演精度的前提下大幅降低了模型规模；通过小波包分解方法对实测全波形和正演观测空间信号时域全波形进行处理，采用M层分解小波包分量波形即可将时域全波形数据量缩减到原始信号的1/2M，在指数级减少波形数据量的同时滤除了高频噪声信号；利用中间层厚度区间极小的特点，提出厚度遍历反演方法；通过上述三项措施，实现了中间薄层的高效高精度反演测量。
2、本发明提供一种基于变步长网格模型小波包分量波形的厚度遍历反演法，选择L1范数作为厚度反演的目标函数，鲁棒性更强。
附图说明
图1为本发明提供一种基于变步长网格模型小波包分量波形的厚度遍历反演法的流程图；
图2为本发明提供的变步长网格介质交界处两种离散情况示意图；
图3为本发明提供的变步长交错网格有限差分模型示意图；
图4为本发明提供的多层结构正演模型示意图；
图5为本发明提供的初始时域全波形信号示意图；
图6为本发明提供的小波包分解示意图；
图7(a)为本发明提供的粘接层厚度为0.4mm情况下的多层结构超声波信号示意图；
图7(b)为本发明提供的粘接层厚度为0.3mm情况下的多层结构超声波信号示意图；
图7(c)为本发明提供的粘接层厚度为0.2mm情况下的多层结构超声波信号示意图；
图8为本发明提供的全求解域遍历厚度反演与目标函数的变化示意图。
具体实施方式
为了使本技术领域的人员更好地理解本申请方案，下面将结合本申请实施例中的附图，对本申请实施例中的技术方案进行清楚、完整地描述。
实施例一
如图1所示，一种基于变步长网格模型小波包分量波形的厚度遍历反演法，待测粘接层位于上层结构和下层结构之间，形成三层结构，待测粘接层的厚度遍历反演法包括以下步骤：
S1：将激励信号从三层结构的上表面垂直入射，得到三层结构的实测回波信号dobs。
需要说明的是，在向三层结构入射激励信号时，先清洁三层结构，并在上层结构后涂抹耦合剂，再选用具有合适中心频率的超声纵波探头激励超声脉冲，从三层结构的上表面垂直入射，得到三层结构的实测回波信号dobs。
S2：构建三层结构的三维网格模型，设定三层结构各层的纵波声速、横波声速、密度以及厚度，同时设定三维网格模型各网格节点处的速度场和应力场的初始值，其中，待测粘接层的网格步长小于上层结构和下层结构的网格步长，且待测粘接层的厚度为可能达到的厚度极大值，网格节点包括网格的顶点以及网格边界的中点。
需要说明的是，在三层结构中，由于上层和下层介质层厚度较大，且与待测位置的物性参数无关，在满足收敛性条件下选择较为稀疏的网格。粘接层为待测区域且厚度较薄，为了提高厚度反演精度，选取细密的网格。三层结构各层的纵波声速、横波声速、密度以及厚度不完全相同，应力场、速度场的初始条件均为零。
S3：将激励模拟信号从三维网格模型的上表面垂直入射，得到三层结构的正演回波信号dcal，并从正演回波信号dcal提取上层结构与待测粘接层之间的界面的观测回波信号S1m。
S4：将观测回波信号S1m与实测回波信号dobs分别进行三次以上的小波包分解，对应得到小波包分量观测波形与小波包分量实测波形。
需要说明的是，在进行小波分解时，每次只对观测回波信号S1m与实测回波信号dobs的低频部分进行小波包分解，保证多次小波包分解后的低频部分保持了观测信号的主要信息，并且过滤了高频噪声信号。
S5：将小波包分量观测波形与小波包分量实测波形之间残差的L1范数作为目标函数值，判断所述目标函数值是否小于设定阈值或迭代次数是否达到设定上限，若满足其中一个条件，则待测粘接层当前设定的厚度为所求的厚度值，结束迭代；若均不满足，则进入步骤S6。
也就是说，本步骤是利用小波包分量波形构造目标函数来衡量正演的回波信号S1m与实际波形dobs之间的残差来反映正演模型与实际模型之间的差异，用于后续反演迭代过程。
在二维笛卡尔坐标系下的L1范数目标函数的表达式如下：
fL1＝∑|Dcal(x,z,t)-Dobs(x,z,t)|
其中，Dcal(x,z,t)表示小波包分量观测波形，Dobs(x,z,t)表示小波包分量实测波形。
S6：按照设定步长，如待测粘接层的网格厚度减小三维网格模型中待测粘接层的厚度，并重新划分网格和设定三维网格模型各网格节点处的速度场和应力场，实现三维网格模型的更新，再采用更新后的三维网格模型重新执行步骤S3～S5，直到得到待测粘接层的厚度值。
也就是说，本发明每次遍历更新的步长为单位粘接层网格厚度，厚度从所预设的厚度(厚度极大值)初值开始，每次反演迭代减少单位粘接层网格厚度，对二维变步长交错网格系统的速度场及压力场进行参数更新，获得更新后的正演模型——知道各网格节点处的速度场和应力场后的三维网格模型，以线性规律迭代直到达到终止条件，输出粘接层厚度值，结束迭代；否则重复S3～S5，对粘接层重新划分网格并更新差分算子的线性系数，再次计算目标函数值。
需要说明的是，在本发明提供的反演方法中，三层结构中间薄层厚度为待测物理量，目标函数为正演时域波形与实际模型实验所得时域波形残差的L1范数。通过正演计算可以得到正演模型对应的观测空间，结合实际模型数据空间求残差后计算目标函数值，当目标函数值收敛至足够小或达到最大迭代次数时终止算法得到最终正演模型，否则根据遍历法求得正演模型的修正量，更新模型进行新一轮正演。
进一步地，以三层结构厚度方向为笛卡尔坐标系x方向，以三层结构长度方向为笛卡尔坐标系z方向，所述速度场包括x方向的速度分量与z方向的速度分量，所述应力场包括切应力、x方向的正应力以及z方向的正应力，步骤S6中所述三维网格模型各网格节点处的速度场和应力场的设定方法为：
Figure BDA0003102082360000101
Figure BDA0003102082360000102
Figure BDA0003102082360000103
Figure BDA0003102082360000104
Figure BDA0003102082360000105
其中，i为网格顶点在x方向上的编号，j为网格顶点在z方向的编号，
Figure BDA0003102082360000106
为坐标为(i,j)的网格节点处x方向的空间差分算子，
Figure BDA0003102082360000107
为坐标为(i,j)的网格节点处z方向的空间差分算子，
Figure BDA0003102082360000108
为坐标为(i,j)的网格节点处的时间差分算子，ρi,j+1/2为坐标为(i,j+1/2)的网格节点处的密度，
Figure BDA0003102082360000109
为坐标为(i,j+1/2)的网格节点处x方向的速度分量，
Figure BDA00031020823600001010
为坐标为(i,j+1/2)的网格节点处x方向的正应力，
Figure BDA00031020823600001011
为坐标为(i,j+1/2)的网格节点处的切应力，ρi+1/2,j为坐标为(i+1/2,j)的网格节点处的密度，
Figure BDA00031020823600001012
为坐标为(i+1/2,j)的网格节点处z方向的速度分量，
Figure BDA00031020823600001013
为坐标为(i+1/2,j)的网格节点处z方向的正应力，
Figure BDA00031020823600001014
为坐标为(i+1/2,j)的网格节点处的切应力，
Figure BDA00031020823600001015
为坐标为(i+1/2,j+1/2)的网格节点处x方向的正应力，
Figure BDA00031020823600001016
为坐标为(i+1/2,j+1/2)的网格节点处x方向的速度分量，λi+1/2,j+1/2为坐标为(i+1/2,j+1/2)的网格节点处的一阶拉梅常数，μi+1/2,j+1/2为坐标为(i+1/2,j+1/2)的网格节点处的二阶拉梅常数，
Figure BDA00031020823600001017
为坐标为(i+1/2,j+1/2)的网格节点处z方向的速度分量，
Figure BDA00031020823600001018
为坐标为(i+1/2,j+1/2)的网格节点处z方向的正应力，
Figure BDA00031020823600001019
为坐标为(i,j)的网格节点处的切应力，μi,j为坐标为(i,j)的网格节点处的二阶拉梅常数，
Figure BDA0003102082360000111
为坐标为(i,j)的网格节点处x方向的速度分量，
Figure BDA0003102082360000112
为坐标为(i,j)的网格节点处z方向的速度分量。
需要说明的是，拉梅常数有一阶和二阶两个，一阶拉梅常数λ表示材料的压缩性，等价与体弹性模量或者Young氏模量，二阶拉梅常数μ表示材料的剪切模量，大致与G相当。也就是说，两个拉梅常数只和材料有关，并不是换一个节点就变了一个值，待测粘接层、上层结构以及下层结构的一阶拉梅常数λ和二阶拉梅常数μ均不同，只是为了表述统一，速度场和应力场的计算公式中均加上网格的坐标。
需要说明的是，i，j为网格顶点，但是网格节点除了包括网格顶点，还包括网格边界的中点，因此用i+1/2，j+1/2来表示网格边界的中点，以上也统一使用节点来进行坐标位置的解释。
此外，由以上速度场和应力场的计算公式可知，每过一个单位时间，可以从应力场推导得到相邻网格节点的速度场，再过一个单位时间，可以从速度场推导到相邻节点的应力场，从而对正演过程中应力场和速度场的参数更新，实现每次迭代过程中的正演模型更新。
进一步地，空间差分算子
Figure BDA0003102082360000113
空间差分算子
Figure BDA0003102082360000114
以及时间差分算子
Figure BDA0003102082360000115
的获取方法为：将激励模拟信号在各网格顶点处的波动进行泰勒展开，得到各网格顶点处的波场值，再根据波场值获取各网格顶点处的空间差分算子
Figure BDA0003102082360000116
空间差分算子
Figure BDA0003102082360000117
以及时间差分算子
Figure BDA0003102082360000118
具体的，上层结构与待测粘接层之间的界面为上交界面，则上交界面处的网格顶点的空间差分算子
Figure BDA0003102082360000119
和
Figure BDA00031020823600001110
的计算方法为：
依次在x方向上取四个网格顶点执行空间差分算子获取操作，且每次选出的四个网格顶点所在的三个网格单元中，至少有一个、至多有两个属于待测粘接层，如图2所示，圆形表示属于上层结构的网格顶点，三角形表示属于待测粘接层的网格顶点；其中，所述空间差分算子获取操作为：
将四个网格顶点沿上层结构至待测粘接层的方向分别记为(m-1,j)、(m,j)、(m+1,j)以及(m+2,j)，顶点(m-1,j)和顶点(m,j)所在的网格记为第一网格，顶点(m,j)和顶点(m+1,j)所在的网格记为第二网格，顶点(m+1,j)和顶点(m+2,j)所在的网格记为第三网格，则有：
Figure BDA0003102082360000121
Figure BDA0003102082360000122
Figure BDA0003102082360000123
Figure BDA0003102082360000124
Figure BDA0003102082360000125
其中，
Figure BDA0003102082360000126
为nΔt时刻顶点(m-1,j)处的波场值，
Figure BDA0003102082360000127
为nΔt时刻顶点(m,j)处的波场值，
Figure BDA0003102082360000128
为nΔt时刻顶点(m+2,j)处的波场值，
Figure BDA0003102082360000129
为nΔt时刻顶点(m-1,j)处的波场值，Δt为设定的时间步长，n为设定的步长倍数，Δx为x方向的设定步长，η1～η4为待定系数，c1为第一网格与第二网格的长度比值，c2为第三网格与第二网格的长度比值。
需要说明的是，若选出的四个网格顶点所在的三个网格单元中只有一个属于待测粘接层，则第二网格和第三网格之间的交界面为上交界面；若选出的四个网格顶点所在的三个网格单元中有两个属于待测粘接层，则第一网格和第二网格之间的交界面为上交界面。
进一步地，如图3所示，假设节点m和m+1的中点为节点i顶点，(m-1,j)，(m,j)，(m+1,j)和(m+2,j)处的波动以x方向进行泰勒展开如下：
Figure BDA0003102082360000131
Figure BDA0003102082360000132
Figure BDA0003102082360000133
Figure BDA0003102082360000134
需要说明的是，为了简洁，以上泰勒展开省略了上标n。将Δ以Δx替换，由上述泰勒展开推导可得变步长交错网格系统对于x方向的空间差分算子：
同时，所述空间差分算子
Figure BDA0003102082360000135
的计算公式为：
Figure BDA0003102082360000136
其中，Δz为z方向的设定步长；
进一步地，顶点(m-1,j)，(m,j)，(m+1,j)和(m+2,j)处的波动以z方向进行泰勒展开如下：
Figure BDA0003102082360000137
Figure BDA0003102082360000138
Figure BDA0003102082360000141
Figure BDA0003102082360000142
则将Δ以Δz替换即可以得到：
Figure BDA0003102082360000143
需要说明的是，为了简洁，以上泰勒展开省略了上标n。同时，横向单位步长Δz在三层结构正演模型中是固定值，就是不改变的，这个左边的Δ可以写成Δz，但是它不一定和x方向的Δ一样大。
上交界面处的网格顶点的时间差分算子
Figure BDA0003102082360000144
的计算方法为：
Figure BDA0003102082360000145
其中，
Figure BDA0003102082360000146
为
Figure BDA0003102082360000147
时刻网格顶点(i,j)处的波场值，
Figure BDA0003102082360000148
为
Figure BDA0003102082360000149
时刻网格顶点(i,j)处的波场值。
实施例二
如图4所示，所述的三层结构的中间层位于上层结构和下层结构之间，形成三层结构，其中，本发明方法尤其适用于中间层的厚度远远小于上层结构和下层结构厚度，即薄层厚度与上层结构和下层结构厚度相差一个量级的三层结构；可选的，薄层为硅橡胶，上层结构为金属，下层结构为高分子材料，三层结构的物理特性参数如表1。
表1多层结构试样物性参数表
Figure BDA0003102082360000151
本发明使用由超声收发仪、示波器和探头组成超声检测系统，清洁三层结构的上层结构表面后涂抹耦合剂，使用5MHz纵波接触式探头测量三层结构的上表面中央位置，使超声脉冲垂直入射，得到三层结构的实验回波信号。
观测时域全波形信号如图5所示，将时域全波形信号进行小波包分解，可以得到高低频两部分信号。初始信号为20000个时域采样点，将其进行三次小波包分解，其中每次只对低频部分进行小波包分解。每次小波包分解得到的高频细节部分和三次小波包分解后得到的低频主体部分如图6所示。三次小波包分解后的低频部分保持了观测信号的主要信息，同时将采样点减少为初始观测信号的1/8，大大缩小了反演数据处理量并且滤除了高频噪声信号。
粘接层厚度测量一般通过脉冲回波法进行，通过计算界面回波声时差与声速的乘积即可。但在厚度较薄的情况下，界面回波发生混叠，使得回波信号无法被清晰地分辨出来。如图7(a)、图7(b)、图7(c)所示，分别表示粘接层厚度为0.4mm、0.3mm和0.4mm的声时图像，其中0～0.75μs为激励波，1.5μs～2.25μs和3μs～3.75μs分别为上层结构-粘接层界面的一次回波和二次回波，一次和二次回波之间的小波为粘接层-下层结构界面的回波。此时已经无法从时域上观察出硅橡胶-高分子材料的回波声时位置，脉冲回波法测厚已不再适用。
采用全波形反演方法对厚度小于等于0.3mm粘接层进行测量。在反演过程中，为了减少反演计算量，提高反演计算效率，使用变步长网格进行厚度反演。设置铝合金和高分子材料介质层网格步长为0.05mm，粘接层硅橡胶介质网格步长为0.025mm。设置厚度反演的初始值为0.5mm，对真实值为0.2mm厚粘接层进行反演，每次遍历之间的步长为单位粘接层网格厚度，即0.025mm，全求解域为20次反演迭代。厚度从0.5mm开始每次反演迭代减少0.025mm，终止条件为目标函数小于5。当满足终止条件时，迭代终止，输出中间薄层的厚度值。否则重复S3～S7，对粘接层重新划分网格并更新差分算子的线性系数，再次计算目标函数值。
全求解域遍历反演方法厚度与目标函数的变化如图8所示。厚度从0.5mm开始每次反演迭代减少0.025mm，以线性规律迭代13次达到真实值0.2mm。
至此，根据本实施例提出的基于变步长网格模型小波包分量波形的高效厚度遍历反演法，可以在已知三层结构各层纵波声速、横波声速和密度，及上、下层的厚度的情况下，对厚度小于等于0.3mm粘接层进行测量，通过建立变步长网格模型与利用小波包分量波形，指数级地减少了反演计算量，有效提高了反演计算效能。
当然，本发明还可有其他多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员当然可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。

Claims (6)

1.一种基于变步长网格模型小波包分量波形的厚度遍历反演法，待测粘接层位于上层结构和下层结构之间，形成三层结构，其特征在于，待测粘接层的厚度遍历反演法包括以下步骤：
S1：将激励信号从三层结构的上表面垂直入射，得到三层结构的实测回波信号dobs；
S2：构建三层结构的三维网格模型，设定三层结构各层的纵波声速、横波声速、密度以及厚度，同时设定三维网格模型各网格节点处的速度场和应力场的初始值，其中，待测粘接层的网格步长小于上层结构和下层结构的网格步长，且待测粘接层的厚度为可能达到的厚度极大值，网格节点包括网格的顶点以及网格边界的中点；
S3：将激励模拟信号从三维网格模型的上表面垂直入射，得到三层结构的正演回波信号dcal，并从正演回波信号dcal提取上层结构与待测粘接层之间的界面的观测回波信号S1m；
S4：将观测回波信号S1m与实测回波信号dobs分别进行三次以上的小波包分解，对应得到小波包分量观测波形与小波包分量实测波形；
S5：将小波包分量观测波形与小波包分量实测波形之间残差的L1范数作为目标函数值，判断所述目标函数值是否小于设定阈值或迭代次数是否达到设定上限，若满足其中一个条件，则待测粘接层当前设定的厚度为所求的厚度值，若均不满足，则进入步骤S6；
S6：按照设定步长减小三维网格模型中待测粘接层的厚度，并重新划分网格和设定三维网格模型各网格节点处的速度场和应力场，实现三维网格模型的更新，再采用更新后的三维网格模型重新执行步骤S3～S5，直到得到待测粘接层的厚度值。
2.如权利要求1所述的一种基于变步长网格模型小波包分量波形的厚度遍历反演法，其特征在于，以三层结构厚度方向为笛卡尔坐标系x方向，以三层结构长度方向为笛卡尔坐标系z方向，所述速度场包括x方向的速度分量与z方向的速度分量，所述应力场包括切应力、x方向的正应力以及z方向的正应力，步骤S6中所述三维网格模型各网格节点处的速度场和应力场的设定方法为：
Figure FDA0003102082350000021
Figure FDA0003102082350000022
Figure FDA0003102082350000023
Figure FDA0003102082350000024
Figure FDA0003102082350000025
其中，i为网格顶点在x方向上的编号，j为网格顶点在z方向的编号，
Figure FDA0003102082350000026
为坐标为(i,j)的网格节点处x方向的空间差分算子，
Figure FDA0003102082350000027
为坐标为(i,j)的网格节点处z方向的空间差分算子，
Figure FDA0003102082350000028
为坐标为(i,j)的网格节点处的时间差分算子，ρi,j+1/2为坐标为(i,j+1/2)的网格节点处的密度，
Figure FDA0003102082350000029
为坐标为(i,j+1/2)的网格节点处x方向的速度分量，
Figure FDA00031020823500000210
为坐标为(i,j+1/2)的网格节点处x方向的正应力，
Figure FDA00031020823500000211
为坐标为(i,j+1/2)的网格节点处的切应力，ρi+1/2,j为坐标为(i+1/2,j)的网格节点处的密度，
Figure FDA00031020823500000212
为坐标为(i+1/2,j)的网格节点处z方向的速度分量，
Figure FDA00031020823500000213
为坐标为(i+1/2,j)的网格节点处z方向的正应力，
Figure FDA00031020823500000214
为坐标为(i+1/2,j)的网格节点处的切应力，
Figure FDA00031020823500000215
为坐标为(i+1/2,j+1/2)的网格节点处x方向的正应力，
Figure FDA00031020823500000216
为坐标为(i+1/2,j+1/2)的网格节点处x方向的速度分量，λi+1/2,j+1/2为坐标为(i+1/2,j+1/2)的网格节点处的一阶拉梅常数，μi+1/2,j+1/2为坐标为(i+1/2,j+1/2)的网格节点处的二阶拉梅常数，
Figure FDA0003102082350000031
为坐标为(i+1/2,j+1/2)的网格节点处z方向的速度分量，
Figure FDA0003102082350000032
为坐标为(i+1/2,j+1/2)的网格节点处z方向的正应力，
Figure FDA0003102082350000033
为坐标为(i,j)的网格节点处的切应力，μi,j为坐标为(i,j)的网格节点处的二阶拉梅常数，
Figure FDA0003102082350000034
为坐标为(i,j)的网格节点处x方向的速度分量，
Figure FDA0003102082350000035
为坐标为(i,j)的网格节点处z方向的速度分量。
3.如权利要求2所述的一种基于变步长网格模型小波包分量波形的厚度遍历反演法，其特征在于，空间差分算子
Figure FDA0003102082350000036
空间差分算子
Figure FDA0003102082350000037
以及时间差分算子
Figure FDA0003102082350000038
的获取方法为：
将激励模拟信号在各网格顶点处的波动进行泰勒展开，得到各网格顶点处的波场值，再根据波场值获取各网格顶点处的空间差分算子
Figure FDA0003102082350000039
空间差分算子
Figure FDA00031020823500000310
以及时间差分算子
Figure FDA00031020823500000311
4.如权利要求3所述的一种基于变步长网格模型小波包分量波形的厚度遍历反演法，其特征在于，上层结构与待测粘接层之间的界面为上交界面，则上交界面处的网格顶点的空间差分算子
Figure FDA00031020823500000312
和
Figure FDA00031020823500000313
的计算方法为：
依次在x方向上取四个网格顶点执行空间差分算子获取操作，且每次选出的四个网格顶点所在的三个网格单元中，至少有一个、至多有两个属于待测粘接层；其中，所述空间差分算子获取操作为：
将四个网格顶点沿上层结构至待测粘接层的方向分别记为(m-1,j)、(m,j)、(m+1,j)以及(m+2,j)，顶点(m-1,j)和顶点(m,j)所在的网格记为第一网格，顶点(m,j)和顶点(m+1,j)所在的网格记为第二网格，顶点(m+1,j)和顶点(m+2,j)所在的网格记为第三网格，则有：
Figure FDA0003102082350000041
Figure FDA0003102082350000042
Figure FDA0003102082350000043
Figure FDA0003102082350000044
Figure FDA0003102082350000045
其中，
Figure FDA0003102082350000046
为nΔt时刻顶点(m-1,j)处的波场值，
Figure FDA0003102082350000047
为nΔt时刻顶点(m,j)处的波场值，
Figure FDA0003102082350000048
为nΔt时刻顶点(m+2,j)处的波场值，
Figure FDA0003102082350000049
为nΔt时刻顶点(m-1,j)处的波场值，Δt为设定的时间步长，n为设定的步长倍数，Δx为x方向的设定步长，η1～η4为待定系数，c1为第一网格与第二网格的长度比值，c2为第三网格与第二网格的长度比值；
同时，所述空间差分算子
Figure FDA00031020823500000410
的计算公式为：
Figure FDA00031020823500000411
其中，Δz为z方向的设定步长；
上交界面处的网格顶点的时间差分算子
Figure FDA00031020823500000412
的计算方法为：
Figure FDA00031020823500000413
其中，
Figure FDA00031020823500000414
为
Figure FDA00031020823500000415
时刻网格顶点(i,j)处的波场值，
Figure FDA00031020823500000416
为
Figure FDA00031020823500000417
时刻网格顶点(i,j)处的波场值。
5.如权利要求1～4任一权利要求所述的一种基于变步长网格模型小波包分量波形的厚度遍历反演法，其特征在于，步骤S4中将观测回波信号S1m与实测回波信号dobs分别进行三次以上的小波包分解时，每次只对观测回波信号S1m与实测回波信号dobs的低频部分进行小波包分解。
6.如权利要求1～4任一权利要求所述的一种基于变步长网格模型小波包分量波形的厚度遍历反演法，其特征在于，步骤S6中按照设定步长减小三维网格模型中待测粘接层的厚度时，所述设定步长为待测粘接层划分的网格长度。
CN202110629436.7A 2021-06-04 2021-06-04 基于变步长网格模型小波包分量波形的厚度遍历反演法 Active CN113376256B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110629436.7A CN113376256B (zh) 2021-06-04 2021-06-04 基于变步长网格模型小波包分量波形的厚度遍历反演法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110629436.7A CN113376256B (zh) 2021-06-04 2021-06-04 基于变步长网格模型小波包分量波形的厚度遍历反演法

Publications (2)

Publication Number Publication Date
CN113376256A true CN113376256A (zh) 2021-09-10
CN113376256B CN113376256B (zh) 2023-03-17

Family

ID=77576227

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110629436.7A Active CN113376256B (zh) 2021-06-04 2021-06-04 基于变步长网格模型小波包分量波形的厚度遍历反演法

Country Status (1)

Country Link
CN (1) CN113376256B (zh)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115952625A (zh) * 2023-03-15 2023-04-11 南京航空航天大学 一种行波旋转型超声电机动力学仿真方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140372043A1 (en) * 2013-06-17 2014-12-18 Wenyi Hu Full Waveform Inversion Using Perfectly Reflectionless Subgridding
CN104977607A (zh) * 2014-04-09 2015-10-14 中国石油集团东方地球物理勘探有限责任公司 利用变步长网格声波波场模拟的时间域全波形反演方法
CN108051855A (zh) * 2017-12-13 2018-05-18 国家深海基地管理中心 一种基于拟空间域声波方程的有限差分计算方法
CN109188519A (zh) * 2018-09-17 2019-01-11 中国石油大学(华东) 一种极坐标下的弹性波纵横波速度反演系统及方法
CN110308204A (zh) * 2019-07-05 2019-10-08 北京理工大学 一种三层结构中间薄层物理特性参数测量方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20140372043A1 (en) * 2013-06-17 2014-12-18 Wenyi Hu Full Waveform Inversion Using Perfectly Reflectionless Subgridding
CN104977607A (zh) * 2014-04-09 2015-10-14 中国石油集团东方地球物理勘探有限责任公司 利用变步长网格声波波场模拟的时间域全波形反演方法
CN108051855A (zh) * 2017-12-13 2018-05-18 国家深海基地管理中心 一种基于拟空间域声波方程的有限差分计算方法
CN109188519A (zh) * 2018-09-17 2019-01-11 中国石油大学(华东) 一种极坐标下的弹性波纵横波速度反演系统及方法
CN110308204A (zh) * 2019-07-05 2019-10-08 北京理工大学 一种三层结构中间薄层物理特性参数测量方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
张文生等: "时间域两网格全波形反演", 《数值计算与计算机应用》 *
曲英铭等: "基于多尺度双变网格的时间域全波形反演", 《石油物探》 *
胡光辉等: "三维声波全波形反演的实现与验证", 《石油物探》 *
魏妮娜: "变网格有限差分波场数值模拟的研究", 《中国优秀硕士学位论文全文数据库》 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115952625A (zh) * 2023-03-15 2023-04-11 南京航空航天大学 一种行波旋转型超声电机动力学仿真方法

Also Published As

Publication number Publication date
CN113376256B (zh) 2023-03-17

Similar Documents

Publication Publication Date Title
Zhao et al. Ultrasonic Lamb wave tomography in structural health monitoring
Malyarenko et al. Ultrasonic Lamb wave diffraction tomography
Lemistre et al. Structural health monitoring system based on diffracted Lamb wave analysis by multiresolution processing
CN110308204B (zh) 一种三层结构中间薄层物理特性参数测量方法
CN104688224B (zh) 一种应用于声学非均匀媒介磁声耦合成像重建方法
CN110286155B (zh) 一种多层复合材料的损伤检测方法及系统
CN113376256B (zh) 基于变步长网格模型小波包分量波形的厚度遍历反演法
Rahani et al. Modeling of transient ultrasonic wave propagation using the distributed point
Yang et al. Parametric study on interply tracking in multilayer composites by analytic-signal technology
Hameed et al. Transverse damage localization and quantitative size estimation for composite laminates based on Lamb waves
Xue et al. Damage localization and robust diagnostics in guided-wave testing using multitask complex hierarchical sparse Bayesian learning
Singh et al. The role of digital signal processing in NDT
Karagiannis et al. Three-dimensional nondestructive “sampling” of art objects using acoustic microscopy and time–frequency analysis
Zhang et al. Electromagnetic ultrasonic signal processing and imaging for debonding detection of bonded structures
CN117233266A (zh) 一种基于循环神经网络的全波形反演导波层析成像方法
Wang et al. Health monitoring of plate structures based on tomography with combination of guided wave transmission and reflection
CN113376257B (zh) 一种基于共轭梯度优化方法反演的多参数测量方法
Diligent et al. Reflection and scattering of the S 0 Lamb mode from circular defects in plates
Imperiale et al. Smart numerical tools for the modelling of ultrasonic testing on curved composite structures
Huang et al. Theory and methodology of electromagnetic ultrasonic guided wave imaging
Chiou et al. Modeling of ultrasonic signals from weak inclusions
Gust Verbesserung der Signalauswertung für die Ultraschallmikroskopie
Wang et al. Full Waveform Inversion Guided Wave Tomography Based on Recurrent Neural Network
Chen et al. Sequential hypothesis testing for reflected signal recognition of time-of-flight estimation in ultrasonic defect detection
ANWAR et al. Ultrasonic Multi-Hole Imaging Using Full Waveform Inversion

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant