CN112697269A - 一种基于模拟声场先验信息的声强稀疏测量高分辨率成像方法 - Google Patents

Abstract

本发明涉及一种基于模拟声场先验信息的声强稀疏测量高分辨率成像方法，首先，对目标声源元件进行结构噪声辐射分析，获取其模拟声强分布的可视化图像；随后提取模拟声场的先验信息包括：声强强点的位置、稀疏矩阵以及稀疏度K；再结合声强强点位置和混合高斯分布模型，完成声强测点的位置设计，并根据设计好的测点位置推导出观测矩阵的数学模型；最后基于稀疏矩阵、观测矩阵和贪婪算法所构成的压缩感知数据处理框架，对声强稀疏测量值进行算法重构，实现高分辨率成像。

Description

一种基于模拟声场先验信息的声强稀疏测量高分辨率成像 方法

技术领域

本发明涉及声强测量成像的技术领域，特别是一种基于模拟声场先验信息的声强稀疏测量高分辨率成像方法。

背景技术

声强测量成像技术以图像形式直观展示振动体外部的声强场能分布，主要应用于车辆动力系统和冰箱制冷系统等系统级产品的噪声源定位，特别地，包含声源细节信息的图像纹理特征还可用于系统的故障诊断，该方法便捷且辨识度高，但其面向小尺寸元件成像时还存在分辨率不足的问题。

声强成像的分辨率受限源于两个方面的问题：声强测量固有的系统误差和瑞利判据的影响。声强测量固有的系统误差包括：有限差分误差、相位匹配误差和近场误差。其中有限差分误差与相位匹配误差分别影响了声强测量的高频及低频精度，进一步限制了双声传感器的间距；近场误差认为面对多源相干声场，声强探头应与声源面的距离应大于3倍的双声传感器间距；瑞利判据则对声波信号的测试过程中声传感器应尽可能靠近测试对象，以捕捉近场消逝波成分，这种消逝波包含辐射声场的细节信息。由上述分析可知，声强成像的分辨率不足体现在两个方面：近场误差与瑞利判据的矛盾以及声强频带宽度的不足。

针对上述问题，可分别从两个方面入手解决：提高声场采集信号的信噪比和突破奈氏采样的束缚以拓展信号采样的频宽。压缩感知理论指出只要信号的稀疏或可压缩的，就能通过少量观测数据高概率恢复原始信号。压缩感知不仅能够有效提高采集信号的带宽，同时这种亚采样髙恢复的方法能够提高同种采样频率下的信噪比，从而使得重构信号突破瑞利分辨率的限制。

发明内容

有鉴于此，本发明将压缩感知理论与声强成像技术相结合提供一种基于模拟声场先验信息的声强稀疏测量高分辨率成像方法，解决原有声强成像技术无法满足振源密集型元件声源辨识精度问题。

本发明采用以下方案实现：一种基于模拟声场先验信息的声强稀疏测量高分辨率成像方法，包括以下步骤：

步骤S1：对目标声源元件进行结构噪声辐射分析，获取目标元件的模拟声场图像；

步骤S2：对步骤S1中获取的模拟声场图像进行稀疏表达，确定能够还原完整声场的稀疏表达矩阵及稀疏度；

步骤S3：根据模拟声场图像设计声强测点分布，并推导出对应观测矩阵数学式；

步骤S4：基于贪婪算法，分别输入步骤S2中确定的稀疏表达矩阵及稀疏度值，步骤S3中确定的观测矩阵，输出目标元件模拟声场的稀疏系数，进一步还原出模拟声场的可视化图像，用以实现目标元件声场稀疏测点信息的完整声场重构，确立压缩感知的算法框架；

步骤S5：基于步骤S4中所建立的压缩感知算法框架下，根据混合高斯分布的测点生成方法，进行随机测点的位置布置，并对各测点进行声强离散测量；采用压缩感知算法程序对稀疏离散测点数据进行后处理重构，最终实现声强稀疏测量的高分辨率成像。

进一步地，步骤S1中所述获取目标元件的模拟声场图像通过三个部分，分别为：目标元件内部激振信息的数值模拟、振动激励在结构体内部至壳体的传递和外表面质点振动引起的空气噪声辐射。

进一步地，所述步骤S3的具体内容为：

基于混合高斯分布模型设计声强测点分布：首先根据模拟声场图像中模拟声场目标区域的声强分布，提取声强高值区域中心点的位置坐标作为程序初始值；随后结合声强区域的大小，并设置与传统声强测量相同的测点数量；最后调用高斯混合模型生成随机测点，并完成测点绘制；

由声强稀疏测点位置向观测矩阵数学式的推导过程如下：

假设模拟声强的完整分布由21×21即441个数据点构成，将其二维声强分布数据用矩阵的形式进行表示如下，记为A；

将矩阵A取每一列元素排成一维列向量为x＝[x0,x1,x2,x3,…,x437,x438,x439,x440]T；混合高斯分布的随机测点位置为100个，根据这些测点信息在矩阵A中的位置，同样取每一列逐个进行排列，假设随机测点向量y＝[x1,x2,x9,…,x433,x440]T；记观测矩阵为B＝(b1,b2,b3,…,b99,b100)T，则y＝Φx,即：x1＝Φ1x,x2＝Φ2x，x9＝Φ3x，…，x433＝Φ99x,x440＝Φ100x，由此推导出观测矩阵如下：

与现有技术相比，本发明具有以下有益效果：

本发明大地减少了声强测量的工作量，有效提高了信噪比，降低对环境噪声的要求，同时拓展了声强测量频宽，其分辨率的提升将有助于对振源密集型元件级产品的声源辨识。

附图说明

图1为本发明实施例的方法流程图。

图2为本发明实施例的DCT冗余字典构造的程序流程图。

图3为本发明实施例的模拟声场与不同稀疏度的压缩还原声场，其中，图3(a)为模拟声场，图3(b)为稀疏度35的压缩还原声场图，图3(c)为稀疏度45的压缩还原声场图，图3(d)为稀疏度55的压缩还原声场图。

图4为本发明实施例的基于混合高斯分布模型的声强测点设计流程图。

图5为本发明实施例的贪婪算法的流程图。

图6为本发明实施例的基于模拟声场先验信息的声强稀疏测量的实验方案。

图7为本发明实施例的应用于斜轴式轴向柱塞马达声强成像的结果示例图。

具体实施方式

下面结合附图及实施例对本发明做进一步说明。

应该指出，以下详细说明都是例示性的，旨在对本申请提供进一步的说明。除非另有指明，本文使用的所有技术和科学术语具有与本申请所属技术领域的普通技术人员通常理解的相同含义。

需要注意的是，这里所使用的术语仅是为了描述具体实施方式，而非意图限制根据本申请的示例性实施方式。如在这里所使用的，除非上下文另外明确指出，否则单数形式也意图包括复数形式，此外，还应当理解的是，当在本说明书中使用术语“包含”和/或“包括”时，其指明存在特征、步骤、操作、器件、组件和/或它们的组合。

本实施例提供一种基于模拟声场先验信息的声强稀疏测量高分辨率成像方法，包括以下步骤：

步骤S1：对目标声源元件进行结构噪声辐射分析，获取目标元件的模拟声场图像；

步骤S2：对步骤S1中获取的模拟声场图像进行稀疏表达，确定能够还原完整声场的稀疏表达矩阵及稀疏度；

步骤S3：根据模拟声场图像设计声强测点分布，并推导出对应观测矩阵数学式；

步骤S4：基于贪婪算法，分别输入步骤S2中确定的稀疏表达矩阵及稀疏度值，步骤S3中确定的观测矩阵，输出目标元件模拟声场的稀疏系数，进一步还原出模拟声场的可视化图像，用以实现目标元件声场稀疏测点信息的完整声场重构，确立压缩感知的算法框架；

步骤S5：基于步骤S4中所建立的压缩感知算法框架下，根据混合高斯分布的测点生成方法，进行随机测点的位置布置，并对各测点进行声强离散测量；采用压缩感知算法程序对稀疏离散测点数据进行后处理重构，最终实现声强稀疏测量的高分辨率成像。

在本实施例中，步骤S1中所述获取目标元件的模拟声场图像通过三个部分，分别为：目标元件内部激振信息的数值模拟、振动激励在结构体内部至壳体的传递和外表面质点振动引起的空气噪声辐射。

在本实施例中，所述步骤2中的稀疏表达的过程在于构造能够对模拟声场进行稀疏分解的字典矩阵ψ，稀疏表达的结果则是将数据量巨大的模拟声场x用稀疏系数α来表示，以减少有效信息的获取量。

结合DCT良好的去相关性和冗余字典的易于稀疏特性，构造基于DCT的冗余字典矩阵，DCT冗余字典的构造具体分为以下三个步骤：

步骤SA：根据待稀疏表达信号的维数，确定DCT冗余字典的行数，令模拟声强的完整分布由441(21×21)个数据点构成，由此确定其视图对应的DCT字典行数为21；

步骤SB：对字典矩阵进行分频采样，DCT变换只针对奇数点进行计算，而其冗余字典变换还对偶数点进行计算，因此确定分频采样后的列数为42；

步骤SC：对分频采样后的矩阵进行调制，完成冗余字典构造；具体的程序内容如图2所示。在完成DCT冗余字典的设计后，即可设定几组稀疏度，进行稀疏表达后再还原各组图像，对比各组图像的还原结果如图3所示，选择能完整还原图像强点并尽可能小的45作为稀疏度值。

在本实施例中，所述步骤S3的具体内容为：

基于混合高斯分布模型设计声强测点分布：首先根据模拟声场图像中模拟声场目标区域的声强分布，提取声强高值区域中心点的位置坐标作为程序初始值；随后结合声强区域的大小，并设置与传统声强测量相同的测点数量；最后调用高斯混合模型生成随机测点，并完成测点绘制；

由声强稀疏测点位置向观测矩阵数学式的推导过程如下：

假设模拟声强的完整分布由21×21即441个数据点构成，将其二维声强分布数据用矩阵的形式进行表示如下，记为A；

将矩阵A取每一列元素排成一维列向量为x＝[x0,x1,x2,x3,…,x437,x438,x439,x440]T；混合高斯分布的随机测点位置为100个，根据这些测点信息在矩阵A中的位置，同样取每一列逐个进行排列，假设随机测点向量y＝[x1,x2,x9,…,x433,x440]T；记观测矩阵为B＝(b1,b2,b3,…,b99,b100)T，则y＝Φx,即：x1＝Φ1x,x2＝Φ2x，x9＝Φ3x，…，x433＝Φ99x,x440＝Φ100x，由此推导出观测矩阵如下：

在本实施例中，所述步骤4中的信号重构算法本例选用贪婪算法，指求解凸优化问题的一列算法，用于低维稀疏信息向高维模拟声场的欠定映射问题。

较佳的，如图1所示，本实施例具体的可分为三个步骤执行：

步骤1，以结构辐射声场分析方法获取目标元件模拟声强分布的可视化图像，为压缩感知应用提供先验信息。通过建立元件的动力学联合仿真模型，使用ADMAS仿真平台和AMESim软件联合仿真获得内部激振信息；以此振源信息，结合有限元动力学响应分析方法，将激振信息导入ansys软件对外表面结构体进行振动响应分析，获取马达外表面振动位移特性；将外表面振动数据作为边界条件导入LMS.Virtual软件，采用边界元方法完成整体结构声辐射表征，并获得元件的模拟声强分布云图。

步骤2，研究压缩感知的算法原理，理论上建立声强的稀疏测量成像方法。构造DCT(Discrete Cosine Transform)冗余字典的稀疏矩阵，对斜轴式马达进行模拟声强分布云图的稀疏表达，确定稀疏度(能够还原完整声场的稀疏信息个数)；根据马达模拟声场的声强分布特点，结合混合高斯分布模型，指导设计随机测点的位置，并推导出对应观测矩阵数学式；基于贪婪重构算法，对随机测点信息进行完整的声场重构。

步骤2中，稀疏表达的过程在于构造能够对模拟声场进行稀疏分解的字典矩阵ψ，稀疏表达的结果则是将数据量巨大的模拟声场x用稀疏系数α来表示，以减少有效信息的获取量。本例旨在结合DCT良好的去相关性和冗余字典的易于稀疏特性，构造基于DCT的冗余字典矩阵，DCT冗余字典的构造具体分为以下三个步骤：根据待稀疏表达信号的维数，确定DCT冗余字典的行数，假设模拟声强的完整分布由441(21×21)个数据点构成，由此确定其视图对应的DCT字典行数为21；对字典矩阵进行分频采样，DCT变换只针对奇数点进行计算，而其冗余字典变换还对偶数点进行计算，因此确定分频采样后的列数为42；进一步对分频采样后的矩阵进行调制，完成冗余字典构造。具体的程序内容如图2所示。

DCT冗余字典的构造具体分为以下三个步骤：

S1:根据待稀疏表达信号的维数，确定DCT字典的行数。前述章节中确定的斜轴式柱塞马达结构噪声辐射计算的上视与侧视方向声强分布的数值矩阵维数为M×M，由此分别确定其视图对用的DCT字典行数分别为M。

S2:对字典矩阵进行分频采样。DCT变换只针对奇数点进行计算，而其冗余字典变换还对偶数点进行计算，因此确定分频采样后的列数为2M。

S3:进一步对分频采样后的矩阵进行调制，完成冗余字典构造。

在完成DCT冗余字典的设计后，即可设定几组稀疏度，进行稀疏表达后再还原各组图像，对比各组图像的还原结果如图3所示，选择能完整还原图像强点并尽可能小的45作为稀疏度值。

图4展示本实施例步骤2中基于混合高斯分布模型的声强测点设计流程，首先根据模拟声场目标区域的声强分布，提取声强高值区域中心点的坐标作为程序初始值；随后结合声强区域的大小，并设置与传统声强测量相同的测点数量；最后调用高斯混合模型生成随机测点，并完成测点绘制。

如图4所示，混合高斯分布的声强测点设计的具体步骤如下：

S1:对斜轴式马达模拟声场特点进行区域划分，具体的区分方式是围绕区域声强的强弱，声强大小的波动以及理论声源位置进行区分。

S2:提取各区域的核心坐标作为(u_x,u_y)，同时按照区域大小设计协方差值δx与δy。

S3:调用混合高斯模型，产生随机的分散点坐标。

步骤2中由声强稀疏测点位置向观测矩阵数学式的推导过程如下：

假设模拟声强的完整分布由441(21×21)个数据点构成，将其二维声强分布数据用矩阵的形式进行表示如下，记为A。

将矩阵A取每一列元素排成一维列向量为x＝[x0,x1,x2,x3,…,x437,x438,x439,x440]T。

混合高斯分布的随机测点位置为100个，根据这些测点信息在矩阵A中的位置，同样取每一列逐个进行排列，假设随机测点向量y＝[x1,x2,x9,…,x433,x440]T。记观测矩阵为B＝(b1,b2,b3,…,b99,b100)T，则y＝Φx,即：x1＝Φ1x,x2＝Φ2x，x9＝Φ3x，…，x433＝Φ99x,x440＝Φ100x，由此推导出观测矩阵如下：

图5展示了步骤2中贪婪算法的流程。其算法步骤简述如下：

输入：①100×1764的传感矩阵A＝Φψ，其中Φ为基于声强测点位置的观测矩阵，ψ为DCT冗余字典矩阵；②声强测点位置下的观测向量y；③稀疏度值45。

输出：①模拟声场信息的稀疏系数theta。②残差residual＝y-A*theta。

初始化残差e0＝y，y代表观测信号；

选择与e0内积绝对值最大的原子，表示为φ1；

进一步建立矩阵Φt＝[φ1,φ2,…],定义其列空间的Hermitian矩阵为：

从e0减去其正交投影量产生新的残差:

e₁＝e₀-Pe₀＝(I-P)e₀

重新进行上述步骤的残差迭代:

e_m+1＝e_m-Pe_m＝(I-P)e_m

步骤3，基于前序建立的压缩感知声强场可视化理论，进行斜轴式柱塞马达局部声场的实验辨识与对比分析。根据步骤2所确定的随机测点位置分布，进行随机测点的位置布置，并对各测点进行声强离散测量；采用压缩感知算法程序对稀疏离散测点数据进行后处理重构，最终实现声强稀疏测量的高分辨率成像。

图6展示了步骤3中基于模拟声场先验信息的声强稀疏测量的实验方案。首先根据目标区域声场的模拟声强云图，结合混合高斯模型，对随机测点进行布置，再进一步完成对应测点的声强测量，由此获取马达目标区域声场的稀疏测点信息y；随后由测点位置推导出相应的观测矩阵Φ，与DCT冗余字典确定的稀疏矩阵ψ组成压缩感知传感矩阵A，并结合稀疏度值K和OMP算法程序CS_OMP(y,A,K)，构造随机测点数据的后处理程序；最后将稀疏测点信息带入数据后处理程序，完成斜轴式马达声场的压缩感知实验重构。

如图6所示，压缩感知实验方案流程如下：

S1:根据目标区域声场的模拟声强云图，结合混合高斯模型，对随机测点进行布置(布置方法为图4)，进一步完成随机测点的声强测量，获取马达声场的稀疏测点信息y.

S2:根据S1确定的测点位置推导出相应的观测矩阵Φ数学式，同时由DCT冗余字典(构造过程为图2)分别确定稀疏矩阵ψ和稀疏度值K，至此完成压缩感知传感矩阵A的构造，其中A＝Φψ.

S3:根据S2确定的传感矩阵A和稀疏度指K，结合OMP算法程序CS_OMP(y,A,K)(程序构造为图5)，组成了随机测点数据的后处理程序。

S4:最后将稀疏测点的声强数据，输入步骤S3确定的压缩感知后处理程序，进而构造出完整声场。

以斜轴式轴向柱塞马达为目标元件，进行基于模拟声场先验信息的声强稀疏成像的案例操作。图7展示了应用于斜轴式马达的传统声强成像以及基于模拟声场信息的声强稀疏成像对比，原图为彩色图，仅做去色处理，数字标识区域皆为深色区其结果显示：上视方向上，传统声强测量的可视化云图辨识出高值区域①，根据网格个数可以看出其半径为70mm，而压缩感知声强测量的可视化云图从中辨识出高值区域②和③，且两处中心的距离为30mm，并且能够明显看两处区域的声强强度。显然压缩感知算法对于声强测量的可视化分辨率具有明显的提升。侧视方向上，传统声强测量的可视化云图识别出的高值区域①，而压缩感知声强测量的可视化云图不仅辨识出高值区域②，还在靠近壳体下部位置辨识出区域③和④，这主要是由于压缩感知理论能够提高信号的信噪比，从而产生更多的声场细节信息。

以上所述仅为本发明的较佳实施例，凡依本发明申请专利范围所做的均等变化与修饰，皆应属本发明的涵盖范围。

Claims (3)

1.一种基于模拟声场先验信息的声强稀疏测量高分辨率成像方法，其特征在于：包括以下步骤：

步骤S1：对目标声源元件进行结构噪声辐射分析，获取目标元件的模拟声场图像；

步骤S2：对步骤S1中获取的模拟声场图像进行稀疏表达，确定能够还原完整声场的稀疏表达矩阵及稀疏度；

步骤S3：根据模拟声场图像设计声强测点分布，并推导出对应观测矩阵数学式；

步骤S4：基于贪婪算法，分别输入步骤S2中确定的稀疏表达矩阵及稀疏度值，步骤S3中确定的观测矩阵，输出目标元件模拟声场的稀疏系数，进一步还原出模拟声场的可视化图像，用以实现目标元件声场稀疏测点信息的完整声场重构，确立压缩感知的算法框架；所述压缩感知的算法框架即为步骤S2、步骤S3确立的稀疏表达矩阵、稀疏度、观测矩阵数学式以及所述完整声场重构的过程；

步骤S5：基于步骤S4中所建立的压缩感知算法框架下，根据混合高斯分布的测点生成方法，进行随机测点的位置布置，并对各测点进行声强离散测量；采用压缩感知算法程序对稀疏离散测点数据进行后处理重构，最终实现声强稀疏测量的高分辨率成像。

2.根据权利要求1所述的一种基于模拟声场先验信息的声强稀疏测量高分辨率成像方法，其特征在于：步骤S1中所述获取目标元件的模拟声场图像通过三个部分，分别为：目标元件内部激振信息的数值模拟、振动激励在结构体内部至壳体的传递和外表面质点振动引起的空气噪声辐射。

3.根据权利要求1所述的一种基于模拟声场先验信息的声强稀疏测量高分辨率成像方法，其特征在于：所述步骤S3的具体内容为：

基于混合高斯分布模型设计声强测点分布：首先根据模拟声场图像中模拟声场目标区域的声强分布，提取声强高值区域中心点的位置坐标作为程序初始值；随后结合声强区域的大小，并设置与传统声强测量相同的测点数量；最后调用高斯混合模型生成随机测点，并完成测点绘制；

由声强稀疏测点位置向观测矩阵数学式的推导过程如下：

假设模拟声强的完整分布由21×21即441个数据点构成，将其二维声强分布数据用矩阵的形式进行表示如下，记为A；

将矩阵A取每一列元素排成一维列向量为x＝[x0,x1,x2,x3,…,x437,x438,x439,x440]T；混合高斯分布的随机测点位置为100个，根据这些测点信息在矩阵A中的位置，同样取每一列逐个进行排列，假设随机测点向量y＝[x1,x2,x9,…,x433,x440]T；记观测矩阵为B＝(b1,b2,b3,…,b99,b100)T，则y＝Φx,即：x1＝Φ1x,x2＝Φ2x，x9＝Φ3x，…，x433＝Φ99x,x440＝Φ100x，由此推导出观测矩阵如下：

Images