CN109885945A - 一种半空间环境下的边界元法近场声全息变换方法 - Google Patents

Abstract

本发明公开了一种半空间环境下的边界元法近场声全息变换方法，包括步骤一：建立半空间环境下的近场声全息测量模型，获取真实声源以及镜像虚源的源面上的节点坐标信息；步骤二：建立基于边界元理论的亥姆霍兹‑基尔霍夫(Helmholtz‑Kirchhoff)积分方程；步骤三：建立相应的传递矩阵；步骤四：建立全息变换关系式；步骤五：全息重建过程中的正则化处理。本发明考虑半空间测试环境，对传统半空间边界元法进行了改进，考虑了半空间中界面对源面声压及振速场的影响，建立了半空间环境下的全息变换矩阵，从而提高重构精度。

Description

一种半空间环境下的边界元法近场声全息变换方法

技术领域

本发明涉及一种半空间环境下的边界元法近场声全息变换方法，属于声全息领域中提高结构源面声压重构精度的变换方法。

背景技术

近场声全息技术已经被引入声源或振动体辐射场的测量与分析、噪声源识别定位及结构表面源强度评价等领域，成为当前声学领域一大研究热点。近场声全息技术的基本思想是对声源近场域内的复声压信息进行测量，并将测量得到的复声压信息按照声场变换算法进行重构，进而实现整个三维辐射声场中的声压、质点振速、声强和辐射声功率等声学量的预测，从而进行噪声源精确定位以及整个声场的可视化。目前，已发展了多项声全息空间变换技术，依据全息重建算法的类型主要分为两大类：共形面的声全息变换技术以及非共形面的声全息变换技术。其中，基于边界元法(BEM)的声全息变换技术凭借其在任意形面的全息变换技术中的重建优势，较为广泛的适用于实际工程应用中。其基本思想是将边界上的积分方程进行离散化处理，根据全息面与源面上声学参量的变换关系，建立全息变换矩阵，利用近场区域内的全息声压测试结果，实现声源表面的声学量重建，以及对源外部的声场预报。

在CNKI库中公开报道有：1.《半空间声全息法水下噪声源识别技术研究》，该论文建立了半空间近场声全息理论模型,提出了基于快速傅利叶变换的半空间声场全息重建技术。2.《运动声源的边界元声全息识别方法研究》，该文章主要提出了半空间格林函数，并进行了半空间环境中的边界元法非共形全息变换方法仿真及试验，验证了将声源关于水面镜面对称处理的有效性。3.《海底反射对声全息测试精度的影响研究》，针对声全息技术在实艇测试中存在的海底声反射问题，结合镜像法和半空间声全息理论，通过仿真计算，研究了某海区海底声反射对声全息测试结果的影响。以上文献中仅考虑了虚源对声场的影响，并没有考虑水面对结构的影响。在实际环境中，当声源近水面放置时，噪声源外部的辐射声场经水面反射后，会对源表面产生影响，因此，采用以上方法的重建结果必然存在较大误差。

发明内容

针对上述现有技术，本发明要解决的技术问题是提供一种考虑了半空间中界面对源面声压及振速场的影响，有效提高重构精度的半空间环境下的边界元法近场声全息变换方法。

为解决上述技术问题，本发明提供一种半空间环境下的边界元法近场声全息变换方法，包括以下步骤：

步骤一：建立半空间环境下的近场声全息测量模型，获取真实声源以及镜像虚源的源面上的节点坐标信息：

结合待分析的结构模型以及半空间界面位置，建立真实声源与镜像虚源的源面模型，并进行网格划分，确定真实声源以及镜像虚源的源面上的各节点的坐标信息，采用双平面作为全息测量面，建立源面模型，按照一个波长内不小于6个点的规则，采用8节点单元进进行网格划分；

步骤二：建立基于边界元理论的亥姆霍兹-基尔霍夫积分方程，具体为：

设弹性结构体处于无限流体介质中，S表示弹性结构体表面，D_i表示弹性结构体内部区域，弹性结构体外部区域的D_o充满密度为ρ，声速为c的流体，弹性结构体外部的场点表示为p，弹性结构体上的点表示为q，将弹性结构体表面边界S进行离散化处理，分为M个小单元，共N个节点，依据有限元和边界元理论，建立基于边界元理论的亥姆霍兹-基尔霍夫积分方程：

式中，分别表示从坐标中心指向对应场点的矢量，N_l(ξ)表示局部直角坐标系下单元的第l号节点处的插值形函数；J(ξ)表示雅克比系数，用于主坐标系和局部坐标系之间的变换；表示在m号单元l号节点上的声压值；表示在m号单元l号节点上的法向振速，代表格林函数，表示弹性结构体表面外法向偏导数，

当S为光滑面时，α(p)取值由相应场点所在的位置决定：

若S是非光滑面，当场点位于S面上时取值满足：

针对真实声源，表面振速声压与场点中声压之间的格林函数取自由场中格林函数，法向导数具体表达式满足：

虚拟声源表面振速声压与场点中声压之间满足半空间格林函数：

其中，是场点p与源面上点q的距离；是场点p与虚源面上点q'的距离；C_ri表示边界的反射系数，满足-1≤C_ri≤1，具体取值满足：

其中，θ表示对应声波的入射角；β表示反射界面的比导纳，设全息测量面上有i个测点，则对应测点的声压反射系数可表示为矩阵的形式C_r＝diag(C_r1,C_r2,···,C_ri)；

步骤三：建立相应的传递矩阵，具体为：

将节点进行统一的排序并相邻的重复节点进行合并，步骤二所述转换为：

其中：

假设在辐射声场中取L_o个测点、声源表面和声源内部分别取L、L_i个考察点，全息面上各个测点的声压用来表示，上式可分解成以下三个矩阵形式的方程：

[C]_L×N[P_s]_N×1＝[D]_L×N[V_n]_N×1 r∈S₀

其中，{P_H}表示D_o区域内某场点的声压；[P_s]表示弹性结构体源面上的声压；[V_n]表示弹性结构体源面上的法向振速；[P_n]表示D_i区域内场点的声压；相应的传递矩阵计算公式满足：

式中，当场点取节点i'处，i'与i重合时δ(i-i')取1；两者不重合时δ(i-i')取0；

步骤四：建立全息变换关系式，具体为：

将声源S关于对称界面进行对称，假设在对称位置处有一个与声源各处相位相反、幅值大小完全一致的“虚拟源”S*，虚源上的声压和振速分别为{P_S*}和{V_S*}，则考虑虚源影响下，真实的源面声压可表示为：

{P′_S}＝{P_S}+[W_SS*]{V_S*}＝[W_S]{V_S}+[W_SS*]{V_S*}

式中，{P_S}、{V_S}分别表示为理想声源的表面声压和振速；[W_S]表示理想声源表面声压与振速之间的传递矩阵，[W_SS*]表示虚拟声源表面振速与理想声源表面声压之间的传递矩阵，实际的辐射声场是虚拟声源和考虑虚源后的真实声源共同影响的叠加声场，则全息测量面实际声压为：

{P_h}＝{P_hd}+{P_hr}

其中，{P_hd}表示真实声源产生的直达声；{P_hr}表示虚拟声源产生的反射声；

将虚源影响下的真实源面声压代入源面振速计算公式可得：

{V′_S}＝[W_S]^-1{P′_S}＝{V_S}+[W_S]^-1[W_SS*]{V_S*}

式中，P_S'、V_S'分别表示虚源影响下真实源面的声压与振速，结合声压的反射系数，则全息面上实际测量的声压公式{P_h}＝{P_hd}+{P_hr}转化为下式：

{P_h}＝[W_d]{V_S'}+[C_r][W_r]{V_S}

式中，[W_d]代表真实源面在全息面上产生的辐射声压和真实源面振速的传递矩阵，[W_r]代表虚拟源面在全息面上产生的辐射声压和虚拟源面振速的传递矩阵；

将{V′_S}＝[W_S]^-1{P′_S}＝{V_S}+[W_S]^-1[W_SS*]{V_S*}代入{P_h}＝[W_d]{V_S'}+[C_r][W_r]{V_S}中，{P_h}＝[W_d]{V_S'}+[C_r][W_r]{V_S}可转化为：

{P_h}＝[W_d]{V_S'}+[C_r][W_r]{V_S}＝{[W_d]+[C_r][W_r]+[W_d][W_S]^-1[W_SS*]}{V_S}＝{W_sh}{V_S}

其中{W_sh}＝{[W_d]+[C_r][W_r]+[W_d][W_S]^-1[W_SS*]}即为半空间的源面到全息面的传递矩阵；

将上式中的传递矩阵采用奇异分解法进行求逆计算，得到半空间边界元法声全息的声源振速重建公式：

{V_S}＝{W_sh}^-1{P_h}＝[U][∑]^-1[V]^H{P_h}

步骤五：全息重建过程中的正则化处理，具体为：

采用Tikhonov正则化方法，通过施加残余范数和单边约束形成的加权组合对{V_S}＝{W_sh}^-1{P_h}＝[U][∑]^-1[V]^H{P_h}取最小的约束条件，即{V_S}_reg满足以下条件：

其中，λ取值大于零，代表正则化参数；(V_S)^*代表源面法向振速的初始估计值；L代表罚矩阵；

对进行求解，得到广义Tikhonov正则化的解，满足：

{V_S}_reg＝(W^HW+λ²L^TL)^-1W^H{P_h}

由上式可知，利用(W^HW+λ²L^TL)^-1W^H取代源面振速重建公式中的病态矩阵的逆W^-1；

当L＝I_n，(V_S)^*＝0时，即标准正则化方法，源面的法向振速可表示为：

标准Tikhonov正则化方法的滤波系数表达式满足：

正则化参数λ的选取方法为L曲线准则，L曲线的拐角处值为正则化参数。

当L≠I_n，且为可逆方阵时，将一般形式Tikhonov正则化化为标准形式，首先将{V_S}、W逆进行转换，得到新的形式将转化为标准形式：

标准正则化下的源面法向振速通过将转换后的矩阵进行奇异值分解求得，并进行转换，即可得到非标准形式的源面法向振速

作为本发明的一种优选方案，步骤五中L为正定矩阵。

作为本发明另一种优选方案，步骤五中正则化参数λ的选取方法为广义交叉检验法，具体为：

正则化参数λ通过方程确定，其中，W^I代表正则化解{V_S}_reg对应的传递矩阵；trace()指传递矩阵的迹；当方程求最小值时，对应的λ即为所求的正则化参数。

本发明有益效果：提出半空间环境下的边界元法近场声全息变换方法，在全息变换时，将镜像虚源看作为声源，将真实声源的源面看作为声场中的一个封闭包面，将由虚源产生的在源面所在位置处的声学量与声源本身源面上的声学量相叠加，得到在半空间环境下，受到界面干涉后，源表面上的振速及声压，进而可以得到源面上实际上由真实声源产生的声场，减少界面反射对源表面产生的影响。本发明考虑半空间测试环境，对传统半空间边界元法进行了改进，考虑了半空间中界面对源面声压及振速场的影响，建立了半空间环境下的全息变换矩阵，从而提高重构精度。

附图说明

图1是本发明的方法流程图；

图2为本发明步骤1中建立的近场声全息测量模型；

图3为本发明步骤2中所建立的无限理想流体介质弹性结构体模型示意图；

图4为本发明步骤4中的半空间辐射声场示意图；

图5为本发明步骤5中的L曲线示意图；

图6为本发明实施例1中的双平面全息面示意图；

图7(a)为本发明实施例1中采用传统半空间全息变换方法得到的球壳表面声压幅值重构结果；

图7(b)为本发明实施例1中采用发明中改进的半空间变换方法得到的球壳表面声压幅值重构结果；

图7(c)为本发明实施例1中球壳表面声压幅值的理论值；

图8(a)为本发明实施例1中采用本发明方法得到的界面距离模型中心0.8m时源面声压全息重构误差；

图8(b)为本发明实施例1中采用本发明方法得到的界面距离模型中心3.0m时源面声压全息重构误差；

图8(c)为本发明实施例1中采用本发明方法得到的界面距离模型中心15m时源面声压全息重构误差；

图8(d)为本发明实施例1中采用本发明方法得到的界面距离模型中心37.5m时源面声压全息重构误差；

图9为本发明实施例2中壳体实物及其内部激励源；

图10(a)为本发明实施例2中传统半空间变换方法源面声压的重构结果；

图10(b)为本发明实施例2本发明中方法源面声压的重构结果。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明是一种半空间环境下的边界元法近场声全息变换方法，方法流程如图1所示，包括以下步骤：

步骤一：建立半空间环境下的近场声全息测量模型，获取真实声源以及镜像虚源的源面上的节点坐标信息；

结合待分析的水下结构模型以及半空间界面位置，建立真实声源与镜像虚源的源面模型，并进行网格划分，确定真实声源以及镜像虚源的源面上的各节点的坐标信息。为提高精度，采用双平面作为全息测量面，以球壳模型为例，所建立的近场声全息测量模型，如图2所示；建立网格模型时，可以采用边界元软件ANSYS建立源面模型，并按照一个波长内不小于6个点的规则，采用8节点单元进进行网格划分。

步骤二：建立基于边界元理论的亥姆霍兹-基尔霍夫(Helmholtz-Kirchhoff)积分方程

设弹性结构体处于无限流体介质中，S表示弹性结构体表面，D_i表示弹性结构体内部区域，弹性结构体外部区域的D_o充满密度为ρ，声速为c的流体，弹性结构体外部的场点表示为p，弹性结构体上的点表示为q，如图3所示。

将弹性结构体表面边界S进行离散化处理，分为M个小单元，共N个节点，依据有限元和边界元理论，建立基于边界元理论的亥姆霍兹-基尔霍夫(Helmholtz-Kirchhoff)积分方程：

式中，分别表示从坐标中心指向对应场点的矢量，N_l(ξ)表示局部直角坐标系下单元的第l号节点处的插值形函数；J(ξ)表示雅克比系数，用于主坐标系和局部坐标系之间的变换；表示在m号单元l号节点上的声压值；表示在m号单元l号节点上的法向振速。代表格林函数，表示弹性结构体表面外法向偏导数。

当S为光滑面时，α(p)取值由相应场点所在的位置决定：

若S是非光滑面，当场点位于S面上时取值为：

针对真实声源，其表面振速声压与场点中声压之间的格林函数应取自由场中格林函数，其法向导数具体表达式如下：

由于在半空间环境下界面的反射效应，存在着镜像虚源，对于虚拟声源而言，其表面振速声压与场点中声压之间应满足半空间格林函数：

其中，是场点p与源面上点q的距离；是场点p与虚源面上点q'的距离；C_ri表示边界的反射系数，满足-1≤C_ri≤1，具体取值如下：

其中，θ表示对应声波的入射角；β表示反射界面的比导纳。因此，设全息测量面上有i个测点，则对应测点的声压反射系数可表示为矩阵的形式C_r＝diag(C_r1,C_r2,···,C_ri)。

步骤三：建立相应的传递矩阵；

由于每个单元与其相邻的单元有共同的节点，将节点进行统一的排序并相邻的重复节点进行合并，(1)式转换为：

其中：

假设在辐射声场中取L_o个测点、声源表面和声源内部分别取L、L_i个考察点，全息面上各个测点的声压用来表示，上式可分解成以下三个矩阵形式的方程：

[C]_L×N[P_s]_N×1＝[D]_L×N[V_n]_N×1 r∈S₀ (11)

其中，{P_H}表示D_o区域内某场点的声压；[P_s]表示弹性结构体源面上的声压；[V_n]表示弹性结构体源面上的法向振速；[P_n]表示D_i区域内场点的声压；相应的传递矩阵计算公式如下：

式中，当场点取节点i'处，i'与i重合时δ(i-i')取1；两者不重合时δ(i-i')取0。

利用以上计算方法，即可根据不同位置关系，建立不同的传递矩阵，具体包括：(1)建立实际声源表面振速与表面声压之间的传递矩阵、镜像虚源表面振速与实际声源表面声压之间的传递矩阵；(2)结合各全息测点与真实声源的位置关系，建立实际声源表面振速与全息声压之间的传递矩阵；(3)结合各全息测点与镜像虚源的位置关系，建立虚拟声源的表面振速与全息声压之间的传递矩阵。

步骤四：建立全息变换关系式

在半空间环境中，全息面上的测量结果为实际声源产生的辐射声压与界面的反射声的叠加结果，采用镜像法，将界面反射声等效为对称虚源的辐射噪声，从而实现半空间中的声源重建。根据镜像原理，虚源上任意节点位置分布情况都和实际声源完全一致，而且虚源源面上每个节点的表面法向振速和实际源面的表面法向振速对应相等；考虑水面对声源的影响，实际声源应看作理想声源和虚源对理想声源的作用两者叠加而成。

根据镜像法的原理建立半空间辐射声场理论模型，如图4所示。以水面为对称界面，将声源S关于对称界面进行对称，假设在对称位置处有一个与声源各处相位相反、幅值大小完全一致的“虚拟源”S*，虚源上的声压和振速分别为{P_S*}和{V_S*}，则考虑虚源影响下，真实的源面声压可表示为：

{P′_S}＝{P_S}+[W_SS*]{V_S*}＝[W_S]{V_S}+[W_SS*]{V_S*} (19)

式中，{P_S}、{V_S}分别表示为理想声源的表面声压和振速；[W_S]表示理想声源表面声压与振速之间的传递矩阵。[W_SS*]表示虚拟声源表面振速与理想声源表面声压之间的传递矩阵。而此时实际的辐射声场是虚拟声源和考虑虚源后的真实声源共同影响的叠加声场。因此，全息测量面实际声压为：

{P_h}＝{P_hd}+{P_hr} (20)

其中，{P_hd}表示真实声源产生的直达声；{P_hr}表示虚拟声源产生的反射声。

将虚源影响下的真实源面声压代入源面振速计算公式(9)可得：

{V′_S}＝[W_S]^-1{P′_S}＝{V_S}+[W_S]^-1[W_SS*]{V_S*} (21)

式中，P_S'、V_S'分别表示虚源影响下真实源面的声压与振速。结合声压的反射系数，则全息面上实际测量的声压公式(20)转化为下式：

{P_h}＝[W_d]{V_S'}+[C_r][W_r]{V_S} (22)

式中，[W_d]代表真实源面在全息面上产生的辐射声压和真实源面振速的传递矩阵，[W_r]代表虚拟源面在全息面上产生的辐射声压和虚拟源面振速的传递矩阵。

将式(21)代入(22)中，上式即可转化为：

{P_h}＝[W_d]{V_S'}+[C_r][W_r]{V_S}＝{[W_d]+[C_r][W_r]+[W_d][W_S]^-1[W_SS*]}{V_S}＝{W_sh}{V_S} (23)

其中{W_sh}＝{[W_d]+[C_r][W_r]+[W_d][W_S]^-1[W_SS*]}即为半空间的源面到全息面的传递矩阵。

如上，即可实现半空间环境下的全息重建。将上式中的传递矩阵采用奇异分解法进行求逆计算，得到：

{V_S}＝{W_sh}^-1{P_h}＝[U][∑]^-1[V]^H{P_h} (24)

上式即为半空间边界元法声全息的声源振速重建公式。

步骤五：全息重建过程中的正则化处理

当求解病态矩阵过程中产生大量的离散值时，正则化处理较为适用就是Tikhonov正则化方法。通过原有方程中，施加残余范数和单边约束形成的加权组合取最小的约束条件，即{V_S}_reg须满足以下条件：

其中，λ取值大于零，代表正则化参数；(V_S)^*代表源面法向振速的初始估计值；L代表罚矩阵，一般为正定矩阵，相当于要求未知解满足一定的光滑性约束条件。

对上述方程进行求解，得到广义Tikhonov正则化的解，

{V_S}_reg＝(W^HW+λ²L^TL)^-1W^H{P_h} (26)

由上式可知，利用(W^HW+λ²L^TL)^-1W^H取代源面振速重建公式中的病态矩阵的逆W^-1。

当L＝I_n，(V_S)^*＝0时，即标准正则化方法，其源面的法向振速可表示为：

为了减少测量误差的影响，正则化系数的取值尤为关键，以下给出标准Tikhonov正则化方法的滤波系数表达式：

其中正则化参数λ最为常见的选取方法包括广义交叉检验法和L曲线准则。广义交叉验证法主要通过下述方程来确定：

其中，W^I代表正则化解{V_S}_reg对应的传递矩阵；trace()指传递矩阵的迹；可以看出，当方程求最小值时，对应的λ即为所求的正则化参数。

以L曲线为基础的L曲线准则(L-Curve Criterion)与广义交叉验证法相同，对误差范数是否已知不作要求。其中，L曲线描述的是利用所有可行的正则化参数λ计算并以对数尺度绘制残余范数与解范数的关系。然后，通过对残余范数与解范数的对比选取有效的正则化参数。图5为L曲线示意图，图中拐角的左半边部分，由于λ太小，出现欠正则化现象，虽然采用正则化处理，但仍存在某些奇异值较小的分量导致未对误差实现有效的控制；而在L曲线拐角的右半边部分，由于λ太大，出现过正则化现象，滤波后的奇异值太少导致有效信息的丢失。因此，L曲线的拐角处才是能有效实现滤波的正则化参数，并且L曲线越是满足Picard条件，其拐角越突出。

当L≠I_n，且为可逆方阵时，需要将一般形式Tikhonov正则化化为标准形式。首先将{V_S}、W逆进行转换，得到新的形式并代入原来方程中，将式(25)转化为标准形式：

标准正则化下的源面法向振速通过将转换后的矩阵进行奇异值分解求得，并进行转换，即可得到非标准形式的源面法向振速

实施例：

对本发明的有益效果如下方式得以验证：

(1)在Matlab仿真条件下，对该方法进行仿真实验：

仿真参数如下：

模型几何参数：球壳，半径0.75m，令中心位于坐标原点，划分网格后，形成218个节点及216个四边形单元；

模型内部声源参数：位于(-0.1,-0.1,0.1)，点声源，幅值为1Pa；

流体介质参数：水，密度1000kg/m³,声速1500m/s

半空间界面：模型处于水中，半空间界面无限大，界面外为空气，可认为界面反射系数R＝-1，界面距离球壳中心Hm；

全息测量面参数：双矩形全息面，尺寸为4m×4m，测点间隔为0.2m，距离球壳中心为0.8m；

计算频率：f＝200Hz。

仿真模型如图6所示。

利用发明所述方法，得到球壳表面声压幅值重构结果如图7(a)-图7(c)所示，其中图7(a)表示传统半空间全息变换方法的球壳表面声压幅值重构结果，图7(b)表示本发明中改进的半空间变换方法的重构结果，图7(c)表示球壳表面声压幅值的理论值。结果表明本发明中考虑了水面干扰以后的全息变换算法则相对精度较高。

图8(a)-图8(d)所示为模型位于不同水深情况下的源面声压全息重构误差，其中图8(a)、图8(b)、图8(c)、图8(d)分别对应界面距离模型中心0.8m、3.0m、15m、37.5m时的源面声压幅值重构误差。结果表明采用本发明中的变换方法，入水深度H越大时，改进后的半空间全息变换方法与传统的半空间边界元变换方法的重构精度相差越小，这也符合入水深度越大，水面对源面的干扰越小的实际情况。而当模型离水面较近时，采用改进的变换方法可以大幅度提高重构精度。

(2)水池中的近场声全息变换方法试验验证

试验模型：单层圆柱壳体，半径为0.2m，壁厚0.005m，长度0.8m，壳体实物及其内部激励源见图9，激励源距离模型上端面约为25cm。

信号发射系统：由信号源(Agilent33522A)控制圆柱壳体内的激励源，经功率放大器(YF5887)实现单频信号激励，频率2000Hz。

水听器接收系统：30个水听器(B&K8103)组成一条垂直阵，阵元间距为0.06m，阵长1.87m，水听器接收信号由数据采集器(B&K3660D)采集。

圆柱壳体竖直放入水中，上表面距离水面为3.5cm，水听器阵竖直放入水中，最上端的水听器距离水面3cm。利用实验水池中的升降回转装置，实现两个1.87m×1.87m的矩形平面的全息扫描。

利用全息测试数据对模型表面的辐射声压进行重构，采用奇异值截断法进行正则化法处理，将信噪比设定为20dB。图10(a)和图10(b)为源面声压的重构结果，分别为传统半空间变换方法与采用本发明中方法的重构结果。结果表明采用本发明中的半空间全息变换模型，可以明显的提高全息重构精度，可以清晰的看出源所在的位置，与实际模型中激励源位置吻合。

本发明具体实施方式还包括：

本发明方法包括以下步骤：

步骤一：建立半空间环境下的近场声全息测量模型，获取真实声源以及镜像虚源的源面上的节点坐标信息；

步骤二：建立基于边界元理论的亥姆霍兹-基尔霍夫(Helmholtz-Kirchhoff)积分方程；

步骤三：建立相应的传递矩阵；

步骤四：建立全息变换关系式；

步骤五：全息重建过程中的正则化处理。

Claims (3)

1.一种半空间环境下的边界元法近场声全息变换方法，其特征在于，包括以下步骤：

步骤一：建立半空间环境下的近场声全息测量模型，获取真实声源以及镜像虚源的源面上的节点坐标信息：

结合待分析的结构模型以及半空间界面位置，建立真实声源与镜像虚源的源面模型，并进行网格划分，确定真实声源以及镜像虚源的源面上的各节点的坐标信息，采用双平面作为全息测量面，建立源面模型，按照一个波长内不小于6个点的规则，采用8节点单元进进行网格划分；

步骤二：建立基于边界元理论的亥姆霍兹-基尔霍夫积分方程，具体为：

设弹性结构体处于无限流体介质中，S表示弹性结构体表面，D_i表示弹性结构体内部区域，弹性结构体外部区域的D_o充满密度为ρ，声速为c的流体，弹性结构体外部的场点表示为p，弹性结构体上的点表示为q，将弹性结构体表面边界S进行离散化处理，分为M个小单元，共N个节点，依据有限元和边界元理论，建立基于边界元理论的亥姆霍兹-基尔霍夫积分方程：

式中，分别表示从坐标中心指向对应场点的矢量，N_l(ξ)表示局部直角坐标系下单元的第l号节点处的插值形函数；J(ξ)表示雅克比系数，用于主坐标系和局部坐标系之间的变换；表示在m号单元l号节点上的声压值；表示在m号单元l号节点上的法向振速，

代表格林函数，表示弹性结构体表面外法向偏导数，

当S为光滑面时，α(p)取值由相应场点所在的位置决定：

若S是非光滑面，当场点位于S面上时取值满足：

针对真实声源，表面振速声压与场点中声压之间的格林函数取自由场中格林函数，法向导数具体表达式满足：

虚拟声源表面振速声压与场点中声压之间满足半空间格林函数：

其中，是场点p与源面上点q的距离；是场点p与虚源面上点q'的距离；C_ri表示边界的反射系数，满足-1≤C_ri≤1，具体取值满足：

其中，θ表示对应声波的入射角；β表示反射界面的比导纳，设全息测量面上有i个测点，则对应测点的声压反射系数可表示为矩阵的形式C_r＝diag(C_r1,C_r2,···,C_ri)；

步骤三：建立相应的传递矩阵，具体为：

将节点进行统一的排序并相邻的重复节点进行合并，步骤二所述转换为：

其中：

假设在辐射声场中取L_o个测点、声源表面和声源内部分别取L、L_i个考察点，全息面上各个测点的声压用来表示，上式可分解成以下三个矩阵形式的方程：

[C]_L×N[P_s]_N×1＝[D]_L×N[V_n]_N×1 r∈S₀

其中，{P_H}表示D_o区域内某场点的声压；[P_s]表示弹性结构体源面上的声压；[V_n]表示弹性结构体源面上的法向振速；[P_n]表示D_i区域内场点的声压；相应的传递矩阵计算公式满足：

式中，当场点取节点i'处，i'与i重合时δ(i-i')取1；两者不重合时δ(i-i')取0；

步骤四：建立全息变换关系式，具体为：

将声源S关于对称界面进行对称，假设在对称位置处有一个与声源各处相位相反、幅值大小完全一致的“虚拟源”S*，虚源上的声压和振速分别为{P_S*}和{V_S*}，则考虑虚源影响下，真实的源面声压可表示为：

式中，{P_S}、{V_S}分别表示为理想声源的表面声压和振速；[W_S]表示理想声源表面声压与振速之间的传递矩阵，表示虚拟声源表面振速与理想声源表面声压之间的传递矩阵，实际的辐射声场是虚拟声源和考虑虚源后的真实声源共同影响的叠加声场，则全息测量面实际声压为：

{P_h}＝{P_hd}+{P_hr}

其中，{P_hd}表示真实声源产生的直达声；{P_hr}表示虚拟声源产生的反射声；

将虚源影响下的真实源面声压代入源面振速计算公式可得：

式中，P_S'、V_S'分别表示虚源影响下真实源面的声压与振速，结合声压的反射系数，则全息面上实际测量的声压公式{P_h}＝{P_hd}+{P_hr}转化为下式：

{P_h}＝[W_d]{V_S'}+[C_r][W_r]{V_S}

式中，[W_d]代表真实源面在全息面上产生的辐射声压和真实源面振速的传递矩阵，[W_r]代表虚拟源面在全息面上产生的辐射声压和虚拟源面振速的传递矩阵；

将代入{P_h}＝[W_d]{V_S'}+[C_r][W_r]{V_S}中，{P_h}＝[W_d]{V_S'}+[C_r][W_r]{V_S}可转化为：

其中即为半空间的源面到全息面的传递矩阵；

将上式中的传递矩阵采用奇异分解法进行求逆计算，得到半空间边界元法声全息的声源振速重建公式：

{V_S}＝{W_sh}^-1{P_h}＝[U][∑]^-1[V]^H{P_h}

步骤五：全息重建过程中的正则化处理，具体为：

采用Tikhonov正则化方法，通过施加残余范数和单边约束形成的加权组合对{V_S}＝{W_sh}^-1{P_h}＝[U][∑]^-1[V]^H{P_h}取最小的约束条件，即{V_S}_reg满足以下条件：

其中，λ取值大于零，代表正则化参数；(V_S)^*代表源面法向振速的初始估计值；L代表罚矩阵；

对进行求解，得到广义Tikhonov正则化的解，满足：

{V_S}_reg＝(W^HW+λ²L^TL)^-1W^H{P_h}

由上式可知，利用(W^HW+λ²L^TL)^-1W^H取代源面振速重建公式中的病态矩阵的逆W^-1；

当L＝I_n，(V_S)^*＝0时，即标准正则化方法，源面的法向振速可表示为：

标准Tikhonov正则化方法的滤波系数表达式满足：

正则化参数λ的选取方法为L曲线准则，L曲线的拐角处值为正则化参数。

当L≠I_n，且为可逆方阵时，将一般形式Tikhonov正则化化为标准形式，首先将{V_S}、W逆进行转换，得到新的形式将转化为标准形式：

标准正则化下的源面法向振速通过将转换后的矩阵进行奇异值分解求得，并进行转换，即可得到非标准形式的源面法向振速

2.根据权利要求1所述的一种半空间环境下的边界元法近场声全息变换方法，其特征在于：步骤五所述L为正定矩阵。

3.根据权利要求1所述的一种半空间环境下的边界元法近场声全息变换方法，其特征在于：步骤五所述正则化参数λ的选取方法为广义交叉检验法，具体为：

正则化参数λ通过方程确定，其中，W^I代表正则化解{V_S}_reg对应的传递矩阵；trace()指传递矩阵的迹；当方程求最小值时，对应的λ即为所求的正则化参数。