CN112083067A - 一种曲面纤维增强树脂基复合材料多向板r区声线示踪算法 - Google Patents

Abstract

一种曲面纤维增强树脂基复合材料多向板R区声线示踪算法，属于复合材料超声检测技术领域。该算法包括将R区沿周向和厚度方向网格化；建立FRP复合材料多向板R区声速关于空间位置和声传播方向角变化的函数关系式，即声速模型；根据起始和目标点之间的相对位置关系，定义网格节点之间的连接关系和连接权重值；根据Fermat最短声时原理，使用Dijkstra最短路径搜索算法从网格节点中搜索计算超声波在起始和目标点之间的传播路径和时间。该算法可实现FRP复合材料多向板R区任意两点间声传播路径和时间的精确快速反演计算，为定量评价材料对声传播行为的影响规律奠定基础，为相控阵超声成像检测反演大量声线路径提供工具。

Description

一种曲面纤维增强树脂基复合材料多向板R区声线示踪算法

技术领域

本发明涉及一种曲面纤维增强树脂基复合材料(Fiber Reinforced Plastic,FRP)多向板R区声线示踪算法，属于复合材料超声检测技术领域。

背景技术

纤维增强树脂基复合材料(Fiber Reinforced Plastics,CFRP)因其轻质、高强和可设计性强等特点，越来越多地代替金属用于飞机、火箭和风力发电机叶片等的承力构件。这些构件普遍尺寸大且形状复杂，存在大量拐角区域(即R区)，是制造和服役阶段的薄弱环节。为实现FRP构件的质量控制，必须借助无损检测技术及时检出R区内部缺陷，其中以超声检测技术应用最为广泛，同时也最具发展潜力。然而，与各向同性的均质金属材料相比，复杂形状FRP构件R区声学特性的特殊性给其超声检测带来了挑战，具体表现为：声速不再是恒定值，而是随传播方向在约3000～9000m/s的范围内变化。工程中所用FRP构件普遍是多向板，相邻二铺层所含纤维的方向不同，在二铺层之间的界面处声学特性失配，倾斜入射至该处的超声波将会因反、折射改变传播方向和声速。同时，由于纤维方向还沿R区周向连续变化，导致声速分布还随位置变化，超声波在单铺层内部也不再保持直线传播，而是呈曲线路径。上述复杂几何形状、多层结构和弹性各向异性三因素间的耦合作用，导致声传播路径复杂难以预测，R区声传播机理不清。为便于操作，工程中普遍是将被检FRP复合材料构件简化为各向同性材料，据此计算所得声传播路径信息失准，进而导致FRP构件超声成像质量差，缺陷定位和定量误差大，甚至缺陷漏检。因此，研究反演计算FRP多向板R区复杂声传播路径的声线示踪方法，是实现声传播行为定量描述，提升缺陷成像质量和缺陷检出能力的迫切需求。

利用Snell定律迭代计算超声在经过界面时的反、折射行为，是计算声传播路径常用的声线示踪解决思路。对于各向同性材料的R区，可直接将R区表面轮廓信息代入至Snell定律公式求解声传播路径。对于CFRP单向板R区，其弹性特性随曲面形状连续变化，内部不存在由声学特性显著差异的区域构成的界面。针对这一特点，Journiac等提出两种声线示踪方案：第一种是将R区沿周向离散化为一系列局域纤维方向恒定的扇形微元，单微元可视作弹性各向异性均质材料，声线在单微元内以直线路径传播，并在经过微元间虚拟界面时发生折射，所以仍可通过使用Snell定律迭代运算计算声传播路径。第二种方法是在借鉴地震学射线理论，将Eikonal和Christoffel方程相结合，针对材料弹性特性沿R区周向连续变化的特点，提出了声线路径及传播方向关于时间的微分公式，实现了初始传播方向已知声线路径的连续化描述。然而，以上两种声线示踪方法若用于后处理重建超声图像，需执行大量迭代运算，效率低，应用受限。而且，对于声学特性更加复杂的FRP多向板R区，声线路径除了在单铺层内部连续改变传播方向外，还会在经过铺层界面时发生偏折，声线路径较FRP单向板R区更加复杂。

Fermat原理指出，超声波总是沿着用时最短的路径在空间两定点之间传播。因此，布里斯托大学Nower等提出了利用计算机领域中的Dijkstra最短路径搜索算法来解决超声检测领域声线示踪问题的新思路。该方法无需迭代计算声线在界面处的反、折射行为，并具备最短路径搜索算法执行高精度快速计算的优势，非常适合于解决几何形状和声学特性复杂介质的全局复杂声线路径快速反演求解问题。林莉等已将该方法成功用于CFRP多向平板声线示踪。因此，开展Dijkstra最短路径搜索算法用于FRP复合材料多向板R区声线示踪的研究，具有十分重要的现实意义。

发明内容

本发明提出一种曲面FRP复合材料多向板R区声线示踪算法，可实现任意两点之间复杂声传播路径的精确、快速反演计算，为定量评价材料对声传播行为的影响规律奠定基础，为相控阵超声成像检测反演大量声线路径提供工具。

本发明采用的技术方案是一种曲面纤维增强树脂基复合材料多向板R区声线示踪算法，其特征在于：将R区沿周向和厚度方向网格化；建立FRP复合材料多向板R区声速关于空间位置和声传播方向角变化的函数关系式(即声速模型)；根据起始和目标点之间的相对位置关系，定义网格节点之间的连接关系和连接权重值；根据Fermat最短声时原理，使用Dijkstra最短路径搜索算法从网格节点中搜索出超声波在起始和目标点之间的传播路径和时间。当声源位于R区凸面侧表面上方区域时，超声波由声源传播至R区不同位置目标点的路径，可通过执行如下步骤计算得到：

a.FRP复合材料多向板R区网格化

以FRP铺层界面为周向网格化边界，厚度方向间隔为△d，周向边界序号记为I，其中I为正整数，从上往下记录；

在此基础上，厚度方向边界以Δα的周向角度间隔将R区网格化为扇形微元，以Δl的周向长度间隔将R区左右两侧平板划分为矩形单元，厚度方向边界序号记为J，其中J为正整数，从左往右记录；

当Δα和Δl足够小时，单个扇形微元或矩形单元内部纤维方向恒定，可等效为各向异性的均质材料，其内部声传播路径为直线，同时声速随方向变化；Δα和Δl越小，声线示踪计算精度越高；

b.建立声速关于空间位置和声传播方向角变化的函数关系式，即声速模型

根据FRP复合材料多向板的几何尺寸、密度、铺层数量、铺放顺序、单铺层厚度和弹性刚度矩阵，计算声速关于空间位置和声传播方向角θ的变化情况，即V(x,y,θ)，其中x和y分别表示网格节点空间位置的x和y轴坐标；

c.判断声线路径类型，并指定网格节点之间的连接关系

记起始点和目标点分别为s和e，记s点所处周向和厚度方向的边界序号分别为I_s和J_s，记e点所处周向和厚度方向的边界序号分别为I_e和J_e。根据s和e之间的相对位置关系，声线路径可分为4类，具体包括：

(1)第1类，满足条件I_s＝0&I_e＝1&(x_s-x_e)x_e+(y_s-y_e)y_e≥0，此时目标点e位于R区凸面侧表面，同时起始点s发射的声线可直接经由探头与R区之间的声耦合层直线传播至e，声线路径为连接s和e的直线；

(2)第2类，满足条件I_s＝0&I_e＝1&(x_s-x_e)x_e+(y_s-y_e)y_e<0，此时目标点e同样是位于R区凸面侧表面，但起始点s发射的声线会因R区凸面侧表面的遮挡，无法直线传播至e，而是先大角度倾斜入射至表面某点，折射后近乎沿着表面轮廓继续向前传播，直至到达e。声线路径可分为两部分：第一部分位于耦合层内部，为连接s和表面入射点的直线；第二部分是CFRP多向板表面轮廓位于表面入射点和e之间的分段。因此，声线路径由满足如下条件的可直接相连网格节点组合p和q构成：

①若I_p＝0&I_q＝1&(x_p-x_q)cos(α_q)+(y_p-y_q)sin(α_q)≥0，即p位于起始点，q位于FRP复合材料多向板表面，其中α_q表示经过q的表面轮廓法线的周向角度；

②若x_ey_s-x_sy_e>0，即e位于s右侧，则声线路径必定位于表面轮廓分段I＝1&J≤J_e内，网格节点p和q可直接相连的条件是I_q＝I_p＝1&J_p＝J_q-1；

③若x_ey_s-x_sy_e<0，即e位于s左侧，则声线路径必定位于表面轮廓分段I＝1&J≥J_e内，网格节点p和q可直接相连的条件是I_p＝I_q＝1&J_p＝J_q+1；

(3)第3类，满足条件I_s＝0&I_e≥1&(x_s-x_e)y_e-(y_s-y_e)x_e＝0，此时目标点e位于R区内部，并与起始点s位于同一径向。由于不同纤维方向FRP铺层沿厚度方向的声速相同，沿此方向传播的声线不改变方向。所以，s发射的声线可沿直线传播至e，声线路径为连接s和e的直线。声线路径可分为两部分：第一部分位于耦合层内部，为连接s和表面入射点的直线；第二部分位于FRP多向板内部，为连接表面入射点和e的直线；

(4)第4类，满足条件I_s＝0&I_e>1&(x_s-x_e)y_e-(y_s-y_e)x_e≠0，此时目标点e位于R区内部，而且不与起始点s位于同一径向，声线将偏离厚度方向在两点之间传播，在弯曲单铺层内部连续改变传播方向和声速，在声学特性失配铺层间界面处也会因折射改变传播方向和声速。所以，声线路径由满足如下条件的可直接相连网格节点组合p和q构成：

①若I_p＝0&I_q＝1&(x_p-x_q)cos(α_q)+(y_p-y_q)sin(α_q)≥0，即p位于起始点，q位于FRP复合材料多向板表面；

②若x_ey_s-x_sy_e>0，即e位于s右侧，则声线路径必定位于区域1≤I≤I_e&J≤J_e内，网格节点p和q可直接相连的条件是I_p≥1&I_q＝I_p+1&J_p<J_q；

③若x_ey_s-x_sy_e<0，即e位于s左侧，则声线路径必定位于区域1≤I≤I_e&J≥J_e内，网格节点p和q可直接相连的条件是I_p≥1&I_q＝I_p+1&J_p>J_q。

d.计算直接相连网格节点p和q之间的声线路径和声时

具体包括以下6种情况：

(1)对于满足第1类声线路径的情况，网格节点p和q分别位于s和e点，相应的声时为两点之间的直线距离与耦合层声速的比值，因此声时t_pq(即t_se)可直接由如下公式计算得到：

(2)对于满足第3类声线路径的情况，网格节点p和q分别位于s和e点，声时t_pq(即t_se)可直接由如下公式计算得到：

其中，R表示R区凸面侧表面曲率半径；

(3)对于第2和4类声线路径满足条件①的情况，p和q之间的声线路径和声时计算方法与第1类声线路径情况相同；

(4)对于第2类声线路径满足条件②或③的情况，p和q位于同一均质的网格单元，所以声线路径为连接二者的直线，声时t_pq可直接由公式(1)计算得到；

(5)对于第4类声线路径满足条件②或③的情况，p和q所处单铺层可视为“FRP复合材料单向板R区”，将p和q所处单铺层进行离散化，并将单铺层等效为由纤维方向恒定的矩形和扇形微元沿周向铺叠而成的多层结构，单微元可视为弹性各向异性均匀介质；指定只位于同一微元的新节点之间可直接连接，并以超声波在两点间直线传播所需时间为连接权重值；利用Dijkstra算法从所有直接连接新节点组合中搜索出p和q之间具有全局最短声时的声线路径，相应声时为构成声线路径的新节点组合连接权重值之和；

(6)指定网格节点之间的连接关系是单向的，即若p可直接与q相连，则q就不能与p直接相连，而且t_qp＝∞；

e.利用Dijkstra算法搜索起始和目标点之间的声线路径

在获知可直接相连网格节点p和q之间的声线路径和声时的基础上，以Fermat原理，两点之间实际声线路径具有全局最短声时为准则，利用Dijkstra最短路径搜索算法计算得到第2和4类s至e声线路径和声时。具体步骤如下：

(1)计算区域内的所有网格节点根据其与起始点之间具有全局最短声时的声线路径是已知还是未知，分别放置于集合S和E；然后，指定超声波由起始点s传播至计算区域内任意节点k所需声时为t_sk＝∞且t_ss＝0，即起初只有s与其自身之间的全局最短声时路径是已知的；

(2)将所有节点分配至集合E，这时集合S是空集；

(3)找出当前集合E中具有最短声时的节点k，对于集合E中的其余任意节点l，根据k和l之间的连接关系，重建出s至l的新路径，相应的声时等于t_sk+t_kl；若s至l当前路径的声时大于新路径的声时，即t_sl>t_sk+t_kl，说明新路径比当前路径更加接近于具有全局最短声时的实际路径，故以新路径替换当前路径，并指定t_sl＝t_sk+t_kl；

(4)对于当前集合E中具有最短声时的节点k，已不存在经由其它节点且声时更短的新路径，所以当前路径即为s至k的全局最短声时路径，将k由集合E转移至集合S；

(5)重复步骤(3)和(4)，直至集合E中的所有节点都被转移至集合S，得到起始和目标点之间具有全局最短声时的实际声线路径。

“FRP复合材料单向板R区”中p和q之间的声线路径和声时可使用“一种曲面纤维增强树脂基复合材料单向板R区声线示踪算法”中提出的方法进行处理。

本发明的有益效果是：该算法将R区沿周向和厚度方向网格化；建立曲面FRP复合材料多向板R区声速关于空间位置和声传播方向角变化的函数关系式(即声速模型)；根据起始和目标点之间的相对位置关系，定义网格节点之间的连接关系和连接权重值；根据Fermat最短声时原理，使用Dijkstra最短路径搜索算法从网格节点中搜索出超声波在起始和目标点之间的传播路径和时间。该算法可实现FRP多向板R区任意两点之间复杂声传播路径的精确、快速反演计算，为定量评价材料对声传播行为的影响规律奠定基础，为相控阵超声成像检测反演大量声线路径提供工具。

附图说明

下面结合附图，以线性阵列相控阵超声探头检测CFRP多向板R区为实施例对本发明作进一步说明。

图1是CFRP多向板R区相控阵超声检测实验条件示意图。

图2是CFRP多向板R区网格化原理示意图。

图3是CFRP多向板R区典型声线路径示意图。

图4是CFRP多向板R区(0,17.08)点声速分布曲线。

图5是CFRP多向板R区声线路径计算结果。

具体实施方式

如图1所示，实验使用5L64-NW1型线性阵列探头(5MHz中心频率，共64阵元，单阵元尺寸为0.92mm×7mm，相邻阵元边缘距0.08mm，相邻阵元中心距p＝1mm)1～31号阵元作为活动阵元，探头水平放置于CFRP多向板试样R区凸面侧上方(即探头与R区右侧平板表面之间的夹角δ＝π/4)，与R区凸面侧表面之间的竖直距离D＝10mm，二者之间以水填充(声速1483.40m/s)。CFRP多向板试样R区凸面侧表面曲率半径为R＝8.80mm，左右两侧平板长度分别为10mm。定义x轴水平向右，y轴竖直向上，坐标系原点o位于R区曲率圆心，R区左右两侧厚度方向边界对应的周向角度分别为α_l＝3π/4和α_r＝π/4。铺放顺序为[-45/45/0/45/90/0/90/-45/90/90/-45/45/0/45/90/0/90/-45/90/90/45/90/-45/90]s，共48层，单铺层厚度为0.18mm，密度为1.592g/cm³，0°铺层的弹性刚度矩阵C如下：

(1)FRP复合材料多向板R区网格化

如图2所示，网格化以“探头/水”界面和FRP铺层界面为周向网格化边界，记周向边界序号记为I，其中I为正整数(1≤I≤50)，从上往下记录。在此基础上，厚度方向边界以Δα＝0.1°的周向角度间隔将R区网格化为扇形微元，以Δl＝0.1mm的周向长度间隔将R区左右两侧平板划分为矩形单元，厚度方向边界序号记为J，其中J为正整数(1≤J≤1101)，从左往右记录；

(2)建立声速关于空间位置和声传播方向角变化的函数关系式(声速模型)

根据试样的几何尺寸、密度、铺层数量、铺放顺序、单铺层厚度和弹性刚度矩阵，通过Christoffel方程计算声速关于空间位置和声传播方向角θ的变化情况，即V(x,y,θ)。以位于(0,17.08)的网格节点为例，局域声速空间分布曲线如图3所示；

(3)根据起始和目标点之间的相对位置关系，判断声线路径类型，明确声线路径可能经过的区域，并定义网格节点之间的连接关系

选择线性阵列探头1～31号阵元为起始点，坐标为(n-16,23.89)，其中n＝1,2,...,31，以位于R区内部的A(-5.533,16.53)、B(-3.712,17.04)和M(0,13.12)点为目标点进行声线示踪。起始和目标点之间的声传播路径及所需时间可通过执行如下步骤获得：

1)当选择以第7号阵元(-9,23.89)为起始点(I＝1,J＝7)，以位于R区凸面侧表面的点A(-5.533,16.53)为目标点(I＝2,J＝716)时，为第1类声线路径(图4)，此时超声波完全在水中传播，相应的声线路径为连接二者的直线，声时可直接由公式(1)计算得到：

表1 7号阵元与A点之间的声线路径和声时

起始点	目标点	水中传播距离/mm	水声速度/m·s<sup>-1</sup>	路径	声时/μs
						7号阵元	A	8.1357	1483.4	I＝1,2&J＝7,716	5.48

2)当选择以第31号阵元(15,23.89)为起始点(I＝1,J＝31)，以位于R区凸面侧表面的点B(-3.712,17.04)为目标点(I＝2,J＝428)时，为第2类声线路径(图4)，可分为两部分：第一部分，声线经过水耦合层以临界角度倾斜入射至R区凸面侧表面点(8.34,15.32)，声线路径为连接s和表面入射点的直线，声时可直接由公式(1)计算得到；第二部分，折射声线近乎沿着R区凸面侧表面轮廓继续向前传播，直至到达e，声线路径为R区凸面侧表面轮廓位于入射点和e点之间的分段，声时等于弧形路径长度与CFRP多向板R区沿此方向声速的比值：

表2 30号阵元与B点之间的声线路径和声时

3)当选择以第16号阵元(0,23.89)为起始点(I＝1,J＝7)，以位于R区凸面侧表面的点M(0,13.12)为目标点(I＝26,J＝551)时，由于其与目标点A位于R区同一径向，为第3类声线路径(图4)，相应的声线路径为连接二者的直线，声时可直接由公式(2)计算得到：

表3 16号阵元与M点之间的声线路径和声时

4)当选择以第1号阵元(-15,23.89)为起始点(I＝1,J＝1)，以位于R区凸面侧表面的点M(0,13.12)为目标点(I＝26,J＝551)时，为第4类声线路径(图4)，需使用Dijkstra算法搜索声线路径，具体执行如下步骤：

①由于x_ey_s-x_sy_e>0，即e位于s右侧，则声线路径必定位于区域1≤I≤I_e&J≤J_e内，网格节点p和q可直接相连的条件是I_p≥1&I_q＝I_p+1&J_p<J_q，或I_p＝0&I_q＝1&(x_p-x_q)cos(α_q)+(y_p-y_q)sin(α_q)≥0，即p位于起始点，q位于FRP复合材料多向板表面；

②计算区域内的所有网格节点根据其与起始点之间具有全局最短声时的声线路径是已知还是未知，分别放置于集合S和E。然后，指定超声波由s点(即第1号阵元)传播至计算区域内任意节点k所需声时为t_sk＝∞且t_ss＝0，即起初只有s点与其自身之间的全局最短声时路径是已知的；

③将所有节点分配至集合E，这时集合S是空集；

④找出当前集合E中具有最短声时的节点k，那么对于集合E中的其余任意节点l，这时便有可能根据k和l之间的连接关系，重建出s至l的新路径，相应的声时等于t_sk+t_kl。若s至l当前路径的声时大于新路径的声时，即t_sl>t_sk+t_kl，说明新路径比当前路径更加接近于具有全局最短声时的实际路径，故以新路径替换当前路径，并指定t_sl＝t_sk+t_kl；

⑤对于当前集合E中具有最短声时的节点k，已不可能存在经由其它节点且声时更短的新路径，所以当前路径即为s至k的全局最短声时路径，并需将k由集合E转移至集合S。

起初具有最短声时的节点k即为s点，与其直接相连的节点l位于R区凸面侧表面，且满足条件I_l＝2&J_l≤J_e＝551。以其中分别满足I＝2&J＝421和I＝2&J＝551条件的节点为例，其声线路径、声时和当前路径是否为最短声时路径的情况如表4所示：

表4利用Dijkstra算法搜索1号阵元与其邻近直接相连节点之间的声线路径

⑥重复步骤④和⑤，直至集合Q中的所有节点都被转移至集合P，这时便可得到起始和目标点之间具有全局最短声时的实际声线路径，结果如图5所示，相应的声时为9.43μs。

Claims (1)

1.一种曲面纤维增强树脂基复合材料多向板R区声线示踪算法，其特征在于：将R区沿周向和厚度方向网格化；建立FRP复合材料多向板R区声速关于空间位置和声传播方向角变化的函数关系式，即声速模型；根据起始和目标点之间的相对位置关系，定义网格节点之间的连接关系和连接权重值；根据Fermat最短声时原理，使用Dijkstra最短路径搜索算法从网格节点中搜索计算超声波在起始和目标点之间的传播路径和时间；具体执行如下步骤：

a.FRP复合材料多向板R区网格化

以FRP铺层界面为周向网格化边界，厚度方向间隔为△d，周向边界序号记为I，其中I为正整数，从上往下记录；

在此基础上，厚度方向边界以Δα的周向角度间隔将R区网格化为扇形微元，以Δl的周向长度间隔将R区左右两侧平板划分为矩形单元，厚度方向边界序号记为J，其中J为正整数，从左往右记录；

当Δα和Δl足够小时，单个扇形微元或矩形单元内部纤维方向恒定，等效为各向异性的均质材料，其内部声传播路径为直线，同时声速随方向变化；Δα和Δl越小，声线示踪计算精度越高；

b.建立声速关于空间位置和声传播方向角变化的函数关系式，即声速模型

根据FRP复合材料多向板的几何尺寸、密度、铺层数量、铺放顺序、单铺层厚度和弹性刚度矩阵，计算声速关于空间位置和声传播方向角θ的变化情况，即V(x,y,θ)，其中x和y分别表示网格节点空间位置的x和y轴坐标；

c.判断声线路径类型，并指定网格节点之间的连接关系

记起始点和目标点分别为s和e，记s点所处周向和厚度方向的边界序号分别为I_s和J_s，记e点所处周向和厚度方向的边界序号分别为I_e和J_e；根据s和e之间的相对位置关系，声线路径可分为4类，具体包括：

(1)第1类，满足条件I_s＝0&I_e＝1&(x_s-x_e)x_e+(y_s-y_e)y_e≥0，此时目标点e位于R区凸面侧表面，同时起始点s发射的声线可直接经由探头与R区之间的声耦合层直线传播至e，声线路径为连接s和e的直线；

(2)第2类，满足条件I_s＝0&I_e＝1&(x_s-x_e)x_e+(y_s-y_e)y_e<0，此时目标点e同样是位于R区凸面侧表面，但起始点s发射的声线会因R区凸面侧表面的遮挡，无法直线传播至e，而是先大角度倾斜入射至表面某点，折射后近乎沿着表面轮廓继续向前传播，直至到达e；声线路径可分为两部分：第一部分位于耦合层内部，为连接s和表面入射点的直线；第二部分是CFRP多向板表面轮廓位于表面入射点和e之间的分段；因此，声线路径由满足如下条件的可直接相连网格节点组合p和q构成：

①若I_p＝0&I_q＝1&(x_p-x_q)cos(α_q)+(y_p-y_q)sin(α_q)≥0，即p位于起始点，q位于FRP复合材料多向板表面，其中α_q表示经过q的表面轮廓法线的周向角度；

②若x_ey_s-x_sy_e>0，即e位于s右侧，则声线路径必定位于表面轮廓分段I＝1&J≤J_e内，网格节点p和q可直接相连的条件是I_q＝I_p＝1&J_p＝J_q-1；

③若x_ey_s-x_sy_e<0，即e位于s左侧，则声线路径必定位于表面轮廓分段I＝1&J≥J_e内，网格节点p和q可直接相连的条件是I_p＝I_q＝1&J_p＝J_q+1；

(3)第3类，满足条件I_s＝0&I_e≥1&(x_s-x_e)y_e-(y_s-y_e)x_e＝0，此时目标点e位于R区内部，并与起始点s位于同一径向；由于不同纤维方向FRP铺层沿厚度方向的声速相同，沿此方向传播的声线不改变方向；所以，s发射的声线可沿直线传播至e，声线路径为连接s和e的直线；声线路径可分为两部分：第一部分位于耦合层内部，为连接s和表面入射点的直线；第二部分位于FRP多向板内部，为连接表面入射点和e的直线；

(4)第4类，满足条件I_s＝0&I_e>1&(x_s-x_e)y_e-(y_s-y_e)x_e≠0，此时目标点e位于R区内部，而且不与起始点s位于同一径向，声线将偏离厚度方向在两点之间传播，在弯曲单铺层内部连续改变传播方向和声速，在声学特性失配铺层间界面处也会因折射改变传播方向和声速；所以声线路径由满足如下条件的可直接相连网格节点组合p和q构成：

①若I_p＝0&I_q＝1&(x_p-x_q)cos(α_q)+(y_p-y_q)sin(α_q)≥0，即p位于起始点，q位于FRP复合材料多向板表面；

②若x_ey_s-x_sy_e>0，即e位于s右侧，则声线路径必定位于区域1≤I≤I_e&J≤J_e内，网格节点p和q可直接相连的条件是I_p≥1&I_q＝I_p+1&J_p<J_q；

③若x_ey_s-x_sy_e<0，即e位于s左侧，则声线路径必定位于区域1≤I≤I_e&J≥J_e内，网格节点p和q可直接相连的条件是I_p≥1&I_q＝I_p+1&J_p>J_q；

d.计算直接相连网格节点p和q之间的声线路径和声时

具体包括以下6种情况：

(1)对于满足第1类声线路径的情况，网格节点p和q分别位于s和e点，相应的声时为两点之间的直线距离与耦合层声速的比值，因此声时t_pq(即t_se)可直接由如下公式计算得到：

(2)对于满足第3类声线路径的情况，网格节点p和q分别位于s和e点，声时t_pq(即t_se)可直接由如下公式计算得到：

其中，R表示R区凸面侧表面曲率半径；

(3)对于第2和4类声线路径满足条件①的情况，p和q之间的声线路径和声时计算方法与第1类声线路径情况相同；

(4)对于第2类声线路径满足条件②或③的情况，p和q位于同一均质的网格单元，所以声线路径为连接二者的直线，声时t_pq可直接由公式(1)计算得到；

(5)对于第4类声线路径满足条件②或③的情况，p和q所处单铺层视为“FRP复合材料单向板R区”，首先将p和q所处单铺层进行离散化，并将单铺层等效为由纤维方向恒定的矩形和扇形微元沿周向铺叠而成的多层结构，单微元可视为弹性各向异性均匀介质；指定只位于同一微元的新节点之间直接连接，并以超声波在两点间直线传播所需时间为连接权重值；利用Dijkstra算法从所有直接连接新节点组合中搜索出p和q之间具有全局最短声时的声线路径，相应声时为构成声线路径的新节点组合连接权重值之和；

(6)指定网格节点之间的连接关系是单向的，即若p可直接与q相连，则q就不能与p直接相连，而且t_qp＝∞；

e.利用Dijkstra算法搜索起始和目标点之间的声线路径

在获知可直接相连网格节点p和q之间的声线路径和声时的基础上，以Fermat原理，两点之间实际声线路径具有全局最短声时为准则，利用Dijkstra最短路径搜索算法计算得到第2和4类s至e声线路径和声时，具体步骤如下：

(1)计算区域内的所有网格节点根据其与起始点之间具有全局最短声时的声线路径是已知还是未知，分别放置于集合S和E；然后，指定超声波由起始点s传播至计算区域内任意节点k所需声时为t_sk＝∞且t_ss＝0，即起初只有s与其自身之间的全局最短声时路径是已知的；

(2)将所有节点分配至集合E，这时集合S是空集；

(3)找出当前集合E中具有最短声时的节点k，对于集合E中的其余任意节点l，根据k和l之间的连接关系，重建出s至l的新路径，相应的声时等于t_sk+t_kl；若s至l当前路径的声时大于新路径的声时，即t_sl>t_sk+t_kl，说明新路径比当前路径更加接近于具有全局最短声时的实际路径，故以新路径替换当前路径，并指定t_sl＝t_sk+t_kl；

(4)对于当前集合E中具有最短声时的节点k，已不存在经由其它节点且声时更短的新路径，所以当前路径即为s至k的全局最短声时路径，将k由集合E转移至集合S；

(5)重复步骤(3)和(4)，直至集合E中的所有节点都被转移至集合S，得到起始和目标点之间具有全局最短声时的实际声线路径。

Images