[go: up one dir, main page]
More Web Proxy on the site http://driver.im/

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

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

Info

Publication number
CN109885945B
CN109885945B CN201910141127.8A CN201910141127A CN109885945B CN 109885945 B CN109885945 B CN 109885945B CN 201910141127 A CN201910141127 A CN 201910141127A CN 109885945 B CN109885945 B CN 109885945B
Authority
CN
China
Prior art keywords
source
sound
sound pressure
holographic
field
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.)
Active
Application number
CN201910141127.8A
Other languages
English (en)
Other versions
CN109885945A (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.)
Harbin Engineering University
Original Assignee
Harbin Engineering University
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 Harbin Engineering University filed Critical Harbin Engineering University
Priority to CN201910141127.8A priority Critical patent/CN109885945B/zh
Publication of CN109885945A publication Critical patent/CN109885945A/zh
Application granted granted Critical
Publication of CN109885945B publication Critical patent/CN109885945B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Landscapes

  • Measurement Of Mechanical Vibrations Or Ultrasonic Waves (AREA)
  • Holo Graphy (AREA)

Abstract

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

Description

一种半空间环境下的边界元法近场声全息变换方法
技术领域
本发明涉及一种半空间环境下的边界元法近场声全息变换方法，属于声全息领域中提高结构源面声压重构精度的变换方法。
背景技术
近场声全息技术已经被引入声源或振动体辐射场的测量与分析、噪声源识别定位及结构表面源强度评价等领域，成为当前声学领域一大研究热点。近场声全息技术的基本思想是对声源近场域内的复声压信息进行测量，并将测量得到的复声压信息按照声场变换算法进行重构，进而实现整个三维辐射声场中的声压、质点振速、声强和辐射声功率等声学量的预测，从而进行噪声源精确定位以及整个声场的可视化。目前，已发展了多项声全息空间变换技术，依据全息重建算法的类型主要分为两大类：共形面的声全息变换技术以及非共形面的声全息变换技术。其中，基于边界元法(BEM)的声全息变换技术凭借其在任意形面的全息变换技术中的重建优势，较为广泛的适用于实际工程应用中。其基本思想是将边界上的积分方程进行离散化处理，根据全息面与源面上声学参量的变换关系，建立全息变换矩阵，利用近场区域内的全息声压测试结果，实现声源表面的声学量重建，以及对源外部的声场预报。
在CNKI库中公开报道有：1.《半空间声全息法水下噪声源识别技术研究》，该论文建立了半空间近场声全息理论模型,提出了基于快速傅利叶变换的半空间声场全息重建技术。2.《运动声源的边界元声全息识别方法研究》，该文章主要提出了半空间格林函数，并进行了半空间环境中的边界元法非共形全息变换方法仿真及试验，验证了将声源关于水面镜面对称处理的有效性。3.《海底反射对声全息测试精度的影响研究》，针对声全息技术在实艇测试中存在的海底声反射问题，结合镜像法和半空间声全息理论，通过仿真计算，研究了某海区海底声反射对声全息测试结果的影响。以上文献中仅考虑了虚源对声场的影响，并没有考虑水面对结构的影响。在实际环境中，当声源近水面放置时，噪声源外部的辐射声场经水面反射后，会对源表面产生影响，因此，采用以上方法的重建结果必然存在较大误差。
发明内容
针对上述现有技术，本发明要解决的技术问题是提供一种考虑了半空间中界面对源面声压及振速场的影响，有效提高重构精度的半空间环境下的边界元法近场声全息变换方法。
为解决上述技术问题，本发明提供一种半空间环境下的边界元法近场声全息变换方法，包括以下步骤：
步骤一：建立半空间环境下的近场声全息测量模型，获取真实声源以及镜像虚源的源面上的节点坐标信息：
结合待分析的结构模型以及半空间界面位置，建立真实声源与镜像虚源的源面模型，并进行网格划分，确定真实声源以及镜像虚源的源面上的各节点的坐标信息，采用双平面作为全息测量面，建立源面模型，按照一个波长内不小于6个点的规则，采用8节点单元进进行网格划分；
步骤二：建立基于边界元理论的亥姆霍兹-基尔霍夫积分方程，具体为：
设弹性结构体处于无限流体介质中，S表示弹性结构体表面，Di表示弹性结构体内部区域，弹性结构体外部区域的Do充满密度为ρ，声速为c的流体，弹性结构体外部的场点表示为p，弹性结构体上的点表示为q，将弹性结构体表面边界S进行离散化处理，分为M个小单元，共N个节点，依据有限元和边界元理论，建立基于边界元理论的亥姆霍兹-基尔霍夫积分方程：
Figure BDA0001978583220000021
式中，
Figure BDA0001978583220000022
分别表示从坐标中心指向对应场点的矢量，Nl(ξ)表示局部直角坐标系下单元的第l号节点处的插值形函数；J(ξ)表示雅克比系数，用于主坐标系和局部坐标系之间的变换；
Figure BDA0001978583220000023
表示在m号单元l号节点上的声压值；
Figure BDA0001978583220000024
表示在m号单元l号节点上的法向振速，
Figure BDA0001978583220000025
代表格林函数，
Figure BDA0001978583220000026
表示弹性结构体表面外法向偏导数，
当S为光滑面时，α(p)取值由相应场点所在的位置决定：
Figure BDA0001978583220000027
若S是非光滑面，当场点位于S面上时取值满足：
Figure BDA0001978583220000028
针对真实声源，表面振速声压与场点中声压之间的格林函数取自由场中格林函数，法向导数具体表达式满足：
Figure BDA0001978583220000029
Figure BDA0001978583220000031
虚拟声源表面振速声压与场点中声压之间满足半空间格林函数：
Figure BDA0001978583220000032
其中，
Figure BDA0001978583220000033
是场点p与源面上点q的距离；
Figure BDA0001978583220000034
是场点p与虚源面上点q'的距离；Cri表示边界的反射系数，满足-1≤Cri≤1，具体取值满足：
Figure BDA0001978583220000035
其中，θ表示对应声波的入射角；β表示反射界面的比导纳，设全息测量面上有i个测点，则对应测点的声压反射系数可表示为矩阵的形式Cr＝diag(Cr1,Cr2,···,Cri)；
步骤三：建立相应的传递矩阵，具体为：
将节点进行统一的排序并相邻的重复节点进行合并，步骤二所述
Figure BDA0001978583220000036
转换为：
Figure BDA0001978583220000037
其中：
Figure BDA0001978583220000038
假设在辐射声场中取Lo个测点、声源表面和声源内部分别取L、Li个考察点，全息面上各个测点的声压用
Figure BDA0001978583220000039
来表示，
Figure BDA00019785832200000310
上式可分解成以下三个矩阵形式的方程：
Figure BDA00019785832200000311
[C]L×N[Ps]N×1＝[D]L×N[Vn]N×1 r∈S0
Figure BDA0001978583220000041
其中，{PH}表示Do区域内某场点的声压；[Ps]表示弹性结构体源面上的声压；[Vn]表示弹性结构体源面上的法向振速；[Pn]表示Di区域内场点的声压；相应的传递矩阵计算公式满足：
Figure BDA0001978583220000042
Figure BDA0001978583220000043
Figure BDA0001978583220000044
式中，当场点取节点i'处，i'与i重合时δ(i-i')取1；两者不重合时δ(i-i')取0；
Figure BDA0001978583220000045
Figure BDA0001978583220000046
Figure BDA0001978583220000047
步骤四：建立全息变换关系式，具体为：
将声源S关于对称界面进行对称，假设在对称位置处有一个与声源各处相位相反、幅值大小完全一致的“虚拟源”S*，虚源上的声压和振速分别为{PS*}和{VS*}，则考虑虚源影响下，真实的源面声压可表示为：
{P′S}＝{PS}+[WSS*]{VS*}＝[WS]{VS}+[WSS*]{VS*}
式中，{PS}、{VS}分别表示为理想声源的表面声压和振速；[WS]表示理想声源表面声压与振速之间的传递矩阵，[WSS*]表示虚拟声源表面振速与理想声源表面声压之间的传递矩阵，实际的辐射声场是虚拟声源和考虑虚源后的真实声源共同影响的叠加声场，则全息测量面实际声压为：
{Ph}＝{Phd}+{Phr}
其中，{Phd}表示真实声源产生的直达声；{Phr}表示虚拟声源产生的反射声；
将虚源影响下的真实源面声压代入源面振速计算公式
Figure BDA0001978583220000051
可得：
{V′S}＝[WS]-1{P′S}＝{VS}+[WS]-1[WSS*]{VS*}
式中，PS'、VS'分别表示虚源影响下真实源面的声压与振速，结合声压的反射系数，则全息面上实际测量的声压公式{Ph}＝{Phd}+{Phr}转化为下式：
{Ph}＝[Wd]{VS'}+[Cr][Wr]{VS}
式中，[Wd]代表真实源面在全息面上产生的辐射声压和真实源面振速的传递矩阵，[Wr]代表虚拟源面在全息面上产生的辐射声压和虚拟源面振速的传递矩阵；
将{V′S}＝[WS]-1{P′S}＝{VS}+[WS]-1[WSS*]{VS*}代入{Ph}＝[Wd]{VS'}+[Cr][Wr]{VS}中，{Ph}＝[Wd]{VS'}+[Cr][Wr]{VS}可转化为：
{Ph}＝[Wd]{VS'}+[Cr][Wr]{VS}＝{[Wd]+[Cr][Wr]+[Wd][WS]-1[WSS*]}{VS}＝{Wsh}{VS}
其中{Wsh}＝{[Wd]+[Cr][Wr]+[Wd][WS]-1[WSS*]}即为半空间的源面到全息面的传递矩阵；
将上式中的传递矩阵采用奇异分解法进行求逆计算，得到半空间边界元法声全息的声源振速重建公式：
{VS}＝{Wsh}-1{Ph}＝[U][∑]-1[V]H{Ph}
步骤五：全息重建过程中的正则化处理，具体为：
采用Tikhonov正则化方法，通过施加残余范数和单边约束形成的加权组合对{VS}＝{Wsh}-1{Ph}＝[U][∑]-1[V]H{Ph}取最小的约束条件，即{VS}reg满足以下条件：
Figure BDA0001978583220000052
其中，λ取值大于零，代表正则化参数；(VS)*代表源面法向振速的初始估计值；L代表罚矩阵；
对
Figure BDA0001978583220000053
进行求解，得到广义Tikhonov正则化的解，满足：
{VS}reg＝(WHW+λ2LTL)-1WH{Ph}
由上式可知，利用(WHW+λ2LTL)-1WH取代源面振速重建公式中的病态矩阵的逆W-1；
当L＝In，(VS)*＝0时，即标准正则化方法，源面的法向振速可表示为：
Figure BDA0001978583220000061
标准Tikhonov正则化方法的滤波系数表达式满足：
Figure BDA0001978583220000062
正则化参数λ的选取方法为L曲线准则，L曲线的拐角处值为正则化参数。
当L≠In，且为可逆方阵时，将一般形式Tikhonov正则化化为标准形式，首先将{VS}、W逆进行转换，得到新的形式
Figure BDA0001978583220000063
将
Figure BDA0001978583220000064
转化为标准形式：
Figure BDA0001978583220000065
标准正则化下的源面法向振速通过将转换后的矩阵进行奇异值分解求得，并进行转换，即可得到非标准形式的源面法向振速
Figure BDA0001978583220000066
作为本发明的一种优选方案，步骤五中L为正定矩阵。
作为本发明另一种优选方案，步骤五中正则化参数λ的选取方法为广义交叉检验法，具体为：
正则化参数λ通过方程
Figure BDA0001978583220000067
确定，其中，WI代表正则化解{VS}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表示弹性结构体表面，Di表示弹性结构体内部区域，弹性结构体外部区域的Do充满密度为ρ，声速为c的流体，弹性结构体外部的场点表示为p，弹性结构体上的点表示为q，如图3所示。
将弹性结构体表面边界S进行离散化处理，分为M个小单元，共N个节点，依据有限元和边界元理论，建立基于边界元理论的亥姆霍兹-基尔霍夫(Helmholtz-Kirchhoff)积分方程：
Figure BDA0001978583220000081
式中，
Figure BDA0001978583220000082
分别表示从坐标中心指向对应场点的矢量，Nl(ξ)表示局部直角坐标系下单元的第l号节点处的插值形函数；J(ξ)表示雅克比系数，用于主坐标系和局部坐标系之间的变换；
Figure BDA0001978583220000083
表示在m号单元l号节点上的声压值；
Figure BDA0001978583220000084
表示在m号单元l号节点上的法向振速。
Figure BDA0001978583220000085
代表格林函数，
Figure BDA0001978583220000086
表示弹性结构体表面外法向偏导数。
当S为光滑面时，α(p)取值由相应场点所在的位置决定：
Figure BDA0001978583220000087
若S是非光滑面，当场点位于S面上时取值为：
Figure BDA0001978583220000088
针对真实声源，其表面振速声压与场点中声压之间的格林函数应取自由场中格林函数，其法向导数具体表达式如下：
Figure BDA0001978583220000089
Figure BDA0001978583220000091
由于在半空间环境下界面的反射效应，存在着镜像虚源，对于虚拟声源而言，其表面振速声压与场点中声压之间应满足半空间格林函数：
Figure BDA0001978583220000092
其中，
Figure BDA0001978583220000093
是场点p与源面上点q的距离；
Figure BDA0001978583220000094
是场点p与虚源面上点q'的距离；Cri表示边界的反射系数，满足-1≤Cri≤1，具体取值如下：
Figure BDA0001978583220000095
其中，θ表示对应声波的入射角；β表示反射界面的比导纳。因此，设全息测量面上有i个测点，则对应测点的声压反射系数可表示为矩阵的形式Cr＝diag(Cr1,Cr2,···,Cri)。
步骤三：建立相应的传递矩阵；
由于每个单元与其相邻的单元有共同的节点，将节点进行统一的排序并相邻的重复节点进行合并，(1)式转换为：
Figure BDA0001978583220000096
其中：
Figure BDA0001978583220000097
假设在辐射声场中取Lo个测点、声源表面和声源内部分别取L、Li个考察点，全息面上各个测点的声压用
Figure BDA0001978583220000098
来表示，上式可分解成以下三个矩阵形式的方程：
Figure BDA0001978583220000099
[C]L×N[Ps]N×1＝[D]L×N[Vn]N×1 r∈S0 (11)
Figure BDA00019785832200000910
其中，{PH}表示Do区域内某场点的声压；[Ps]表示弹性结构体源面上的声压；[Vn]表示弹性结构体源面上的法向振速；[Pn]表示Di区域内场点的声压；相应的传递矩阵计算公式如下：
Figure BDA0001978583220000101
Figure BDA0001978583220000102
Figure BDA0001978583220000103
式中，当场点取节点i'处，i'与i重合时δ(i-i')取1；两者不重合时δ(i-i')取0。
Figure BDA0001978583220000104
Figure BDA0001978583220000105
Figure BDA0001978583220000106
利用以上计算方法，即可根据不同位置关系，建立不同的传递矩阵，具体包括：(1)建立实际声源表面振速与表面声压之间的传递矩阵、镜像虚源表面振速与实际声源表面声压之间的传递矩阵；(2)结合各全息测点与真实声源的位置关系，建立实际声源表面振速与全息声压之间的传递矩阵；(3)结合各全息测点与镜像虚源的位置关系，建立虚拟声源的表面振速与全息声压之间的传递矩阵。
步骤四：建立全息变换关系式
在半空间环境中，全息面上的测量结果为实际声源产生的辐射声压与界面的反射声的叠加结果，采用镜像法，将界面反射声等效为对称虚源的辐射噪声，从而实现半空间中的声源重建。根据镜像原理，虚源上任意节点位置分布情况都和实际声源完全一致，而且虚源源面上每个节点的表面法向振速和实际源面的表面法向振速对应相等；考虑水面对声源的影响，实际声源应看作理想声源和虚源对理想声源的作用两者叠加而成。
根据镜像法的原理建立半空间辐射声场理论模型，如图4所示。以水面为对称界面，将声源S关于对称界面进行对称，假设在对称位置处有一个与声源各处相位相反、幅值大小完全一致的“虚拟源”S*，虚源上的声压和振速分别为{PS*}和{VS*}，则考虑虚源影响下，真实的源面声压可表示为：
{P′S}＝{PS}+[WSS*]{VS*}＝[WS]{VS}+[WSS*]{VS*} (19)
式中，{PS}、{VS}分别表示为理想声源的表面声压和振速；[WS]表示理想声源表面声压与振速之间的传递矩阵。[WSS*]表示虚拟声源表面振速与理想声源表面声压之间的传递矩阵。而此时实际的辐射声场是虚拟声源和考虑虚源后的真实声源共同影响的叠加声场。因此，全息测量面实际声压为：
{Ph}＝{Phd}+{Phr} (20)
其中，{Phd}表示真实声源产生的直达声；{Phr}表示虚拟声源产生的反射声。
将虚源影响下的真实源面声压代入源面振速计算公式(9)可得：
{V′S}＝[WS]-1{P′S}＝{VS}+[WS]-1[WSS*]{VS*} (21)
式中，PS'、VS'分别表示虚源影响下真实源面的声压与振速。结合声压的反射系数，则全息面上实际测量的声压公式(20)转化为下式：
{Ph}＝[Wd]{VS'}+[Cr][Wr]{VS} (22)
式中，[Wd]代表真实源面在全息面上产生的辐射声压和真实源面振速的传递矩阵，[Wr]代表虚拟源面在全息面上产生的辐射声压和虚拟源面振速的传递矩阵。
将式(21)代入(22)中，上式即可转化为：
{Ph}＝[Wd]{VS'}+[Cr][Wr]{VS}＝{[Wd]+[Cr][Wr]+[Wd][WS]-1[WSS*]}{VS}＝{Wsh}{VS} (23)
其中{Wsh}＝{[Wd]+[Cr][Wr]+[Wd][WS]-1[WSS*]}即为半空间的源面到全息面的传递矩阵。
如上，即可实现半空间环境下的全息重建。将上式中的传递矩阵采用奇异分解法进行求逆计算，得到：
{VS}＝{Wsh}-1{Ph}＝[U][∑]-1[V]H{Ph} (24)
上式即为半空间边界元法声全息的声源振速重建公式。
步骤五：全息重建过程中的正则化处理
当求解病态矩阵过程中产生大量的离散值时，正则化处理较为适用就是Tikhonov正则化方法。通过原有方程中，施加残余范数和单边约束形成的加权组合取最小的约束条件，即{VS}reg须满足以下条件：
Figure BDA0001978583220000121
其中，λ取值大于零，代表正则化参数；(VS)*代表源面法向振速的初始估计值；L代表罚矩阵，一般为正定矩阵，相当于要求未知解满足一定的光滑性约束条件。
对上述方程进行求解，得到广义Tikhonov正则化的解，
{VS}reg＝(WHW+λ2LTL)-1WH{Ph} (26)
由上式可知，利用(WHW+λ2LTL)-1WH取代源面振速重建公式中的病态矩阵的逆W-1。
当L＝In，(VS)*＝0时，即标准正则化方法，其源面的法向振速可表示为：
Figure BDA0001978583220000122
为了减少测量误差的影响，正则化系数的取值尤为关键，以下给出标准Tikhonov正则化方法的滤波系数表达式：
Figure BDA0001978583220000123
其中正则化参数λ最为常见的选取方法包括广义交叉检验法和L曲线准则。广义交叉验证法主要通过下述方程来确定：
Figure BDA0001978583220000124
其中，WI代表正则化解{VS}reg对应的传递矩阵；trace()指传递矩阵的迹；可以看出，当方程求最小值时，对应的λ即为所求的正则化参数。
以L曲线为基础的L曲线准则(L-Curve Criterion)与广义交叉验证法相同，对误差范数是否已知不作要求。其中，L曲线描述的是利用所有可行的正则化参数λ计算
Figure BDA0001978583220000125
并以对数尺度绘制残余范数与解范数的关系。然后，通过对残余范数与解范数的对比选取有效的正则化参数。图5为L曲线示意图，图中拐角的左半边部分，由于λ太小，出现欠正则化现象，虽然采用正则化处理，但仍存在某些奇异值较小的分量导致未对误差实现有效的控制；而在L曲线拐角的右半边部分，由于λ太大，出现过正则化现象，滤波后的奇异值太少导致有效信息的丢失。因此，L曲线的拐角处才是能有效实现滤波的正则化参数，并且L曲线越是满足Picard条件，其拐角越突出。
当L≠In，且为可逆方阵时，需要将一般形式Tikhonov正则化化为标准形式。首先将{VS}、W逆进行转换，得到新的形式
Figure BDA0001978583220000131
并代入原来方程中，将式(25)转化为标准形式：
Figure BDA0001978583220000132
标准正则化下的源面法向振速通过将转换后的矩阵进行奇异值分解求得，并进行转换，即可得到非标准形式的源面法向振速
Figure BDA0001978583220000133
实施例：
对本发明的有益效果如下方式得以验证：
(1)在Matlab仿真条件下，对该方法进行仿真实验：
仿真参数如下：
模型几何参数：球壳，半径0.75m，令中心位于坐标原点，划分网格后，形成218个节点及216个四边形单元；
模型内部声源参数：位于(-0.1,-0.1,0.1)，点声源，幅值为1Pa；
流体介质参数：水，密度1000kg/m3,声速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表示弹性结构体表面，Di表示弹性结构体内部区域，弹性结构体外部区域的Do充满密度为ρ，声速为c的流体，弹性结构体外部的场点表示为p，弹性结构体上的点表示为q，将弹性结构体表面边界S进行离散化处理，分为M个小单元，共N个节点，依据有限元和边界元理论，建立基于边界元理论的亥姆霍兹-基尔霍夫积分方程：
Figure FDA0003636225310000011
式中，
Figure FDA0003636225310000012
分别表示从坐标中心指向对应场点的矢量，Nl(ξ)表示局部直角坐标系下单元的第l号节点处的插值形函数；J(ξ)表示雅克比系数，用于主坐标系和局部坐标系之间的变换；
Figure FDA0003636225310000013
表示在m号单元l号节点上的声压值；
Figure FDA0003636225310000014
表示在m号单元l号节点上的法向振速，
Figure FDA0003636225310000015
代表格林函数，
Figure FDA0003636225310000016
表示弹性结构体表面外法向偏导数，
当S为光滑面时，α(p)取值由相应场点所在的位置决定：
Figure FDA0003636225310000017
若S是非光滑面，当场点位于S面上时取值满足：
Figure FDA0003636225310000018
针对真实声源，表面振速声压与场点中声压之间的格林函数取自由场中格林函数，法向导数具体表达式满足：
Figure FDA0003636225310000021
Figure FDA0003636225310000022
虚拟声源表面振速声压与场点中声压之间满足半空间格林函数：
Figure FDA0003636225310000023
其中，
Figure FDA0003636225310000024
是场点p与源面上点q的距离；
Figure FDA0003636225310000025
是场点p与虚源面上点q'的距离；Cri表示边界的反射系数，满足-1≤Cri≤1，具体取值满足：
Figure FDA0003636225310000026
其中，θ表示对应声波的入射角；β表示反射界面的比导纳，设全息测量面上有i个测点，则对应测点的声压反射系数可表示为矩阵的形式Cr＝diag(Cr1,Cr2,…,Cri)；
步骤三：建立相应的传递矩阵，具体为：
将节点进行统一的排序并相邻的重复节点进行合并，步骤二所述
Figure FDA0003636225310000027
转换为：
Figure FDA0003636225310000028
其中：
Figure FDA0003636225310000029
假设在辐射声场中取Lo个测点、声源表面和声源内部分别取L、Li个考察点，全息面上各个测点的声压用
Figure FDA00036362253100000210
来表示，
Figure FDA00036362253100000211
上式可分解成以下三个矩阵形式的方程：
Figure FDA0003636225310000031
[C]L×N[Ps]N×1＝[D]L×N[Vn]N×1 r∈S0
Figure FDA0003636225310000032
其中，{PH}表示Do区域内某场点的声压；[Ps]表示弹性结构体源面上的声压；[Vn]表示弹性结构体源面上的法向振速；[Pn]表示Di区域内场点的声压；相应的传递矩阵计算公式满足：
[A]＝[a11,a12,…,a1N],
Figure FDA0003636225310000033
[B]＝[b11,b12,…,b1N],
Figure FDA0003636225310000034
[C]＝[c11,c12,…,c1N],
Figure FDA0003636225310000035
式中，当场点取节点i'处，i'与i重合时δ(i-i')取1；两者不重合时δ(i-i')取0；
[D]＝[d11,d12,…,d1N],
Figure FDA0003636225310000036
[E]＝[e11,e12,…,e1N],
Figure FDA0003636225310000037
[F]＝[f11,f12,…,f1N],
Figure FDA0003636225310000038
步骤四：建立全息变换关系式，具体为：
将声源S关于对称界面进行对称，假设在对称位置处有一个与声源各处相位相反、幅值大小完全一致的“虚拟源”S*，虚源上的声压和振速分别为{PS*}和{VS*}，则考虑虚源影响下，真实的源面声压可表示为：
Figure FDA0003636225310000039
式中，{PS}、{VS}分别表示为理想声源的表面声压和振速；[WS]表示理想声源表面声压与振速之间的传递矩阵，
Figure FDA00036362253100000310
表示虚拟声源表面振速与理想声源表面声压之间的传递矩阵，实际的辐射声场是虚拟声源和考虑虚源后的真实声源共同影响的叠加声场，则全息测量面实际声压为：
{Ph}＝{Phd}+{Phr}
其中，{Phd}表示真实声源产生的直达声；{Phr}表示虚拟声源产生的反射声；
将虚源影响下的真实源面声压代入源面振速计算公式
Figure FDA0003636225310000041
可得：
Figure FDA0003636225310000042
式中，PS'、VS'分别表示虚源影响下真实源面的声压与振速，结合声压的反射系数，则全息面上实际测量的声压公式{Ph}＝{Phd}+{Phr}转化为下式：
{Ph}＝[Wd]{VS'}+[Cr][Wr]{VS}
式中，[Wd]代表真实源面在全息面上产生的辐射声压和真实源面振速的传递矩阵，[Wr]代表虚拟源面在全息面上产生的辐射声压和虚拟源面振速的传递矩阵；
将
Figure FDA0003636225310000043
代入{Ph}＝[Wd]{VS'}+[Cr][Wr]{VS}中，{Ph}＝[Wd]{VS'}+[Cr][Wr]{VS}可转化为：
Figure FDA0003636225310000044
其中
Figure FDA0003636225310000045
即为半空间的源面到全息面的传递矩阵；
将上式中的传递矩阵采用奇异分解法进行求逆计算，得到半空间边界元法声全息的声源振速重建公式：
{VS}＝{Wsh}-1{Ph}＝[U][∑]-1[V]H{Ph}
步骤五：全息重建过程中的正则化处理，具体为：
采用Tikhonov正则化方法，通过施加残余范数和单边约束形成的加权组合对{VS}＝{Wsh}-1{Ph}＝[U][∑]-1[V]H{Ph}取最小的约束条件，即{VS}reg满足以下条件：
Figure FDA0003636225310000046
其中，λ取值大于零，代表正则化参数；(VS)*代表源面法向振速的初始估计值；L代表罚矩阵；
对
Figure FDA0003636225310000051
进行求解，得到广义Tikhonov正则化的解，满足：
{VS}reg＝(WHW+λ2LTL)-1WH{Ph}
由上式可知，利用(WHW+λ2LTL)-1WH取代源面振速重建公式中的病态矩阵的逆W-1；
当L＝In，(VS)*＝0时，即标准正则化方法，源面的法向振速可表示为：
Figure FDA0003636225310000052
标准Tikhonov正则化方法的滤波系数表达式满足：
Figure FDA0003636225310000053
正则化参数λ的选取方法为L曲线准则，L曲线的拐角处值为正则化参数；
当L≠In，且为可逆方阵时，将一般形式Tikhonov正则化化为标准形式，首先将{VS}、W逆进行转换，得到新的形式
Figure FDA0003636225310000054
将
Figure FDA0003636225310000055
转化为标准形式：
Figure FDA0003636225310000056
标准正则化下的源面法向振速通过将转换后的矩阵进行奇异值分解求得，并进行转换，即可得到非标准形式的源面法向振速
Figure FDA0003636225310000057
2.根据权利要求1所述的一种半空间环境下的边界元法近场声全息变换方法，其特征在于：步骤五所述L为正定矩阵。
3.根据权利要求1所述的一种半空间环境下的边界元法近场声全息变换方法，其特征在于：步骤五所述正则化参数λ的选取方法为广义交叉检验法，具体为：
正则化参数λ通过方程
Figure FDA0003636225310000058
确定，其中，WI代表正则化解{VS}reg对应的传递矩阵；trace()指传递矩阵的迹；当方程求最小值时，对应的λ即为所求的正则化参数。
CN201910141127.8A 2019-02-26 2019-02-26 一种半空间环境下的边界元法近场声全息变换方法 Active CN109885945B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201910141127.8A CN109885945B (zh) 2019-02-26 2019-02-26 一种半空间环境下的边界元法近场声全息变换方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201910141127.8A CN109885945B (zh) 2019-02-26 2019-02-26 一种半空间环境下的边界元法近场声全息变换方法

Publications (2)

Publication Number Publication Date
CN109885945A CN109885945A (zh) 2019-06-14
CN109885945B true CN109885945B (zh) 2022-08-02

Family

ID=66929426

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201910141127.8A Active CN109885945B (zh) 2019-02-26 2019-02-26 一种半空间环境下的边界元法近场声全息变换方法

Country Status (1)

Country Link
CN (1) CN109885945B (zh)

Families Citing this family (13)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110401898B (zh) * 2019-07-18 2021-05-07 广州酷狗计算机科技有限公司 输出音频数据的方法、装置、设备和存储介质
CN110554358B (zh) * 2019-09-25 2022-12-13 哈尔滨工程大学 一种基于虚拟球阵列扩展技术的噪声源定位识别方法
CN110763328B (zh) * 2019-11-18 2021-10-22 湖北文理学院 一种半空间声场重建方法和装置
JP7374760B2 (ja) * 2019-12-26 2023-11-07 Ｔｏｙｏ Ｔｉｒｅ株式会社 音源探査方法
CN111060596B (zh) * 2020-01-02 2022-04-05 合肥工业大学 一种吸声材料吸声系数的测量方法
CN111812587B (zh) * 2020-07-06 2023-04-07 上海交通大学 基于机器视觉和全息方法的声场测试分析方法及系统
CN112577592B (zh) * 2020-11-27 2022-10-28 哈尔滨工程大学 基于空间傅里叶变换的有限空间平面近场声全息测量方法
CN112926231B (zh) * 2020-11-27 2023-06-20 哈尔滨工程大学 一种基于等效源法的有限空间中近场声全息测量方法
CN113536554B (zh) * 2021-07-01 2022-08-02 哈尔滨工程大学 一种采用压缩等效源法的封闭空间内声场预测方法
CN113639934B (zh) * 2021-08-31 2024-03-26 郑州大学 气体泄漏位置三维定位方法、定位系统
CN114630260B (zh) * 2022-02-17 2023-11-10 安徽建筑大学 一种半空间均匀流中声场重建方法
CN114543979B (zh) * 2022-02-17 2024-05-03 浙江工业大学 有界空间中基于近场声全息的声源直接辐射远场声学量的预测方法
CN117436293A (zh) * 2023-12-21 2024-01-23 国网浙江省电力有限公司电力科学研究院 基于声场重构的低频变压器测点仿真方法和电子设备

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017107337A1 (zh) * 2015-12-22 2017-06-29 合肥工业大学 一种基于改进的按位替换法求矩阵三角分解的模块及方法
CN107576388A (zh) * 2017-08-22 2018-01-12 哈尔滨工程大学 一种浅海信道下三维结构声源辐射声场预报方法

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2017107337A1 (zh) * 2015-12-22 2017-06-29 合肥工业大学 一种基于改进的按位替换法求矩阵三角分解的模块及方法
CN107576388A (zh) * 2017-08-22 2018-01-12 哈尔滨工程大学 一种浅海信道下三维结构声源辐射声场预报方法

Non-Patent Citations (4)

* Cited by examiner, † Cited by third party
Title
Nearfield acoustic holography based on inverse boundary element method for moving sound source identification;Chen Mengying 等;《Acta Acustica》;20110930;第36卷(第5期);第489-584页 *
半自由场波叠加噪声源识别方法研究;潘汉军等;《中国机械工程》;20060415(第07期);第733-737页 *
基于矢量阵测量的振速重构近场声全息技术;胡博等;《机械工程学报》;20101112;第46卷(第22期);第16-23页 *
运动声源的边界元声全息识别方法研究;陈梦英等;《声学学报》;20110915;第36卷(第05期);第489-495页 *

Also Published As

Publication number Publication date
CN109885945A (zh) 2019-06-14

Similar Documents

Publication Publication Date Title
CN109885945B (zh) 一种半空间环境下的边界元法近场声全息变换方法
Fernandez-Grande et al. A sparse equivalent source method for near-field acoustic holography
CN107885934B (zh) 基于耦合fem-pe的海洋信道下弹性结构声辐射预报方法
CN109883532B (zh) 一种声源识别与声场预报方法
CN105866740B (zh) 一种基于压缩感知的水声匹配场定位方法
CN112577592B (zh) 基于空间傅里叶变换的有限空间平面近场声全息测量方法
CN113063490A (zh) 一种基于声压和质点振速双面测量的声场分离方法
CN112926231B (zh) 一种基于等效源法的有限空间中近场声全息测量方法
Wu et al. A boundary element method based near field acoustic holography in noisy environments
CN114779170A (zh) 一种浅海近场声源定位方法
Lu One-way large range step methods for Helmholtz waveguides
de Hoop et al. On the construction of virtual interior point source travel time distances from the hyperbolic Neumann-to-Dirichlet map
Liang et al. Scattering pattern analysis and true‐amplitude generalized Radon transform migration for acoustic transversely isotropic media with a vertical axis of symmetry
Ladopoulos Non-linear Seismic Wave Motion in Elastodynamics with Application to Real-Time Expert Seismology
CN109164488B (zh) 一种梯形网格有限差分地震波场模拟方法
Wiskin et al. Inverse scattering from arbitrary two-dimensional objects in stratified environments via a Green’s operator
Nemitz et al. Topological sensitivity and FMM-accelerated BEM applied to 3D acoustic inverse scattering
Nadimi et al. An efficient acoustic scattering model based on target surface statistical descriptors for synthetic aperture sonar systems
Fan et al. Prediction of object geometry from acoustic scattering using convolutional neural networks
Wu et al. A novel hybrid superposition method for predicting ship seismic wave in shallow sea
Quijano et al. Geoacoustic inversion for the seabed transition layer using a Bernstein polynomial model
KR20120138702A (ko) 해저면 지형을 반영한 지하구조 영상화 방법
Li et al. Fast frequency-domain forward and inverse methods for acoustic scattering from inhomogeneous objects in layered media
CN118857440A (zh) 一种预报水下结构跨介质声地场的波数积分叠加计算方法
Dai et al. A novel finite-infinite coupling approach for 3D simulation of borehole-to-surface electrical potential

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