基于双平方根算子叠前深度偏移成像技术的结构损伤成像与探测

郭原草，郭少华

(中南大学 土木工程学院，湖南 长沙，410004)

摘要：介绍一种基于双平方根算子的叠前深度偏移方法用于结构损伤探测。首先推导上行波传播时间(走时)的全微分表达式，推导在频率波数域和频率空间域中，用源点和接收点坐标表示的双平方根算子波动方程。然后，用频率空间域的有限差分法及混合法进行双平方根算子的波场外推，推导关于炮点和检波点的差分方程，以及双平方根算子混合法波场外推公式。最后，用Marmousi数据体对结构损伤进行偏移成像。根据双程波动方程分析该方法的数值频散性、奇异性以及波场、位移、振幅、能量、动应力集中系数等力学特征，计算得到波场快照图和剖面图，达到对结构损伤成像的目的。数值算例证明了这种方法的适用性。

关键词：叠前深度偏移；结构损伤成像；双平方根算子；波场外推

中图分类号：TG111.3        文献标志码：A         文章编号：1672-7207(2013)05-2022-09

Structural damage imaging and detection based on prestack depth migration imaging technology of double square root operator

GUO Yuancao, GUO Shaohua

(School of Civil Engineering, Central South University, Changsha 410004, China)

Abstract: A prestack depth migration based on the double square root operator was introduced. Firstly, the total difference expression of wave propagation time was deduced, then, the double square root operator wave equation expressed by the source and receiver coordinate was deduced in two spaces, respectively. Then, the double square root operator wave field extrapolation was executed by the finite difference and mixed method. The finite difference equations about the explosion and receiver were deduced. At last, the migration image was executed based on the Marmousi data.The numerical volatilization, singularity, wave field, displacement, amplitude, energy, dynamic stress concentration coefficiency according to the double distance wave equation, the wave field snapshot and section were ultimately obtained. Numerical examples testify the feasibility of the method.

Key words: prestack depth migration; structural damage image; double square root operator; wave field extrapolation

目前，结构损伤探测领域主要采用基于结构动力学理论的检测方法。这种方法主要分析结构振动中的动力学模态参数的变化。曹晖等^[1]采用单元刚度折减系数应用于简支梁损伤并找出了损伤位置。邹万杰等^[2]采用模型修正法求出刚度联系矩阵，然后，计算刚度损伤量。郭惠勇等^[3]采用频率和位移进行结构二次损伤识别。王根会等^[4]采用单元模态应变能变化作为桥梁损伤定位指标。李学平等^[5]引入损伤分布函数进行混凝土简支梁损伤识别。赵启林等^[6]采用静态信息结构的模式匹配法对桁架和超静定梁进行损伤定位。冉志红等^[7]采用模态柔度法对连续梁桥进行损伤定位研究。焦莉等^[8]应用小波分析对5层剪力墙结构进行结构损伤诊断。Lee等^[9]利用压电陶瓷技术结合机电阻抗法进行结构损伤诊断。叠前偏移成像技术由于其直观性已经从最初的地球物理勘探领域延伸到工程结构损伤识别领域。Hudson^[10]分析了损伤固体中弹性波高阶近似。牛滨华等^[11]分析了裂隙各向异性介质波动分量记录的数值模拟。等^[12]研究了带有平行线性裂缝介质中的弹性波传播。Lin等^[13]首先将叠前逆时偏移技术应用到铝质板的多重损伤探测中。Arkadiusz 等^[14]利用面波对带有裂缝损伤的各向异性板式结构进行损伤成像。张文生等^[15]分析了二维横向各向同性介质的伪谱法模拟。孙卫涛^[16]分析了弹性波动方程有限差分方法在复杂介质探测中的应用。目前，工业上普遍使用的叠前深度偏移方法是Kirchhoff积分法。该方法具有计算速度快、对观测系统适应性强的特点。但是对于复杂的结构损伤，采用积分法难以进行精确成像。炮集叠前偏移法是精度最高的叠前深度成像方法之一。本研究介绍基于双平方根算子的叠前偏移方法。

1  模型建立

考虑1个沿x轴方向传播的平面波，传播方向与轴夹角为θ，传播速度为v，角频率为ω，波场位移为U，经过距离dx后所用时间为dt，则

                   (1)

若考虑上行波，则应满足：

          (2)

            (3)

在频率波数域中利用对应关系可得：

             (4)

其中：i为虚数单位；k_x为x方向的波数；z_r和z_s分别为接收点与震源z方向坐标；x_r和x_s分别为接收点与震源的x方向坐标。当接收点下移距离dz_r时，上行波走时变化满足

            (5)

当震源下移距离dz_s，则上行波走时变化满足

           (6)

当源点和接受点同时下移dz=dz_r=dz_s时，则上行波走时变化应为

            (7)

或者

   (8)

将式(8)返回频率波数域，这时，U(x_r, x_s, t)对应U(k_r, k_s, ω)。

   (9)

其中：k_r和k_s分别为与x_r和x_s对应的波数。在频率空间域中可以表示为：

 (10)

其中：v_r=v(x_r, z)和v_s=v(x_s, z)分别表示接收点与震源处的波速。式(9)和(10)是以源点和接受点坐标表示的双平方根算子方程。双平方根算子方程也可以用中心点x_m和半炮检距坐标x_h来表示。

    (11)

    (12)

频率空间域中以中心点和半炮检距表示的双程波波动方程为

       (13)

同样可以得到相应的频率波数域中的双程波波动方程，将U(x_s, x_r, z, ω)看成U(x_m, x_h, z, ω)可得：

    (14)

    (15)

将式(14)和(15)变回到波数域，得到频率波数域以中心点波数k_m和半炮检距波数k_h表示的双平方根算子方程：

         (16)

其中：k_m和k_n分别对应x_m和x_n的波数，v_s=v(x_s, z)和v_r=v(x_r, z)分别是源点和接收点移动方向波速；δ_s和δ_r分别为炮点与检波点差分算子计算步长。双平方根算子的波场外推法可以用频率空间域的有限差分及混合法。既可以在炮点检波点域中操作，也可以在中心点半炮检距域中操作。在计算中记录U(x_s, x_r, z=0, ω)或U(x_m,x_h,z=0,ω)是已知的。成像原理是。设△s为炮点间距，△r为检波点间距，a和b为系数，取a=1/2，b=1/4。因此，计算只要对记录进行外推即可。这里先考虑频率空间域中的有限差分法，根据单平方根算子近似，考虑炮点的差分方程为

      (17)

关于检波点的差分方程为

       (18)

在中心点-半炮检距域中，频率空间域中所求解的方程为

          (19)

类似单平方根算子的混合法，对双平方根算子，炮点-检波点域中的混合法波场外推公式为

+        (20)

其中：v₀为波场参考速度。在中心点-半偏移距域中，对应的混合法波场外推公式为：

    (21)

2  数值算例

2.1 波场位移与振幅变化

建立结构模型：长×宽为200 cm×200 cm的铝质板式结构，密度ρ=2 710 kg/m³，厚度h=1.6 cm，泊松比=0.3，弹性模量E=72.5 GPa。结构内部损伤为常速介质中的1个水平裂缝(界面)。裂缝长度×宽度为15 mm×1 mm。介质速度设为v=3 200 m/s采样率为△x=△z=1 cm，△t=2 ms，时间样点数为250，初始激励幅值P=1 N，f₀=40 kHz。假设板中只存在裂缝。波场激励和接受装置采用压电薄片，在铝板上沿着裂缝长度按照△x=1 cm布置9个。通过人工激励方法模拟得到损伤位置(10 cm, -41 cm)，(-25 cm, -25 cm)。根据式(10)计算波场位移变化规律。这里取如下参数△t=2 ms，f=30 Hz，k₁=7，k₂=9，k₁和k₂分别为x和z方向的波数。由式(18)可以计算出波场振幅和位移的函数关系，如图1和图2所示。从图1和图2中看出：由于结构中损伤的存在对弹性波场起到阻碍和延缓作用，使位移和振幅随着传播距离和时间的增加逐渐减小，减小的速率也逐渐减小。这种传播特征反映出波场能量在结构损伤中逐渐消耗。

图1  波场振幅随距离变化曲线

Fig.1  Relationship between wave field amplitude variation and horizontal coordinate x

图2  波场位移随距离变化曲线

Fig.2  Relationship between wave field displacement variation and horizontal coordinate x

2.2  波场能量变化

为了进一步研究弹性波在工程结构损伤区域中的传播能量变化规律，设n为传播方向，l，m和n为方向余弦，r为波场任意点的矢径，λ和μ为介质弹性常数，沿水平方向传播的平面简谐纵波的位移满足：

           (22)

其中：；；A为振幅，能量密度。其中应变能密度P_u为

     (23)

动能密度P_k为

       (24)

即

    (25)

于是能量密度ε为

   (26)

能量密度矢量I(能通量密度矢量)可以表示为

   (27)

即

    (28)

取λ=2，μ=0.4，K=6，cosθ=0.4，cosα=0.5根据式(28)和(30)作水平方向应力变化曲线与传播能量变化曲线，如图3和图4所示。从图3和图4可以看出：应力变化规律近似于位移变化规律，而能量传播随着损伤的出现而逐渐消耗并减少。

2.3  动应力集中程度

结构内部损伤的存在将引起动应力集中现象，其严重程度采用动应力集中系数τ表示。它是弹性波在一点上产生的动应力与入射波在同一点上产生的动应力之比。W_I和W_II分别为反射波与入射波位移，裂缝周边的动应力分布系数τ_θz^(G)可以由下式计算：

图3  波场水平方向应力随时间变化曲线

Fig.3  Relationship between wave field horizon stress variation and propagation time t

图4  波场能量密度随距离变化曲线

Fig.4  Relationship between wave field energy density variation and horizon distance x

           (29)

  (30)

 (31)

其中：A和B为待定任意常数；H_n为第一类Hankel函数。入射波引起的切应力最大幅值τ₀为

                 (32)

动应力集中系数为

                (33)

在计算中采用如下映射形式：

            (34)

这里进行如下数学处理：将映射平面上的单位圆映射成复平面上的椭圆。设z，η和m为波场应力复平面参量，k₁和k₂分别为x和z方向的波数，h和b分别为裂缝深度与宽度，k₁a=1，h/b=1.5，k₂b=1.2。根据式(21)作不同波数情况下动应力集中系数极坐标变化图。图5所示为裂缝周围动应力集中系数的分布规律。裂缝的存在对于动应力系数分布规律有很大的影响。当损伤尺寸较小时，这种影响尤其明显。动应力集中系数的最大值和整个分布形状有向损伤靠近的趋势，就像被裂缝吸引一样。

图5  动应力集中系数分布变化曲线

Fig.5  Distribution variation curve of dynamic concentration coefficient

2.4  数值频散和稳定性

上述差分格式对于网格具有时间上的二阶精度和空间上的四阶精度。为了便于分析数值各向异性，引入孔隙度的概念，用有限个孔隙的离散组合近似代替裂缝损伤。为孔隙度，ζ=ρ_f/ρ_s为密度比。引入控制数值频散的无量纲参数ξ=v^P-solid△t/△和网格尺寸无量纲参数H=△/λ；γ₁，γ₂和γ₃为波传播方向与坐标轴夹角，η₁和η₂为四阶差分系数。为便于分析，将损伤结构看作双相各向异性介质。下面给出双相各向异性介质的数值频散和稳定性分析。当结构内部骨架体积模量K_b、剪切模量μ与流体体积模量K_f之间满足K_b，μ＞＞ K_f时(结构内部空洞所含空气看作流体)，快纵波数值频散推导成如下形式^[25]：

   (35)

  (36)

很明显，数值频散无量纲量q与泊松比无关，而且式中出现了孔隙度和孔隙度几何度量的系数项。同理可得慢纵波数值频散和横波数值频散计算式：

   (37)

(38)

其中：q_p-fast，q_p-slow和q_s分别为快纵波、慢纵波和横波的数值频散无量纲量；v_f，v_p-solid和v_shear-solid分别为流体介质中纵波波速、固体介质中纵波波速与固体介质中剪切波波速。

这里考察2个不同传播方向(即γ₁=45°，γ₂=45°，γ₃=90°和γ₁=54.7°，γ₂=54.7°，γ₃=54.7°)，ξ取0.6，孔隙流体与骨架固体介质密度比为0.426，孔隙度为0.13，几何度量参数为a=3.5，骨架固体纵波速度为3 800 m/s，横波速度为2 400 m/s，孔隙流体纵波速度为1 600 m/s，四阶差分系数为η₁=9/8，η₂=-1/20。由于差分网格离散引起的快纵波、横波和慢纵波的数值频散曲线如图6和图7所示。其中横坐标为单位波长内网格点数控制参数，纵坐标为数值频散参数。考察2个不同方向的频散效应。从图6和图7可以看出：波场数值计算速度因传播方向不同而大于或小于真实速度。本研究中交错网格四阶差分格式的数值频散效应比二阶差分格式弱得多。在相同计算精度下，对网格尺寸的要求大大放宽，单位波长内网格点数从至少10个减少到5个。

图6  二阶有限差分快纵波数值频散曲线

Fig.6  Numerical frequency volatilization curve of second order finite difference rapid longitude wave

图7  四阶有限差分横波数值频散曲线

Fig.7  Numerical frequency volatilization curve of fourth order finite difference rapid transverse wave

以上计算表明：用全波波动方程和无层间反射的波动方程可以合成零偏移距记录。对于边界吸收，采用反周期扩展法进行有效处理。模型计算结果表明了方法及算法的正确性和有效性。傅里叶变换采用快速离散傅里叶算法实现，提高了计算速度。

弹性波在离散网格点上沿着不同方向传播的不同频散特性会导致数值各向异性现象出现。弹性波各向异性参数使用无量纲量A表示：

          (39)

其中：v_num为数值计算波场速度；v_true为真实波场速度。图8~10所示为不同孔隙度的纵波和横波各向异性曲面，网格参数H=0.125。网格离散引起的各向异性不仅沿不同传播方向变化，且随着孔隙度增大而增大；当H=0.125时，单位波长网格点是8个，本研究得到的快纵波各向异性最大值不超过1%；即使在孔隙度为0.45时也不超过1.4%；2种孔隙度下的横波各向异性都小于0.18%，慢纵波的各向异性均小于0.14。由于损伤与孔隙的存在，导致数值频散加剧。在相同计算精度下，需要增加单位波长内的网格点数，但是，与二阶差分格式相比减少了对网格点数的要求，提高了正演模拟的速度。

图8  快纵波数值各向异性随方向变化曲面

Fig.8  Rapid longitude wave numerical anisotropy surface with direction variation

图9  慢纵波数值各向异性随方向变化曲面

Fig.9  Slow longitude wave numerical anisotropy surface with direction variation

图10  横波数值各向异性随方向变化曲面

Fig.10  Transverse wave numerical anisotropy surface with direction variation

2.5  波场快照与损伤成像

上面推导出了2个坐标空间中的用不同坐标表示的双平方根算子方程以及使用2种不同方法在2个空间域中的波场外推公式。下面建立如下结构模型：模型是1个铝质板式结构，密度ρ=2 710 kg/m，厚度h=0.15 cm，泊松比=0.34，弹性模量E=750 GPa。结构内部损伤由常速介质中的1个夹杂(点散射体)和1个水平裂缝(界面)组成。数据体采用Marmousi模型数据，介质速度设定为v=2 500 m/s，采样率为△x=△z=20 m，△t=20 ms，时间样点数为250。图11所示为在频率空间域内用有限差分法计算得出的t=0.48 ms时刻波场水平应力叠前偏移记录，该记录由1条双曲线和1条水平曲线构成。图12所示为在频率波数域内用混合法计算得到的t=0.96 ms时刻波场水平应力叠前偏移记录。对应的叠前偏移记录在中间出现了与结构损伤形状相似的水平阴影以及和波场形状相似的弯曲。通过比较可以看到：通过双平方根算子叠前深度偏移，由点散射体产生的散射波得到了较好的成像^[13]。

图11  频率空间域叠前偏移记录

Fig.11  Zero-migration record with boundary absorbing

图12  频率波数域叠前偏移记录

Fig.12  Zero-migration record without boundary absorbing

使用规则网格有限差分方法计算得到不同时刻垂直位移分量波场快照如图13和14所示。从图13和图14可以看出：在结构损伤处有明显的波场阴影。波场在遇到结构损伤后出现了衍射现象，产生了花纹状波场，但波场整体形状没有剧烈的变化。通过这些衍射现象和阴影区基本可以判断结构损伤的存在和 位置^[14]。

对于炮集叠前深度成像，目前通常采用Marmousi模型数据由有限差分法波动方程数值模拟产生的二维叠前数据体。该数据体包含250炮，每炮90道，每道750个采样点，道间距为20 m，时间采样为4 ms，深度采样为4 m。对Marmousi模型分别采用有限差分方法和混合法进行双平方根算子波动方程的偏移成像计算，如图15和16所示。从成像图中可以看到椭圆型结构夹杂中的直达波、反射波和透射波，还可以观察到由于椭圆型夹杂分界面处近似模拟边界引起的复杂散射波花纹。不规则网格模拟不仅能较好地模拟复杂几何结构的弹性波场传播和分布特征，而且对于一定精度要求的复杂界面正演问题，该方法需要的网格数更少，计算速度更快^[15]。

图13  t=153 μs时水平方向应力波场快照

Fig.13  t=153 μs wave field snapshot of horizon stress

图14  t=176 μs水平方向应力波场快照

Fig.14  t=176 μs wave field snapshot of horizon stress

图15  Marmousi模型双平方根算子频率空间域有限差分法叠前深度成像

Fig.15  Marmousi model double square root operator F-W space field finite difference prestack depth image

图16  Marmousi模型双平方根算子频率空间域混合法叠前深度成像

Fig.16  Marmousi model double square root operator F-W space mix method prestack depth image

3  结论

(1) 采用双平方根算子波动方程可以得到较清晰的叠前偏移成像效果；缺省少数炮集或道对结构损伤成像结果影响很小，采用光滑速度与精确速度模型成像结果基本一致；采用傅里叶有限差分法也可以提高计算速度和叠前偏移成像图成像效果。

(2) 有限差分法计算效率较高，但网格步长必须明显小于最短波长的一半，否则会产生严重的数值频散现象，因此，目前趋向于使用高阶交错网格法。

(3) 波场数值计算速度因波传播方向不同大于或小于实际速度，损伤的存在导致数值频散现象加剧。传播算子是波场在时间上的外推，为提高计算效率，常用深度外推法代替时间外推法。

(4) 波场在通过结构损伤时在叠前偏移距记录上会产生阴影，在波场快照图中产生衍射状的花纹；弹性波在通过结构内部损伤时能量会损耗，表现在波场振幅出现一定程度衰减；结构内部损伤对于动应力集中系数分布规律有明显影响。

参考文献：

[1] 曹晖, 张新亮, 曹永红. 采用改进单元刚度折减系数法识别结构损伤[J]. 振动与冲击, 2008, 27(6): 132-139.

CAO Hui, ZHANG Xinliang, CAO Yonghong. Structural damage identification based on the improved element rigid discount coefficient[J]. Journal of Vibration and Shock, 2008, 27(6): 132-139.

[2] 邹万杰, 翟伟廉. 基于残余力向量法德桁架结构损伤直接识别法[J]. 广西工学院学报, 2009, 20(1): 31-34.

ZOU Wanjie, ZHAI Weilian, A direct damaged identification method of truss structures based on the residue force vector[J]. Journal of Guangxi University of Technology, 2009, 20(1): 31-34.

[3] 郭惠勇, 李正良. 基于测量位移和频率的结构损伤二次识别方法[J]. 振动与冲击, 2007, 26(4): 94-96.

GUO Huiyong, LI Zhengliang. The second identification method of the structural damage based on the tested displacement and frequency[J]. Journal of Vibration and Shock, 2007, 26(4): 94-96.

[4] 王根会, 胡良红. 基于单元模态应变能法的桥梁结构损伤识别[J]. 铁道学报, 2006, 28(3): 83-86.

WANG Genhui, HU Lianghong. Study on damage identification of bridge structures based on the method of element modal strain energy[J]. Journal of the China Railway Society, 2006, 28(3): 83-86.

[5] 李学平, 余志武. 基于动力特性的结构损伤识别方法[J]. 动力学与控制学报, 2006, 4(1): 84-87.

LI Xueping, YU Zhiwu. Structural damage identification method based on dynamic property[J]. Journal of Dynamics and Control, 2006, 4(1): 84-87.

[6] 赵启林, 翟可为. 基于静态信息结构损伤定位的模式匹配法[J]. 计算力学学报, 2006, 23(6): 789-793.

ZHAO Qilin, ZHAI Kewei. Structural damage location model matching method based on static information[J]. Chinese Journal of Computational Mechanics, 2006, 23(6): 789-793.

[7] 冉志红, 缪升. 基于模态柔度法的连续梁结构损伤识别[J]. 公路工程, 2009, 34(2): 34-40.

RAN Zhihong, MIU Sheng. Research of damage identification of continuous beam based on mode flexibility matrix method[J]. Highway Engineering, 2009, 34(2): 34-40.

[8] 焦莉, 李宏男. 基于自适应共振理论的结构损伤识别[J]. 沈阳建筑大学学报, 2005, 21(5): 486-490.

JIAO Li, LI Hongnan. Damage identification of structures based on adaptive resonance theory[J]. Journal of Shenyang Architectures University, 2005, 21(5): 486-490.

[9] Lee U, Shin J. A frequency response function based structural damage identification method[J]. Computers and Structures, 2002: 117-132.

[10] Hudson J. A higher order approximation to the wave propagation constants for a cracked solid[J]. Geophysics, 1986, 87(1): 265-274.

[11] 牛滨华, 何樵登, 孙春岩. 裂隙各向异性介质波动多分量记录的数值模拟[J]. 地球物理学报, 1995, 38(4): 519-527.

NIU Binhua, HE Qiaodeng, SUN Chunyan. Wave multi heft record numerical modeling in Anisotropy media with crack[J]. Journal of Geophysics, 1995, 38(4): 519-527.

[12] M, Douna J. Elastic wave propagation in media with parallel fractures and aligned cracks[J]. Geophysical Prospecting, 1988, 36(6): 571-589.

[13] Lin Xin, Yuan F G. Experimental study applying a migration technique in structural health monitoring[J]. Structural Health Monitoring, 2005, 4(4): 341-353.

[14] Arkadiusz Z, Marek. K, Wieslaw O. Propagation of in-plane in an isotropic panel with a crack[J]. Finite Elements in Analysis and Design, 2006, 42: 929-941.

[15] 张文生. 波动方程成像方法及其计算[M]. 北京: 科学出版社, 2009: 38-45.

ZHANG Wensheng. The image and calculation of the wave equation[M]. Beijing: Science Press, 2009: 38-45.

[16] 孙卫涛. 弹性波动方程的有限差分数值方法[M]. 北京: 清华大学出版社, 2009: 106-118.

SUN Weitao. Finite difference numerical method of elastic wave equation[M]. The Tsinghua University Press, 2009: 106-118.

(编辑  何运斌)

收稿日期：2012-04-20；修回日期：2012-06-25

基金项目：国家自然科学基金资助项目(50778179)

通信作者：郭少华(1960-)，男，陕西西安人，教授，博士生导师，从事结构损伤探测研究；电话：13974814025；E-mail: gsh@mail.csu.edu.cn