飞机结冰是指在飞行或停留地面过程中，在特定环境条件下发生的表面结冰现象。一些研究表明，机翼表面的结冰会影响飞机的气动性能、稳定性和操纵性[1]，在进气口等关键区域结冰会导致严重的飞行事故[2]。因此，及早发现并清除结冰对提高飞行效率和安全具有重要意义。然而，现有的除冰技术通常不能长时间工作，如电热法除冰能耗很高，机械法除冰会改变飞机的气动外形，液体除冰会产生一定的环境污染。因此，除冰技术必须依赖于准确地积冰探测，从而选定合适的除冰时间，提高除冰效率，确保飞行安全。到目前为止，基于不同的技术原理，已经开发出了几种在飞机蒙皮上进行结冰检测的方法，包括基于飞行过程中气动参数变化的参数辨识算法，如批量最小二乘辨识算法、卡尔曼滤波(KF)算法、H∞辨识算法、神经网络(NN)算法、神经网络/KF组合算法[3-6]等。此外，还有各种新型传感器的运用给积冰探测技术提供了新的可能，Zou等[7]报道了一种新型的光纤冰传感器，其输出电压会随着玻璃冰厚度的增加而迅速降低。同样，Prasad等[8]提出了使用光纤传感技术检测和监测冰层形成。此外，还有一些独特的积冰探测方法，Dongo等[9]提出了热脉冲冰探测系统的概念，该系统通过记录温度变化并评估表面受到电脉冲产生的热激励后的时间延迟来探测积冰形成。Liu等[10]应用超声脉冲回波法测量了随频率变化的衰减系数，证明了探测冰类型的可行性。Schlegl等[11]提出了一种用于飞机表面的太阳能无线冰层探测和温度测量系统，该系统基于空气、水和冰之间电特性的差异。超声导波是超声波在波导结构中传播的一种特殊形式，主要特点包括传播距离广、能耗小、方向性强等，因此常用于大型结构的损伤检测。过去几十年，基于超声导波的结构损伤检测技术发展很快，并逐渐走向成熟和具体的工业应用。近年来，随着超声导波损伤检测技术的发展和成熟，超声导波在冰层检测中的潜在应用也引起了人们的极大关注。Zhao等[12]使用皮尔逊相关系数来比较结冰前后的导波信号。Munoz等[13]尝试了各种算法来提取导波信号特征，如能量比、自相关等，以识别冰的形状。Jimenez等[14]开发了一种基于导波模式识别检测和分类冰层厚度的方法。类似地，在信号特征分析中，Siavash等[15]提出了一种对冰层敏感的基于信号积分差异的指标参数。虽然这些研究在时域分析或试验中取得了较好的冰检测效果，但这些基于数据处理的方法不能从理论上解释信号变化的原因。本文建立了频域有限元模型，通过仿真定量分析了结冰铝板中单一模式超声导波的模式转换过程，并从声传播机理的角度解释了结冰引起信号变化的原因。然后通过时域分析验证了导波冰层对信号到达时间的影响，验证了模式转换过程分析的正确性。与基于信号能量变化的方法相比，该方法避免了耦合剂、探头连接性等外部因素的影响，有效提高了积冰检测的稳定性和精确性。本研究的主要贡献包括两部分：基于频域有限元定量分析了UGW在结冰铝板中的模式转换；基于群速度变化研究了积冰层的识别和评估。1　理论模型本文旨在研究超声导波探测机翼蒙皮表面的积冰问题，分析超声导波传播经过积冰区域后的变化情况，根据采集到的数据分析积冰参数。铝和冰层的材料属性见表1，对于各向同性材料，Lamb和SH波解耦，其中Lamb仅包含ux和uz，SH波仅包含uy，因此讨论Lamb问题时，可以把模型简化为二维，忽略y方向的位移，如图1所示。10.12050/are20220316.T001表1材料属性[16]Table 1Material properties16]材料弹性模量E/GPa泊松比υ密度ρ/(kg/m3)铝700.332700冰8.30.35190010.12050/are20220316.F001图1三维模型简化为二维模型Fig.13D model is simplified to 2D model超声波在特定的波导结构（如板、管等）中传播时会形成导波，与体波相比，导波能量衰减小、传播距离更远，适合大型结构的监测。如图1所示为冰-铝双层结构的二维简化模型，其中铝板厚度为2mm，冰层厚度为0.5mm。可以用基于部分波法的全局矩阵法计算该结构的频散曲线，并与纯铝板的频散曲线进行对比，如图2所示。10.12050/are20220316.F002图2群速度频散曲线Fig.2Dispersion curves for group velocity如图2中的群速度曲线表示对应导波信号点的能量传播速度，也是波包的传播速度。在实际探测中，常用的探测指标有幅值衰减、信号到达时间、模式截止频率等[16]。其中幅值衰减常用于表面附着物的探测，但是，导致信号幅值变化的因素很多，包括探头接触不良、信号干涉、环境温度等，探测效果容易受到这些因素的干扰。而模式截止频率所能够提供的有效信息很少，不能用来进一步判断积冰层的尺寸参数。与之相比，信号到达时间具有更高的鲁棒性[17]，同时与积冰层的长度具有相关性，可以用来准确判断积冰层的出现和生长。图中区域A、B为选定的最佳区域，其主要特征包括：（1）激励信号模式的群速度最大，并且与同频率下其他模式的群速度有明显差异。保证接收信号的第一个波包为有效信号，避免不同模式的波包叠加。（2）与积冰区域对应的模式群速度存在差异，确保积冰出现后，信号到达时间发生明显变化。2　频域有限元模型超声导波的频散特性导致在波导结构中激励特定频率的信号时，会同时产生多种模式的信号，不同模式的信号分别具有不同的相速度和群速度，相互叠加使信号变得非常复杂。并且，当导波进入积冰结构中以后，信号模式发生转换，信号复杂度更高，不利于信号处理和分析。在仿真分析中，使用时域分析模型不可避免地会遇到上述问题，部分研究者使用相控阵激励的方式，通过相位叠加提高特定模式的激发能量，可以在一定程度上降低其他模式的干扰，但仍然不能完全避免。频域分析能够产生频散曲线中任意一点对应的频率和模式，从而准确分析任意单一模式导波信号的传播过程。在COMSOL®中建立了频域有限元模型，如图3所示。该模型包括2mm厚度的铝层和0.5mm厚度的冰层，分为6个部分：加载区域内施加载荷，以指定的频率-模式激励产生单一模式导波，该区域的长度为100mm；入射区和出射区为单层铝板，该区域的长度为800mm；结冰区为包含冰层和铝层的结构，长度为800mm；两端为完美匹配层（perfectly matched layer），可以吸收反射信号，模拟无限长边界，该区域的长度为100mm。结构整体使用二阶单元划分网格，单元尺寸为0.1mm。10.12050/are20220316.F003图3频域有限元模型Fig.3Frequency domain finite element model加载区体载荷的分布与选定模式的波结构相关。根据式(1)计算得到应力分布σ31z,σ33z和位移分布u1z,u3z，根据各向同性材料的本构方程，可以推导所需的应力ε11z=iku1zε33z=∂u3z∂zσ11z=λε11z+ε33z+2με11zσ13z=σ31z （1）式中：ε, σ, k和 u分别为波数、位移、应变和应力，λ 和 μ为拉梅常数（Lame Constants）。加载区的体载荷分布可表示为f1(x,z)=σ11zexp(i2πx/cp)f3x,z=σ13zexp(i2πx/cp) （2）根据选定的频率和模式计算相应的波结构和相速度并代入式(2)，能够产生所需的单一频率和模式的信号。1201kHz-Mode A1激励信号的位移云图如图4所示，从图4中可以看出，在加载区域，x方向的位移大小逐渐增加，代表波的传播方向为x方向；在入射区域形成了稳定且规律的位移云图，表示只存在一种模式，分析云图中的波长可以看出得到波频率与激励频率一致；在积冰区域，云图发生了变化，表示出现了新的模式；出射区同样发生变化，与入射区和积冰区域都有不同，可以看出是多种模式叠加的结果。10.12050/are20220316.F004图4激励信号为1201kHz-mode A1的x方向位移云图Fig.4Displacement contour of x direction for thecase (1201kHz, mode A1)在此，引入坡印亭矢量（Poynting Vector）来量化分析积冰铝板中的模式转换情况。坡印亭矢量可以表示固体中声波的能流密度[18]，通过计算不同区域内的坡印亭矢量可以量化分析结构中各区域内的不同方向上的能流分布，从而得到实际存在的导波模式数量以及不同模式的能量大小。P=- 12ν*σ （3）ν=v1v2v3, σ=σ11σ12σ13σ21σ22σ23σ31σ32σ33 （4）式（4）为坡印亭矢量的表达式，其中v*,σ分别表示速度矢量矩阵的共轭和应力矢量矩阵。波印亭矢量的具体计算步骤如下：（1）从频域仿真结果中获得指定区域内的速度矢量矩阵vx, z和应力矢量矩阵σx, z；（2）将采样数据沿x方向做快速傅里叶变换（FFT）得到波数域的速度矢量矩阵Fv(k, z)和应力矢量矩阵Fσ(k, z)，再沿z方向求和，得到波数域的度矢量矩阵Fv(k)和应力矢量矩阵Fσ(k)；（3）将应力矢量矩阵代入到坡印亭矢量表达式中，得到结果P(k)；（4）最后通过式（4）将结果转化为相速度域P(cp)。cp=ωk （5）式中：cp，ω和k分别为相速度、角频率和波数。3　频域仿真结果分析本文从频散曲线中300~2000kHz范围选取信号点进行扫频分析，其中包括mode A0,S0,A1,S1，分别计算波结构，并输入到频域模型中，通过仿真计算分析对应信号点在积冰铝板结构中的传播特性。其中，入射区和出射区是只有单层铝板的对称结构，导波模式可以划分为对称模式和反对称模式，分别依次记为S0,S1等和A0,A1等；积冰区为非对称结构，导波模式可依次记为mode 1、mode 2等。区域A主要包括mode S0的低频部分，频率范围在1MHz以下。在低频区域，薄冰层的出现会使相速度和群速度降低，但影响很小。从图5(a)~图5(d)中可以看出，4种低频情况（300kHz-mode S0, 400kHz-mode S0, 500kHz-mode S0, 530kHz-mode S0）下的模式转换情况基本一致。以图5(c)为例，其激励信号为500kHz-mode S0，此时导波在铝板和积冰铝板中传播的群速度分别为5118.63m/s和4894.9m/s，积冰层对导波传播速度的影响较小。同时，在三个区域内均只存在一个主要模式。图5(e)中激励信号为854kHz-mode S0情况下，虽然在积冰区存在较小的mode 3，但明显与主要模式mode 2差距很大。因此，仍可视为在三个区域内均稳定存在一个主要模式，即mode S0-mode 2-mode S0。图5激励信号为低频时三个区域内的波印亭矢量Fig.5Poynting vector of three areas for the low frequency cases10.12050/are20220316.F6a110.12050/are20220316.F6a210.12050/are20220316.F6a310.12050/are20220316.F6a410.12050/are20220316.F6a5 与之相比，区域B主要包括高频部分的模式3，其频率范围在1MHz以上。在此范围内，铝板中可能存在三种模式（A0，S0，A1），积冰铝板中可能存在4种模式。因此，该区域的模式转换过程更加复杂。如图6所示分别为该区域中三种不同频率的激励信号的传播模式，其中，激励频率为1101kHz-mode A1时，如图6(a)所示，经过第一次模式转换，信号在积冰区域中产生两种模式的信号，其中模式3为主要模式。经过第二次模式转换，导波在出射区中存在三个模式，其中A1为主要模式，包含大部分能量。同时，入射区和积冰区存在明显的反射信号，说明导波传播过程中可能会存在明显的信号衰减；当激励频率为1178kHz-mode A1时，如图6(b)所示，导波的模式转换过程与上述情况相似，在积冰区存在两个模式，其中模式 4 为主要模式。在出射区存在三个模式，且S0和A1模式包含的能量相近；当激励频率为1239kHz-mode A1时，如图6(c)所示，导波传播中的模式转换过程明显更加平滑，三个区域均只存在一个主要模式，分别为A1、mode 4和A1，且没有明显的反射信号。10.12050/are20220316.F005图6激励信号为高频时三个区域内的波印亭矢量Fig.6Poynting vector of three areas for the high frequency cases综上所述，低频导波信号的模式转换过程比较简单平滑，在三个区域内均只存在一个主要模式，可以据此确定不同区域内的导波群速度，从而根据信号到达时间的差异探测积冰，准确估计积冰层尺寸。而高频导波信号的模式转换过程更加复杂，在积冰区和出射区会出现多个模式，导致信号更加复杂。但部分信号点的模式转换过程比较稳定（如1239kHz-mode A1），同样可以用于积冰探测。4　时域仿真使用COMSOL®软件建立二维时域有限元模型，分析多种激励信号下的信号传播情况。其中，二维有限元模型的尺寸参数与频域有限元相似，如图7所示，包括1300mm×2mm铝板和0.5mm冰层，冰层长度可变，材料参数见表1。加载方式采用多点延时加载，加载点在激励器所在位置以线性阵列排列，间隔为波长的1/10，通过相位叠加产生单一模式的激励信号。10.12050/are20220316.F007图7时域有限元模型Fig.7Time domain finite element model在模型右端传感器处设置多个连续采样点，对于采样得到的多点时域信号，根据目标模式的相速度以及采样间隔计算时间延迟并叠加得到目标信号。加载点的激励信号表达式为Δt=Δxcpsi=sin(2πfct- iΔt)1- cos(2πfcNt- iΔt)      i=1,2,…,15 （6）式中：cp，fc，Δx，N，i分别为相速度中心频率、加载点间隔、周期数和加载点序号。为了提高激发效率，降低信号频率带宽，加载点个数设为15个，周期数为30。右端出射区获得采样信号，经过叠加得到最终结果Δt(n)=(n-1)dx/cpssuperposition(t)=∑n=1msreceving(n)(t-Δt(n)),    m=51 （7）式中：t为时间，ssuperposition和sreceiving(n)分别为相位叠加后的信号和第(n)个采样点接收到的信号，dx为采样间隔0.5mm，sn表示第n个采样点得到的时域信号。时域仿真中，采样区域为25mm。因此，最大采样点数为51pts。4.1　低频区域A图8(a)、图8(b)对应的激励信号分别为00kHz-mode S0 和400kHz-mode S0，两者的时域分析结果基本相同，即无冰情况下的单点接收信号，多点接收信号的相位叠加以及有冰情况下的叠加信号均存在一个主要波包，对应于mode S0。同时，由于导波在积冰区和无冰区的速度差异小，有冰情况下的信号到达时间并不明显。如图8(c)中，当激励信号为500kHz-mode S0 时，无冰情况下的单点接收信号包含两种模式A0和S0，这是导波的频散现象，通过多点接收信号的相位叠加可以消去不需要的模式，得到纯净的S0模式信号。积冰层出现后接收信号的S0模式波包出现微小的时间延迟，原因是500kHz导波信号在积冰铝板中（mode 2）的群速度略小于在铝板中（mode S0）的速度，但速度差异很小。如图8(d)所示，530kHz-mode S0激励信号得到的仿真结果与图8(c)非常相似，这也符合频域分析的结果。如图8(e)所示，当激励信号为854kHz-mode S0时，信号出现了明显变化，其中无冰情况下的单点接收信号出现多个小波包，但对主要波包没有明显影响，分别对应S0、A0和A1，符合频域有限元的计算结果（见图5(e)）。10.12050/are20220316.F008图8区域A中 5个典型信号的时域仿真结果Fig.8Time domain simulation results of five cases at region A4.2　高频区域B区域B中的激励频率较大，一般超过1MHz。根据上一节的分析也可以明显看出，该区域的导波传播过程更加复杂。如图9(a)所示，当激励信号为1101kHz-mode A1时，无冰情况下的单点接收信号只有一个主要波包，说明该频率下的激励效果很好，能量集中于所需的模式。结冰后的采样信号出现明显的时间延迟。同时，信号能量分布在多个小波包中导致主要波包的幅值明显降低，而经过信号叠加，小波包没有消失，说明这些小波包对应的模式与目标模式具有相同的相速度，即均为A1模式。如图9(b)所示，当激励信号为1178kHz-mode A1，激励效果同样较好，其他模式的波包经过相位叠加处理后完全消失。积冰层出现后，导波在积冰区出现多个模式，分别具有不同的群速度，经过积冰区后，导波信号转换回A1模式信号对应不同的到达时间，形成不同的波包，且不同波包之间相互叠加影响，使第一个波包的宽度明显增加，影响了对基于波峰位置的信号到达时间的判定。如图9(c)所示，当激励信号为1239kHz-mode A1，无冰信号与前述两种情况类似。有冰信号同样存在多余的小波包，但第一个主要波包的宽度未受到影响，因为该激励频率对应的导波在传播过程中的模式转换情况比较平滑（见图6(c)），主要能量均转化为单一的模式A1-mode 4-A1，其他模式包含能量较小，形成的波包可以忽略。通过上述时域结果可以看出，随着导波在积冰区和无冰区的群速度差异增加，积冰层对信号到达时间的影响更加显著。10.12050/are20220316.F009图9区域B中三个典型信号的时域仿真结果Fig.9Time domain simulation results of three cases at region B4.3　时域仿真结果通过对比有冰和无冰情况下叠加信号的到达时间（TOF）计算积冰尺寸。其中信号到达时间以目标模式所对应的波包的峰值时间为主。因此，首先利用希尔伯特变换（Hilbert Transform）对叠加信号取包络线，然后找到其中最大值所在的点对应的时间即为到达时间。最后始于以下公式计算积冰层长度ΔTOF=TOFicing-TOF0,  x=cgAlcgicingcgAl-cgicingΔTOF (8)式中：TOFicing，TOF0，cgicing，cgAl，x分别为有冰情况下的到达时间、无冰情况下的到达时间、积冰铝板的群速度，无冰铝板的群速度、积冰层长度。如图10所示为区域A、B中典型信号的时域仿真积冰尺寸探测结果，其中，黑色线条为积冰层的实际长度（0,100mm,200mm,300mm），与实际长度相比，区域A为低频信号，模式转换过程稳定，没有多余模式的干扰，所以探测结果精度普遍较高。而区域B为高频信号，模式转换过程复杂，多个不同速度的模式的出现对检测结果影响很大，导致探测精度降低，但部分信号点（如1239kHz-mode A1）的模式转换过程相对稳定，仍然能够保持良好的检测精度。10.12050/are20220316.F010图10时域仿真的积冰尺寸计算结果Fig.10Size evaluation results of ice10.12050/are20220316.T002表2比较不同激励信号平均相对误差Table 2Comparison of average relative error of differentexcitation signals激励频率/kHz平均相对误差/%300kHz-mode S08.27400kHz-mode S00.4500kHz-mode S03.90530kHz-mode S018.45854kHz-mode S06.811101kHz-mode A124.701178kHz-mode A144.401239kHz-mode A116.125　结论本文基于频域有限元仿真和波印亭矢量量化分析了超声导波在积冰铝板中的传播情况，根据群速度频散曲线的特征选定适合的信号区域，扫频分析域内各信号点的模式转换情况。使用时域有限元仿真分析8个典型激励信号的传播情况，并与频域分析结果相印证。结果发现，导波的模式转换对检测结果有显著影响，相比之下，能量集中，模式转换过程稳定的信号检测精度更高。与传统的基于幅值衰减的积冰探测方法相比，本方法基于群速度变化，具有更强的稳定性和鲁棒性。未来还能在此基础上进一步研究积冰层厚度、多点结冰等问题，实现更加全面可靠的积冰探测。