Synchronization of non-synchronous measurement data based on Wiener predictor
-
摘要:
基于最小均方差准则, 推导了用于带参考传声器的非同步测量的维纳预测器表达式, 提出了基于维纳预测器的带有参考传声器的传声器阵列非同步测量数据同步化处理方法。维纳预测器通过传递参考传声器与阵列之间的时空相关性, 实现非同步测量数据同步化。分别通过仿真和实验将非同步测量结果和同步测量结果进行比较, 验证了所提基于维纳预测器的非同步测量数据同步化方法的准确性和可靠性, 同时也表明非同步测量方法可以提高阵列孔径和阵元密度。
Abstract:A data synchronization method is proposed to synchronize the data in microphone array non-synchronous measurements by employing the Wiener predictor derived from the minimum mean square error criterion. The Wiener predictor synchronizes the asynchronous measurement signal by harnessing the spatio-temporal correlation between the reference microphones and the array. To validate the proposed data synchronization method, extensive simulations and experiments are conducted to compare the results of non-synchronous and synchronous measurements. Both the simulation results and the experimental results confirm that the proposed method using the Wiener predictor for non-synchronous measurements is robust and reliable. Meanwhile, the results show the non-synchronous measurement method has the advantage of improving array aperture and element density.
-
引言
近场声全息是一种基于传声器阵列测量声场进行重建的经典方法, 可以获得包括声源位置[1]、声源表面声压分布[2]、振速分布[3]等信息, 也可以分析各声源的贡献量[4], 已经在多个领域得到大量应用。常用的近场声全息技术主要有基于空间傅里叶变换的近场声全息技术[5]、基于分布源边界点法的近场声全息技术[6]、基于等效源法的近场声全息技术[7]等, 其中基于等效源法的近场声全息可以避免窗效应、卷绕误差对声场重建、预测带来的影响, 也可以不受声源形状的限制, 同时模型简单、计算简便[8], 因此被广泛使用。
在近场声全息中, 要求测量孔径无限大, 以保证重建质量[9]。但是在实际应用中, 有限的测量孔径只能覆盖部分辐射区域, 在基于等效源法的近场声全息中, 有限的传声器数如果小于等效源数, 传递矩阵不定, 在高频段会产生较大误差[10], 因此随着声源频率的提升, 需要更多的节点以保证重建精度。针对测量孔径尺寸限制、阵元间距过大以及阵元数量不足等问题, 可以用等效源法实现全息面声压的外推计算, 建立基于等效源法的Patch近场声全息技术, 通过等效源线性叠加产生的声场来拟合实际声场, 实现全息面声压外推[11-12], 也可以采用插值[13]、人工智能算法[14-15]、稀疏采样、压缩感知方法[16-17]、压缩模态等效源方法[18]以及块稀疏贝叶斯[2, 19]等方法从算法层面进行改善。但是算法层面存在计算效率低[12]、算法本身误差会影响精度[11]、模型难以泛化[20-21]等不足, 采用稀疏测量方法易因采样量不足而在高频段导致过滤波, 难以满足重建精度要求[17]。
针对阵列孔径过小、阵元间距过大、传声器数量不足以及算法层面存在的问题, 也可以采用非同步测量方式, 即通过多次移动阵列面进行非同步测量, 并进行数据同步化处理, 实现更大测量孔径、更大阵元密度、更高空间采样率的测量效果。然而, 在非同步测量中, 阵列面多次移动会导致相邻两次测量之间的相位关系丢失, 从而影响声场重建、声场预测的实现。非同步测量可以分为不带参考传声器的非同步测量以及带参考传声器的非同步测量两种方式。采用不带参考传声器的非同步测量方法, 可以避免参考传声器数量、位置等对准确采集数据的影响, 但阵列移动有重叠率要求, 基于低秩假设和连续性假设, 通过算法实现相位信息补全[22-23], 或者基于平面波假设进行互质位置非同步测量, 通过算法实现相位同步[24]。尽管该方法目前在近场声全息有一定应用[25], 但其测量实施效率较低。相比之下, 带参考传声器的非同步测量方法不要求阵列面重叠, 可以显著减少测量次数。当前该方法使用广泛, 其进行相位补全的方法有条件谱分析法[26]、扫描互谱法[27-28], 但条件谱分析法将参考相位均置零, 不能将非同步测量数据同步至任一测量时刻[28]。针对带参考传声器的非同步测量方法, 文献[29]基于互谱矩阵的秩可等于子矩阵的秩的前提, 提出了一种计算互谱矩阵的方法。本文从最小均方差准则出发, 基于维纳预测器对该方法给出一种新的解释。
维纳预测器是一种常用的线性预测器, 针对带参考传声器的非同步测量数据同步化处理, 本文基于维纳预测器对非同步测量数据同步化机理进行阐述。首先从最小均方差准则出发, 推导获得维纳预测器公式, 在此基础上从时间、空间相位差角度分析, 得出基于维纳预测器的非同步测量数据同步化步骤, 实现将非同步测量数据同步至任一测量时刻, 然后进一步推导获得互谱矩阵补全方法。通过理论推导、仿真计算和试验, 研究了维纳预测器进行非同步测量数据同步化的机理, 并验证了所提基于维纳预测器数据同步化处理方法的准确性和有效性。
1. 问题描述
在理想流体介质中的小振幅声波满足三维波动方程:
{\nabla ^2}p\left( {{\boldsymbol{r}}, t} \right) = \frac{1}{{c_0^2}}\frac{{{\partial ^2}p\left( {{\boldsymbol{r}}, t} \right)}}{{\partial {t^2}}}, (1) 其中, p(r, t)表示t时刻在r处的声压, c0表示声速,
{{\nabla}}^{{2}} 为拉普拉斯算子。由于波动方程是线性方程, 满足叠加原理, 转换至频域分析:{\nabla ^2}p\left( {\boldsymbol{r}} \right) + {k^2}p\left( {\boldsymbol{r}} \right) = 0, (2) 其中, k = ω/c0为波数, ω为角频率, p(r) = p(r, ω)为频域声压。考虑空间中存在声源Q, 声源体积为V, 辐射声场中点r处的声压可表示为
p\left( {\boldsymbol{r}} \right) = \int_V {Q\left( {{{\boldsymbol{r}}_{{0}}}} \right)g\left( {{\boldsymbol{r, }}{{\boldsymbol{r}}_{{0}}}} \right)} {{{\mathrm{d}}}}V\left( {{{\boldsymbol{r}}_{{0}}}} \right), (3) 其中, g(r, r0)为三维空间中的格林函数。
考虑无噪声干扰条件, 对声源进行离散建模[30], 假定有N个离散声源, 可视为k个不相干声源
{{\boldsymbol{q}}}{=} {\left[{{q}}_{{1}}{, }{{q}}_{{2}}{, \cdots, }{{q}}_{{k}}\right]}^{{{\mathrm{T}}}} , 声场中有带m个传声器的测量阵列, 数量u不小于k的参考传声器。参考传声器声压为{{\boldsymbol{p}}_{{{\mathrm{R}}}}} = {{\boldsymbol{H}}_{{{\mathrm{R}}}}}{\boldsymbol{q}}, (4) 其中,
{{\boldsymbol{p}}_{\rm R}} = {\left[{{p}}_{{1}}\left({{{\boldsymbol{r}}}}_{{1}}\right){, }{{p}}_{{2}}\left({{{\boldsymbol{r}}}}_{{2}}\right){, \cdots , }{{p}}_{{u}}\left({{{\boldsymbol{r}}}}_{{u}}\right)\right]}^{{{\mathrm{T}}}} 表示u个参考传声器的复声压列向量;{{{\boldsymbol{H}}}}_{{{\rm{R}}}}{=}{{}\{{{h}}_{{ij}}\}{{}}}_{{u \times k}}{\in}{{{\mathbb{C}}}}^{{u \times k}} 为声源到参考传声器的前向传播矩阵, hij为前向传播函数。由式(4), 有{\boldsymbol{q}} = {\boldsymbol{H}}_{{{\rm{R}}}}^ + {{\boldsymbol{p}}_{{{\rm{R}}}}}, (5) 其中,
{{{\boldsymbol{H}}}}_{{{\rm{R}}}}^{{ + }} 为HR的广义逆, 因为u ≥ k, 所以q有唯一解, 即q可由{{{\boldsymbol{p}}}}_{{{\rm{R}}}} 线性表示。考虑实际测量中存在误差, 结合离散后的波动方程, 有{\boldsymbol{p = }}\left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{p}}_{{{\rm{R}}}}}} \\ {{{\boldsymbol{p}}_{{{\rm{M}}}}}} \end{array}} \right]{\boldsymbol{ = Hq + n}} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{H}}_{{{\rm{R}}}}}} \\ {{{\boldsymbol{H}}_{{{\rm{M}}}}}} \end{array}} \right]{\boldsymbol{q + n}}, (6) 其中,
{{{\boldsymbol{p}}}}_{{{\rm{M}}}} = {\left[{{{\boldsymbol{p}}}}_{{u + }{1}}\left({{{\boldsymbol{r}}}}_{{u + }{1}}\right){, }{{p}}_{{u + }{2}}\left({{{\boldsymbol{r}}}}_{{u + }{2}}\right){, \cdots, }{{p}}_{{u + m}}\left({{{\boldsymbol{r}}}}_{{u + m}}\right)\right]}^{{{\rm{T}}}} 表示阵列的复声压列向量,{{\boldsymbol{p}}}{=}{{}\{{{p}}_{{i}}\}{{}}}_{{(}{u + m}{) \times 1}}{\in}{{{\mathbb{C}}}}^{{(}{u + m}{) \times 1}} 为u个参考传声器与阵列m个传声器的复声压向量,{{\boldsymbol{H}}}{=} {{}\{{{h}}_{{ij}}\}{{}}}_{{(}{u + m}{) \times }{k}}{\in}{{{\mathbb{C}}}}^{{(}{u + m}{) \times }{k}} 为声源到传声器的前向传播矩阵。单个传声器声压为{p_i} = \sum\limits_{j = 1}^k {{h_{ij}}{q_j}} + {n_i} = {\boldsymbol{H}}_i^{}{\boldsymbol{q}} + {n_i}. (7) {p_i} = {\boldsymbol{H}}_i^{}{\boldsymbol{H}}_{{{\rm{R}}}}^ + {{\boldsymbol{p}}_{{{\rm{R}}}}} + {n_i}, (8) 其中,
{{{\boldsymbol{H}}}}_{{i}}{=}{{}\{{{h}}_{{ij}}\}{{}}}_{{1}{ \times k}}{\in}{{\mathbb{C}}}^{{1 \times }{k}} 表示从声源到第i个传声器的前向传播矩阵, ni为第i个传声器声压的相关噪声。令{\boldsymbol{G}}_i^{} = {\boldsymbol{H}}_i^{}{\boldsymbol{H}}_{{{\rm{R}}}}^ + , (9) 代入式(8), 得第i个传声器所测声压为
{p_i} = {\boldsymbol{G}}_i^{}{{\boldsymbol{p}}_{{{\rm{R}}}}} + {n_i}. (10) 其期望输出为
{\widehat p_i} = {\boldsymbol{\widehat G}}_i^{}{{\boldsymbol{p}}_{{{\rm{R}}}}}, (11) 其中, Gi表示由参考传声器到第i个传声器的传递矩阵,
{\widehat{{\boldsymbol G}}}_{{i}} 为{\widehat{{p}}}_{{i}} 所对应的由参考传声器到第i个传声器的传递矩阵。在测量过程中, 由于非同步测量导致相位丢失, 需要对阵列数据进行同步化处理。2. 非同步数据的维纳预测器同步化处理
2.1 维纳预测器推导
为使
{\widehat{{p}}}_{{i}} 尽可能接近原始信号{{p}}_{{i}} , 采用最小均方误差分析:e = \left|{p_i} - {\widehat p_i} \right|, (12) \mathop {\arg \min }\limits_{{{\widehat p}_i}} \mathbb{E}\left\{ {{e^2}} \right\} = \mathop {\arg \min }\limits_{{{\widehat p}_i}} \mathbb{E}\left\{ {{{({p_i} - {{\widehat p}_i})}^{{{\rm{H}}}}}({p_i} - {{\widehat p}_i})} \right\}. (13) \mathop {\arg \min }\limits_{{\boldsymbol{\widehat G}}_i^{}} \mathbb{E}\left\{ {{{({p_i} - {{\widehat p}_i})}^{{{\rm{H}}}}}({p_i} - {{\widehat p}_i})} \right\}, (14) 式(14)对
{\widehat{{{\boldsymbol{G}}}}}_{{i}} 求导:\begin{split}& \frac{{\partial \mathbb{E}\left\{ {{e^2}} \right\}}}{{\partial {\boldsymbol{\widehat G}}_i^{}}} = \mathbb{E}\left\{ {{{({p_i} - {\boldsymbol{\widehat G}}_i^{}{{\boldsymbol{p}}_{{{\rm{R}}}}})}^{{{\rm{H}}}}}({p_i} - {\boldsymbol{\widehat G}}_i^{}{{\boldsymbol{p}}_{{{\rm{R}}}}})} \right\} =\\&\quad \frac{\partial }{{\partial {\boldsymbol{\widehat G}}_i^{}}}\mathbb{E}\left\{ {p_i^{{{\rm{H}}}}{p_i} - p_i^{{{\rm{H}}}}{\boldsymbol{\widehat G}}_i^{}{{\boldsymbol{p}}_{{{\rm{R}}}}} - {\boldsymbol{p}}_{{{\rm{R}}}}^{{{\rm{H}}}}{{\left( {{\boldsymbol{\widehat G}}_i^{}} \right)}^{{{\rm{H}}}}}{p_i} + {\boldsymbol{p}}_{{{\rm{R}}}}^{{{\rm{H}}}}{{\left( {{\boldsymbol{\widehat G}}_i^{}} \right)}^{{{\rm{H}}}}}{\boldsymbol{\widehat G}}_i^{}{{\boldsymbol{p}}_{{{\rm{R}}}}}} \right\}. \end{split} (15) 根据矩阵变元求导法则以及极值点与导数的关系, 进一步得到
\frac{{\partial \mathbb{E}\left\{ {{e^2}} \right\}}}{{\partial {\boldsymbol{\widehat G}}_i^{}}} = 2\mathbb{E}\left\{ {{\boldsymbol{\widehat G}}_i^{}{{\boldsymbol{p}}_{{{\rm{R}}}}}{\boldsymbol{p}}_{{{\rm{R}}}}^{{{\rm{H}}}} - {p_i}{\boldsymbol{p}}_{{{\rm{R}}}}^{{{\rm{H}}}}} \right\} = 2\left( {{\boldsymbol{\widehat G}}_i^{}\mathbb{E}\left\{ {{{\boldsymbol{p}}_{{{\rm{R}}}}}{\boldsymbol{p}}_{{{\rm{R}}}}^{{{\rm{H}}}}} \right\} - \mathbb{E}\left\{ {{p_i}{\boldsymbol{p}}_{{{\rm{R}}}}^{{{\rm{H}}}}} \right\}} \right) = 0. (16) 由式(16), 得
{\boldsymbol{\widehat G}}_i^{}{{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{{\rm{R}}}}}{{\boldsymbol{p}}_{{{\rm{R}}}}}}} - {{\boldsymbol{S}}_{{p_i}{{\boldsymbol{p}}_{{{\rm{R}}}}}}} = 0, (17) {\boldsymbol{\widehat G}}_i^{} = {{\boldsymbol{S}}_{{p_i}{{\boldsymbol{p}}_{{{\rm{R}}}}}}}{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{{\rm{R}}}}}{{\boldsymbol{p}}_{{{\rm{R}}}}}}^{ - 1}, (18) 其中,
{\widehat{{{\boldsymbol{G}}}}}_{{i}} 即为用于第i个传声器信号预测的维纳预测器;{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{{\rm{R}}}}{{{\boldsymbol{p}}}}_{{{\rm{R}}}}}{=}{{\mathbb{E}}}{{}\{{{{\boldsymbol{p}}}}_{{{\rm{R}}}}{{{\boldsymbol{p}}}}_{{{\rm{R}}}}^{{{\rm{H}}}}\}{}} 表示参考传声器信号的互谱矩阵,{{{\boldsymbol{S}}}}_{{{p}}_{{i}}{{{\boldsymbol{p}}}}_{{{\rm{R}}}}}{=}{\mathbb E}{{}\{{{p}}_{{i}}{{{\boldsymbol{p}}}}_{{{\rm{R}}}}^{{{\rm{H}}}}\}{}} 表示第i个传声器的输出信号与参考传声器信号的互谱矩阵。2.2 数据同步化
假定经过d次非同步测量, 测量数据为p1, p2, ···, pd, 其中
{{\boldsymbol{p}}^i}{\boldsymbol{ = }}\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{p}}_{{{\rm{R}}}}^i} \\[1.5mm] {{\boldsymbol{p}}_{{{\rm{M}}}}^i} \end{array}} \right], (19) 包括第i次测量的参考传声器声压信号
{{{\boldsymbol{p}}}}_{{{\rm{R}}}}^{{i}} 和传声器阵列声压信号{{{\boldsymbol{p}}}}_{{{\rm{M}}}}^{{i}} 。{{{\boldsymbol{p}}}}_{{{\rm{R}}}}^{{i}} 与{{{\boldsymbol{p}}}}_{{{\rm{M}}}}^{{i}} 为同步测量的数据, 两组数据之间仅存在空间相位差;{{{\boldsymbol{p}}}}_{{{\rm{R}}}}^{{i}} 与{{{\boldsymbol{p}}}}_{{{\rm{R}}}}^{{j}} 为非同步测量, 但两次测量空间位置不变, 因此不存在空间相位差, 仅存在时间相位差。每次测量时阵列在空间和时间上均不同步, 故存在相位丢失问题, 需要进行相位补偿。根据式(11)和式(18), 有{\boldsymbol{p}}_{{\rm M}}^i = {\boldsymbol{G}}_{{{{M}}^i}}^{}{\boldsymbol{p}}_{_{{\rm R}}}^i = {{\boldsymbol{S}}_{{\boldsymbol{p}}_{{\rm M}}^i{\boldsymbol{p}}_{{\rm R}}^i}}{\boldsymbol{S}}_{{\boldsymbol{p}}_{{\rm R}}^i{\boldsymbol{p}}_{{\rm R}}^i}^{ - 1}{\boldsymbol{p}}_{{\rm R}}^i, (20) 这表明维纳预测器
{\boldsymbol{G}}_{{{{M}}^i}}^{}{{ = }}{{\boldsymbol{S}}_{{\boldsymbol{p}}_{{\rm M}}^i{\boldsymbol{p}}_{_{{\rm R}}}^i}}{\boldsymbol{S}}_{{\boldsymbol{p}}_{_{{\rm R}}}^i{\boldsymbol{p}}_{_{{\rm R}}}^i}^{ - 1} (21) 具有空间相位差补偿作用。因为测量过程中参考传声器位置不变, 故
{{{\boldsymbol{p}}}}_{{\rm R}}^{{i}} 与{{{\boldsymbol{p}}}}_{{\rm M}}^{{i}} 之间空间相位差等于{{{\boldsymbol{p}}}}_{{\rm R}}^{{j}} 与{{{\boldsymbol{p}}}}_{{\rm M}}^{{i}} 之间的空间相位差, 即{{{\boldsymbol{p}}}}_{{\rm M}}^{{i}} 与{{{\boldsymbol{p}}}}_{{\rm R}}^{{j}} 之间的空间相位差可通过式(21)所示维纳预测器进行补偿, 从而实现数据同步化处理, 将第i次测量的阵列数据与第j次测量的阵列数据进行同步, 有{\boldsymbol{p}}_{{\rm M}}^{i, j}{\boldsymbol{ = }}{{\boldsymbol{S}}_{{\boldsymbol{p}}_{{\rm M}}^i{\boldsymbol{p}}_{{\rm R}}^i}}{\boldsymbol{S}}_{{\boldsymbol{p}}_{{\rm R}}^i{\boldsymbol{p}}_{{\rm R}}^i}^{{{ - 1}}}{\boldsymbol{p}}_{{\rm R}}^j, (22) 式中,
{{{\boldsymbol{p}}}}_{{\rm M}}^{{i, j}} 右上标i表示原始数据序贯序号, j表示将原始数据同步至第j次测量。由式(22)可将所有测量数据相位与第j次测量相位进行同步, 此时所有数据之间不存在相位丢失。从最小均方误差出发, 推导出利用相关性进行数据同步化的维纳预测器, 表明在非同步测量过程中, 维纳预测器可以利用阵列与参考传声器信号之间的互谱关系, 传递声波传播过程中时间、空间的相位关联性, 通过维纳预测器对时空相位差进行补偿, 实现非同步测量数据的同步化处理。
2.3 互谱矩阵计算
为根据所测信息反演声源, 一般需要进行互谱矩阵求解, 文献[29]给出了一种不需要进行相位补全的互谱矩阵计算方法(详见附录A), 这里基于维纳预测器对该方法提供一种新的解释。
{{\boldsymbol{S}}_{{\boldsymbol{pp}}}}{{ = }}\mathbb{E}\left\{ {{\boldsymbol{p}}{{\boldsymbol{p}}^{{\rm H}}}} \right\} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}}^{}}&{{{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}}} \\[2mm] {{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}^{{\rm H}}}&{{{\boldsymbol{S}}_{{{MM}}}}} \end{array}} \right], (23) 其中, Spp表示参考传声器信号与所有非同步测量阵列信号构成的互谱矩阵,
{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm M}}} 表示每次测量中参考传声器信号与阵列信号的互谱矩阵, SMM表示所有非同步测量阵列数据的互谱矩阵。每次非同步测量阵列数据的互谱矩阵{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm M}}^{{i}}{{{\boldsymbol{p}}}}_{{\rm M}}^{{i}}} 可直接计算获得, 但因相位丢失的原因, 所有阵列数据所构成的互谱矩阵存在数据缺失, 如图1所示, 非同步测量阵列数据互谱矩阵SMM不完整, 除了对角块{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm M}}^{{i}}{{{\boldsymbol{p}}}}_{{\rm M}}^{{i}}} 以及每次阵列数据与参考传声器的互谱矩阵块{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm R}}} ,{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm M}}} ,{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm M}}}^{{\rm H}} , 其余数据均缺失, 可初始化为0。利用维纳预测器计算SMM:\begin{split} {{\boldsymbol{S}}_{{{{\rm{MM}}}}}} = & \mathbb{E}\{ {{\boldsymbol{p}}_{{\rm M}}}{\boldsymbol{p}}_{{\rm M}}^{{\rm H}}\} = \mathbb{E}\{ {\boldsymbol{G}}_{{\rm M}}^{}{{\boldsymbol{p}}_{{\rm R}}}{({\boldsymbol{G}}_{{\rm M}}^{}{{\boldsymbol{p}}_{{\rm R}}})^{{\rm H}}}\} = \\ & \mathbb{E}\{ {{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm M}}}{{\boldsymbol{p}}_{{\rm R}}}}}{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}}^{ - 1}{{\boldsymbol{p}}_{{\rm R}}}{{{(}}{{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm M}}}{{\boldsymbol{p}}_{{\rm R}}}}}{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}}^{ - 1}{{\boldsymbol{p}}_{{\rm R}}}{{)}}^{{\rm H}}}\} , \end{split} (24) \begin{split} {{\boldsymbol{S}}_{{{{\rm{MM}}}}}} = & {{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm M}}}{{\boldsymbol{p}}_{{\rm R}}}}}{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}}^{ - 1}\mathbb{E}\{ {{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}^{{\rm H}}\} {\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}}^{ - {1^{{\rm H}}}}{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm M}}}{{\boldsymbol{p}}_{{\rm R}}}}^{{\rm H}} =\\ & {{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm M}}}{{\boldsymbol{p}}_{{\rm R}}}}}{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}}^{ - {1^{{\rm H}}}}{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm M}}}{{\boldsymbol{p}}_{{\rm R}}}}^{{\rm H}}, \end{split} (25) 其中,
{\boldsymbol G}_{\rm M} 表示非同步测量参考传声器信号{{{\boldsymbol{p}}}}_{{\rm R}} 到阵列信号{{{\boldsymbol{p}}}}_{{\rm M}} 的传递矩阵,{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm R}}}^{-1} 为厄尔米特矩阵, 且{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm M}}{{{\boldsymbol{p}}}}_{{\rm R}}}{=}{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm M}}}^{{\rm H}} , 所以有{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm R}}}^{{-1}^{{\rm H}}}={{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm R}}}^{-1} ,{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm M}}{{{\boldsymbol{p}}}}_{{\rm R}}}^{\rm{H}}={{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm M}}} , 式(25)可进一步转换为{{\boldsymbol{S}}_{{{{\rm{MM}}}}}} = {\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}^{{\rm H}}{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}}^{ - 1}{{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}}. (26) 故可在获得互谱矩阵数据
{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm R}}} ,{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm M}}} ,{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm M}}}^{{\rm H}} 后, 对SMM数据进行补全, 从而获得完整互谱矩阵。式(26)与文献[29]所得结果一致, 即从最小均方误差角度出发, 基于维纳预测器对该结论给出了新的推导方法。3. 数值仿真
为验证所提出的基于维纳预测器的带有参考传声器的非同步测量数据同步化方法的准确性, 以稳态双点声源叠加声场以及板声源辐射声场非同步测量为例进行仿真分析, 采用基于等效源的近场声全息技术进行声场重建, 求逆过程采用Tikhonov正则化方法, 其中正则化参数采用L曲线法进行求解。
定义数据同步化误差为
{\rm{Error}}_{{{\rm{Syn}}}}{ = }\frac{{\Vert {{\boldsymbol{p}}}_{{{\rm{syn}}}}-{\displaystyle \underset{i{=1}}{\overset{d}{\cup }}{{\boldsymbol{p}}}_{{\rm M}}^{i, j}}\Vert }_{2}}{{\Vert {{\boldsymbol{p}}}_{{{\rm{syn}}}}\Vert }_{2}}{, } (27) 其中, psyn表示同步测量数据,
\bigcup _{{i=1}}^{{d}}{{\boldsymbol p}}_{{\rm M}}^{{i, j}} 表示所有同步化处理后的阵列数据合并,{\|{\cdot}\|}_{{2}} 表示二范数。定义重建误差为
{\rm Error}_{{\rm{re}}} = \frac{{{{\left\| {{{\boldsymbol{p}}_{{{{\rm{cal}}}}}} - {{\boldsymbol{p}}_{{{{\rm{mea}}}}}}} \right\|}_2}}}{{{{\left\| {{{\boldsymbol{p}}_{{{{\rm{mea}}}}}}} \right\|}_2}}}, (28) 其中, pcal表示计算所得重建值, pmea表示测量所得声压分布, 仿真与实验设定一致。
3.1 仿真一: 双点声源辐射声场
3.1.1 仿真参数设置
在(−0.125, 0, 0.055) m, (0.125, 0, 0.055) m处分别布置一个源强为
{{Q}}_{{0}}{=}{{10}}^{{-4}}\;{{{\rm{m}}}}^{{3}}{/}{{\rm{s}}} 的点声源, 2个点声源之间相位差保持恒定, 等效源面垂直z轴, 在z=0.075 m处, 大小为0.55 m × 0.35 m, 离散点数为12 × 8。传声器阵列尺寸为0.25 m × 0.35 m, 阵元数量为6 × 8个, 阵元间距为0.05 m, 全息面位于xOy平面, 对声场进行3次非同步测量, 阵列中心分别为(−0.125, 0, 0) m, (0, 0, 0) m, (0.175, 0, 0) m, 在(−0.6, 0, 0.25) m处布置一个参考传声器, 在Z=0.025 m处布置一个重建面, 大小为0.55 m × 0.35 m, 离散点数为12 × 8, 如图2所示。仿真过程声源相位差为恒定, 因此2个声源为相干声源, 从统计角度可视为1个声源[26]。3.1.2 仿真结果分析
(1)数据同步化分析
在信噪比为20 dB环境下, 对2010 Hz双点声源声场进行距离为0.055 m的3次非同步测量, 如图3所示, 对测量结果进行同步化处理, 如图4所示。该算例中, 通过式(22)将三次非同步测量数据与第1次测量的时间相位进行匹配, 第1次非同步测量数据与其同步化处理后数据相等, 而第2, 3次非同步测量数据与其同步化处理后的数据存在明显差异。数据同步化处理后相位发生变化, 因此在各个位置复声压幅值不变的情况下, 实声压发生变化, 第2, 3次测量值与同步化处理后的值存在差异。
经同步化处理后, 将3次测量数据根据坐标关系进行合并, 如图5所示, 对应同步测量数据如图6所示。非同步测量数据经同步化处理后能得到与同步测量接近的声压分布, 存在的差异主要由仿真过程中所添加随机噪声干扰所致, 由式(27)可得同步误差为10.29%。与图3和图4对比, 图5表明非同步测量可以获得更大的成像范围, 即通过非同步测量可以合成更大的阵列孔径, 从而覆盖更大的测量空间。
因空间和时间上的差异, 同一次测量不同传声器之间存在相位差异, 不同测量, 同一传声器相位也存在差异, 在信噪比为20 dB、声源频率为2010 Hz工况下第1~3次测量传声器1, 12, 25, 37, 48以及参考传声器的相位如表1所示。通过式(22)将三次非同步测量数据与第1次测量相位进行同步, 从表中可以看出第1次测量6个传声器相位差为零, 第2次、第3次相位变化不为零, 5个传声器相位变化与参考传声器相位变化一致, 表明在参考传声器作用下, 5个传声器进行时间相位调整。其中, 第3次非同步测量数据同步化后, 传声器1, 12的相位变化为0.9849, 参考传声器的相位变化为0.985, 误差为0.01%, 为随机误差干扰所致。在数据同步后, 同一传声器的相位不一致, 为其空间位置不同所致, 而参考传声器位置不变, 因此其同步后相位一致。
表 1 不同传声器在2010 Hz下同步前后的相位及其变化值传声器序号 参考传声器 传声器1 传声器12 传声器25 传声器37 传声器48 第1次
相位 (rad)同步前 0.8162 −0.6142 −2.4326 1.8952 −2.0315 −0.3048 同步后 0.8162 −0.6142 −2.4326 1.8952 −2.0315 −0.3048 相位差 0 0 0 0 0 0 第2次
相位 (rad)同步前 0.2294 1.1106 −1.4217 −0.9299 −1.5229 1.0363 同步后 0.8162 1.6974 −0.8349 −0.3431 −0.9361 1.6231 相位差 0.5868 0.5868 0.5868 0.5868 0.5868 0.5868 第3次
相位 (rad)同步前 −0.1688 −1.2845 −1.8122 −0.4965 1.1880 −3.0281 同步后 0.8162 −0.2996 −0.8273 0.4884 2.1730 −2.0431 相位差 0.9850 0.9849 0.9849 0.9849 0.9850 0.9850 在不同信噪比条件下, 根据式(27)计算三次非同步测量数据的同步化误差, 随声源频率变化如图7所示。从图中可知, 在110~2010 Hz范围内, 当无噪声干扰时, 同步化后的数据与同步测量理论值一致, 即误差为0, 且数据同步化误差不受频率变化影响, 验证了该方法的准确性。当信噪比为30 dB时, 数据同步化误差均维持在6%以下; 当信噪比为20 dB时, 数据同步化误差低于25%。随着信噪比降低, 误差增大, 表明数据同步化误差主要受仿真过程中随机噪声的影响, 信噪比越低, 随机成分越多, 同步误差越大。
在信噪比为30 dB条件下, 考虑不同测量距离对数据同步误差的影响, 依次设置测量距离为0.075 m, 0.065 m, 0.055 m, 0.045 m, 0.035 m, 根据式(27)计算数据同步化误差, 不同测量距离条件下, 数据同步化误差随频率变化曲线如图8所示。在30 dB下, 数据同步化误差随距离无明显变化, 均维持在较低水平, 表明该数据同步化方法在合理距离范围内不受距离影响, 误差主要受测量信号中随机成分的影响。
(2)重建结果分析
在测量距离为0.055 m条件下, 将采用单全息面第2次测量数据对重建面声压进行重建(方法1)以及采用3次非同步测量数据对重建面声压进行重建(方法2)的结果分别与理论值进行对比, 根据式(28)计算不同信噪比条件下方法1和方法2的重建误差, 如图9所示。对比可得, 方法2的重建误差远小于方法1的重建误差, 且在110~2010 Hz频率范围内, 重建误差波动受频率变化影响较小。分析原因: 方法1中, 点声源位于全息面边沿, 且测量点数小于等效源点数, 因此在测量孔径和方程数目上均存在不足, 故重建误差较大; 方法2中, 非同步测量一方面增大了测量范围, 能有效覆盖声源, 另一方面, 非同步测量方法增加了测量点数, 可以提供更多有效的声压数据, 避免方程欠定, 因此方法2能维持较高精度的测量要求。
(3)互谱矩阵计算分析
在测量距离为0.055 m条件下, 采用非同步测量方法获取3次单全息面数据后, 根据式(26)计算互谱矩阵。在不同信噪比条件下, 与理论值相比, 补全误差随频率变化如图10所示, 当无噪声干扰时, 补全误差为0, 表明该方法的准确性。随着信噪比降低, 受噪声干扰影响, 补全误差增大。
3.2 仿真二: 板声源辐射声场
3.2.1 仿真参数设置
在z = 0.055 m处放置一个平行于xOy平面的正方形简支铝板, 尺寸为0.5 m × 0.5 m, 厚度为0.003 m, 泊松比为0.32, 杨氏模量为7.1 × 1010 Pa, 密度为2.7 × 103 kg/m3。激励点位于板的中心, 激励力幅值为2 N, 等效源面垂直z轴, 在z = 0.075 m处, 大小为0.5 m × 0.5 m, 离散点数为19 × 19。传声器阵列尺寸为0.25 m × 0.35 m, 阵元数量为6 × 8个, 阵元间距为0.05 m, 全息面位于xOy平面, 以坐标原点为中心均匀布置测量阵列, 相邻2次测量间距为0.05 m, 对板声源声场进行9次非同步测量, 在(−0.6, 0, 0.25) m处布置一个参考传声器, 在z = 0.025 m处布置一个重建面, 大小为0.55 m × 0.35 m, 离散点数为12 × 8, 如图11所示。
3.2.2 仿真结果分析
(1)数据同步化分析
在信噪比为30 dB环境下, 对频率为2010 Hz的板声源辐射声场进行9次非同步测量, 如图12所示, 对测量结果进行同步化处理, 如图13所示。该算例中, 通过式(22)将9次非同步测量数据与第1次测量的时间相位进行匹配, 第1次非同步测量数据与其同步化处理后数据相等, 而第2~8次非同步测量数据与其同步化处理后的数据存在明显差异, 同步化处理后相位发生变化, 因此在各个位置复声压幅值不变的情况下, 实声压发生变化, 第2~8次测量值与同步化处理后的值存在差异。如图13所示, 经同步化处理的声压数据云图呈对称分布, 符合对称声源的辐射效果。
![]()
图 12 非同步测量结果(声压实部) (a) 第1次测量; (b) 第2次测量; (c) 第3次测量; (d) 第4次测量; (e) 第5次测量; (f) 第6次测量; (g) 第7次测量; (h) 第8次测量; (i) 第9次测量![]()
图 13 非同步测量数据同步化结果(声压实部) (a) 第1次测量; (b) 第2次测量; (c) 第3次测量; (d) 第4次测量; (e) 第5次测量; (f) 第6次测量; (g) 第7次测量; (h) 第8次测量; (i) 第9次测量经同步化处理后, 将9次测量数据根据坐标关系进行合并, 如图14所示, 对应同步测量数据如图15所示。对比可得, 非同步测量数据经同步化处理后能得到与同步测量接近的声压分布, 存在的差异主要由仿真过程中所添加随机噪声干扰所致, 由式(27)可得同步误差为7.30%。与图12和图13对比, 图14表明非同步测量可以获得更大的成像范围, 即通过非同步测量可以合成更大的阵列孔径, 从而覆盖更大的测量空间。
因空间和时间上的差异, 同一次测量不同传声器之间存在相位差异, 不同测量, 同一传声器相位也存在差异, 在信噪比为30 dB、声源频率为2010 Hz工况下第1, 5, 9次测量传声器2, 13, 26, 36, 47以及参考传声器的相位如表2所示。通过式(22)将三次非同步测量数据与第1次测量相位进行同步, 从表中可以看出第1次测量6个传声器相位差为零, 第5, 9次相位变化不为零, 其中, 第5次测量的5个传声器相位变化与参考传声器相位变化一致, 第9次测量的5个传声器相位差与参考传声器的相位差之和为2
\pi , 表明在参考传声器作用下, 5个传声器进行时间相位调整。其中, 第9次非同步测量数据同步化后, 传声器13, 26的相位变化为−2.3465, 参考传声器的相位变化为−2.3464 (3.9368+2.3464 = 2\pi ), 误差为0.004%, 为随机误差干扰所致。在数据同步后, 同一传声器的相位不一致, 为其空间位置不同所致, 而参考传声器位置不变, 因此其同步后相位一致。表 2 不同传声器在2010 Hz下同步前后的相位及其变化值传声器序号 参考传声器 传声器2 传声器13 传声器26 传声器36 传声器47 第1次
相位 (rad)同步前 1.0951 0.6194 −1.4823 −1.8059 −2.9863 0.3467 同步后 1.0951 0.6194 −1.4823 −1.8059 −2.9863 0.3467 相位差 0 0 0 0 0 0 第5次
相位 (rad)同步前 −0.9164 −2.4905 0.6970 0.0671 0.6933 −2.4456 同步后 1.0951 −0.4790 2.7085 2.0786 2.7048 −0.4341 相位差 2.0115 2.0115 2.0115 2.0115 2.0115 2.0115 第9次
相位 (rad)同步前 −2.8417 2.6326 −0.7051 −0.5852 0.9455 2.9820 同步后 1.0951 0.2862 −3.0516 −2.9317 −1.4010 0.6356 相位差 3.9368 −2.3464 −2.3465 −2.3465 −2.3465 −2.3464 在不同信噪比条件下, 根据式(27)计算9次非同步测量数据的同步化误差, 其随声源频率的变化如图16所示。从图中可知, 在110~2010 Hz范围内, 当无噪声干扰时, 同步化后的数据与同步测量理论值一致, 即误差为0, 且数据同步化误差不受频率变化影响, 验证了该方法的准确性, 表明在复杂声源辐射声场中, 该方法仍能满足非同步测量数据同步化处理需求。当信噪比为30 dB时, 除410 Hz工况外, 其余数据同步化误差均维持在10%以下; 当信噪比为20 dB时, 数据同步化误差低于25%。随着信噪比降低, 误差增大, 表明数据同步化误差主要受仿真过程中随机噪声的影响, 信噪比越低, 随机成分越多, 同步误差越大。
对比双点声源辐射声场以及板声源辐射声场的非同步测量数据同步化结果, 从理论上分析, 数据同步化误差不受声源形式影响, 因此在无干扰条件下, 双点声源和板声源辐射声场的数据同步化误差均为0。在信噪比为20 dB和30 dB条件下, 双点声源声场的数据同步化误差整体趋势小于板声源声场的同步化误差, 主要受随机噪声的干扰。在有噪声干扰条件下, 采用同样的阵列, 板声源声场进行了9次非同步测量, 而双点声源声场仅进行3次非同步测量, 因此板声源声场累积的随机干扰成分更多, 故其数据同步化整体误差更大。
(2)重建结果分析
将采用单全息面第5次测量数据对重建面声压进行重建(方法3)以及采用9次非同步测量数据对重建面声压进行重建(方法4)的结果分别与理论值进行对比, 根据式(28)计算不同信噪比条件下方法3和方法4的重建误差, 如图17所示。对比可得, 方法4的重建误差远小于方法3的重建误差, 且在110~2010 Hz频率范围内, 重建误差波动受频率变化影响较小。分析原因: 方法3中, 单次测量全息面尺寸小于板声源尺寸, 小于重建面尺寸, 且测量点数小于等效源点数, 因此在测量孔径和方程数目上均存在不足, 故重建误差较大; 随着频率增大, 误差急剧增大, 表明测量点数不足, 严重影响高频段重建精度; 方法4中, 非同步测量一方面增大了测量范围, 能有效覆盖声源, 另一方面通过非同步测量, 增大了测量点数, 避免方程欠定, 从而可以有效减小重建误差。
(3)互谱矩阵计算分析
采用非同步测量方法获取9次单全息面数据后, 根据式(26)计算互谱矩阵。在不同信噪比条件下, 与理论值相比, 补全误差随频率变化如图18所示, 当无噪声干扰时, 补全误差为0, 表明该方法的准确性。随着信噪比降低, 随机成分增多, 受噪声干扰影响, 补全误差增大。
4. 实验验证
为进一步验证本文所提非同步测量数据同步化方法的有效性, 分别采用单次测量对重建面声压进行重建与三次非同步测量对重建面声压进行重建, 重建面距离声源0.03 m。实验在半消声室中进行, 选取两个相同的音箱作为声源, 播放含有1010 Hz, 1210 Hz, 1410 Hz, 1610 Hz, 1810 Hz, 2010 Hz共6个频率成分的谐波信号, 实验布置和测量过程如图19所示。两个扬声器的声源中心分别位于(−0.125, 0, 0.055) m和(0.125, 0, 0.055) m处, 置于1 m × 2 m吸声棉上, 以减小或消除反射对传声器阵列测量的影响, 从而满足自由场条件。传声器阵列尺寸为0.25 m × 0.35 m, 阵元数量为6 × 8个, 阵元间距为0.05 m, 全息面位于xOy平面, 对声场进行3次非同步测量, 阵列中心分别为(−0.125, 0, 0) m, (0, 0, 0) m, (0.175, 0, 0) m, 在(−0.6, 0, 0.25) m处布置1个参考传声器。实验过程声源相位差恒定, 因此2个声源为相干声源, 从统计角度可视为1个声源[26]。
传声器阵列通道1与参考传声器在3次非同步测量中时域波形如图20所示, 转换至频域, 幅值如图21(a)所示, 声源信号为谐波信号, 所以幅频图在1010 Hz, 1210 Hz, 1410 Hz, 1610 Hz, 1810 Hz, 2010 Hz六个分析频率处存在明显线谱。相应频率处的相位如图21(b)所示, 因为时间和空间位置上的差异, 所以三次非同步测量传声器1的幅值和相位均存在明显差异。因为时间的差异, 参考传声器的相位明显不一致; 幅值理论上应一致, 但因测量误差存在, 存在较小差异。
因空间和时间上的差异, 同一次测量不同传声器之间存在相位差异, 不同测量, 同一传声器相位也存在差异, 第1~3次测量传声器1以及参考传声器的相位如表3和表4所示。通过式(22)将三次非同步测量数据与第1次测量相位进行同步, 从表中可以看出第1次测量2个传声器相位差为零, 第2次、第3次相位变化不为零, 且相同频率下, 传声器1相位变化与参考传声器相位变化一致, 表明在参考传声器作用下, 传声器1进行时间相位调整。在数据同步后, 同一频率下传声器1的相位不一致, 为其空间位置不同所致, 而参考传声器位置不变, 因此其同步后相位一致。
表 3 传声器1不同频率下同步前后的相位及其变化值频率 (Hz) 1010 1210 1410 1610 1810 2010 第1次
相位 (rad)同步前 −1.0265 2.0688 −0.7234 2.8355 −0.5824 2.4639 同步后 −1.0265 2.0688 −0.7234 2.8355 −0.5824 2.4639 相位差 0 0 0 0 0 0 第2次
相位 (rad)同步前 −0.5541 0.4942 1.8092 −2.9202 −0.7091 0.4154 同步后 −0.2406 2.9902 −0.1096 −2.8706 1.3593 −1.8306 相位差 −0.3135 −2.4960 1.9188 −0.0496 −2.0683 2.2459 第3次
相位 (rad)同步前 0.1539 −0.9385 1.0209 0.4902 −0.3722 −0.9339 同步后 −2.3767 0.0884 −0.5297 2.6740 −1.0746 2.2049 相位差 2.5306 −1.0269 1.5506 −2.1838 0.7024 −3.1388 表 4 参考传声器不同频率下同步前后的相位及其变化值频率 (Hz) 1010 1210 1410 1610 1810 2010 第1次
相位 (rad)同步前 −2.5066 −0.1400 1.4457 −3.1069 −0.9973 0.5982 同步后 −2.5066 −0.1400 1.4457 −3.1069 −0.9973 0.5982 相位差 0 0 0 0 0 0 第2次
相位 (rad)同步前 −2.8201 −2.6360 −2.9187 3.1267 −3.0656 2.8441 同步后 −2.5066 −0.1400 1.4457 −3.1069 −0.9973 0.5982 相位差 −0.3135 −2.4960 −4.3644 6.2336 −2.0683 2.2459 第3次
相位 (rad)同步前 0.0239 −1.1669 2.9963 0.9926 −0.2949 −2.5406 同步后 −2.5066 −0.1400 1.4457 −3.1069 −0.9973 0.5982 相位差 2.5305 −1.0269 1.5506 4.0995 0.7024 −3.1388 在完成非同步测量数据同步处理后, 分别采用单全息面2数据和三次非同步测量数据对重建面声压进行重建, 1010 Hz, 1410 Hz, 1810 Hz的重建结果与测量结果分别如图22—图24所示, 采用非同步测量对重建面声压进行重建所得结果更接近测量值。
采用单全息面数据和三次非同步测量数据进行重建面声压重建的误差曲线如图25所示。采用非同步测量数据进行声压重建误差远小于采用单全息面数据进行重建的误差, 采用非同步测量数据进行声场重建的误差均维持在15%以下, 采用单全息面进行声场重建的误差均大于20%, 且随频率变化波动较大。这表明了本文所提非同步测量数据同步化方法的准确性。采用非同步测量方法可以扩大测量阵列孔径和密度, 有利于提高声场重建精度, 且能在更宽的频率范围内提供准确的声场重建结果。从误差趋势分析, 采用单全息面1、单全息面3和非同步测量数据进行声场重建的误差均随频率的增大而增大, 而采用单全息面2数据进行声场重建的误差则先减小后增大, 与仿真趋势不一致, 与其余误差曲线趋势不一致。分析原因: (1)采用单全息2面数据进行重建, 单全息面不能完全覆盖声源区域, 且测量点数少于等效源数目, 因此重建误差较大, 且不稳定, 故出现误差曲线异常; (2)实验仅进行了单次测量, 未能排除为实验过程随机误差所致的误差曲线异常。
5. 结论
从最小均方差准则出发, 通过理论推导给出了用于非同步测量数据同步化处理的维纳预测器公式, 分析维纳预测器的相位补全作用, 结果表明维纳预测器可以通过参考传声器与阵列之间的空间相位关联性以及参考传声器之间的时间相位关联性, 对非同步测量数据进行同步化处理, 从而对非同步测量的机理进行阐述, 且实现将非同步测量数据同步至任一测量时刻。通过数值仿真和实验研究, 计算非同步测量数据同步化误差、互谱矩阵补全误差以及声场重建误差, 验证了该数据同步化方法的准确性和有效性。利用非同步测量数据可以扩大阵列孔径以及提高阵列密度, 在声场重建中提高声场重建精度。
附录A 互谱矩阵计算方法[29]
将复声压向量p分为由参考传声器测量所得的复声压列向量
{{{\boldsymbol{p}}}}_{{\rm R}} 以及由移动阵列测量所得的复声压列向量{{{\boldsymbol{p}}}}_{{\rm M}} 。为了对声源进行重建, 需要求解互谱矩阵Spp:{{\boldsymbol{S}}_{{\boldsymbol{pp}}}}{{ = }}\mathbb{E}\left\{ {{\boldsymbol{p}}{{\boldsymbol{p}}^{{\rm H}}}} \right\} = \left[ {\begin{array}{*{20}{c}} {{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}}^{}}&{{{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}}} \\ {{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}^{{\rm H}}}&{{{\boldsymbol{S}}_{{{MM}}}}} \end{array}} \right], (A1) 其中,
{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm R}}}{=}{{\mathbb{E}}}{{}\{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm R}}^{{\rm H}}\}{}} ,{{{\boldsymbol{S}}}}_{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm M}}}{=}{{\mathbb{E}}}{{}\{{{{\boldsymbol{p}}}}_{{\rm R}}{{{\boldsymbol{p}}}}_{{\rm M}}^{{\rm H}}\}{}} ,{{{\boldsymbol{S}}}}_{{{\mathrm{MM}}}}{=}{{\mathbb{E}}}{{}\{{{{\boldsymbol{p}}}}_{{\rm M}}{{{\boldsymbol{p}}}}_{{\rm M}}^{{\rm H}}\}{}} 。令{{\boldsymbol{S}}_1} = \left[ \begin{gathered} {\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}} \\[1.5mm] {\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}^{{\rm H}} \\ \end{gathered} \right], ~~ {{\boldsymbol{S}}_2} = \left[ \begin{gathered} {{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}} \\[1.5mm] {{\boldsymbol{S}}_{{{{\mathrm{MM}}}}}} \\ \end{gathered} \right], (A2) 则有
{{\boldsymbol{S}}_{{\boldsymbol{pp}}}} = \left[ {{{\boldsymbol{S}}_1}{{ }}{{\boldsymbol{S}}_2}} \right]. (A3) 假定Spp与S1秩相等, 则S2可由S1线性表示, 因此存在矩阵T, 使得
{{\boldsymbol{S}}_2} = {{\boldsymbol{S}}_1}{\boldsymbol{T}}, (A4) {{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}} = {\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}}^{}{\boldsymbol{T}}, {{\boldsymbol{S}}_{{{MM}}}} = {\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}^{{\rm H}}{\boldsymbol{T}}, (A5) 则有
{{\boldsymbol{S}}_{{{{\mathrm{MM}}}}}} = {\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}^{{\rm H}}{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm R}}}}^{ - 1}{{\boldsymbol{S}}_{{{\boldsymbol{p}}_{{\rm R}}}{{\boldsymbol{p}}_{{\rm M}}}}}. (A6) 当且仅当
{{S}}_{{{p}}_{{\rm R}}{{p}}_{{\rm R}}} 为满秩时, 式(A6)成立, 否则须用矩阵广义逆代替矩阵逆。 -
![]()
图 12 非同步测量结果(声压实部) (a) 第1次测量; (b) 第2次测量; (c) 第3次测量; (d) 第4次测量; (e) 第5次测量; (f) 第6次测量; (g) 第7次测量; (h) 第8次测量; (i) 第9次测量
![]()
图 13 非同步测量数据同步化结果(声压实部) (a) 第1次测量; (b) 第2次测量; (c) 第3次测量; (d) 第4次测量; (e) 第5次测量; (f) 第6次测量; (g) 第7次测量; (h) 第8次测量; (i) 第9次测量
表 1 不同传声器在2010 Hz下同步前后的相位及其变化值
传声器序号 参考传声器 传声器1 传声器12 传声器25 传声器37 传声器48 第1次
相位 (rad)同步前 0.8162 −0.6142 −2.4326 1.8952 −2.0315 −0.3048 同步后 0.8162 −0.6142 −2.4326 1.8952 −2.0315 −0.3048 相位差 0 0 0 0 0 0 第2次
相位 (rad)同步前 0.2294 1.1106 −1.4217 −0.9299 −1.5229 1.0363 同步后 0.8162 1.6974 −0.8349 −0.3431 −0.9361 1.6231 相位差 0.5868 0.5868 0.5868 0.5868 0.5868 0.5868 第3次
相位 (rad)同步前 −0.1688 −1.2845 −1.8122 −0.4965 1.1880 −3.0281 同步后 0.8162 −0.2996 −0.8273 0.4884 2.1730 −2.0431 相位差 0.9850 0.9849 0.9849 0.9849 0.9850 0.9850 
下载: 导出CSV
表 2 不同传声器在2010 Hz下同步前后的相位及其变化值
传声器序号 参考传声器 传声器2 传声器13 传声器26 传声器36 传声器47 第1次
相位 (rad)同步前 1.0951 0.6194 −1.4823 −1.8059 −2.9863 0.3467 同步后 1.0951 0.6194 −1.4823 −1.8059 −2.9863 0.3467 相位差 0 0 0 0 0 0 第5次
相位 (rad)同步前 −0.9164 −2.4905 0.6970 0.0671 0.6933 −2.4456 同步后 1.0951 −0.4790 2.7085 2.0786 2.7048 −0.4341 相位差 2.0115 2.0115 2.0115 2.0115 2.0115 2.0115 第9次
相位 (rad)同步前 −2.8417 2.6326 −0.7051 −0.5852 0.9455 2.9820 同步后 1.0951 0.2862 −3.0516 −2.9317 −1.4010 0.6356 相位差 3.9368 −2.3464 −2.3465 −2.3465 −2.3465 −2.3464
下载: 导出CSV
表 3 传声器1不同频率下同步前后的相位及其变化值
频率 (Hz) 1010 1210 1410 1610 1810 2010 第1次
相位 (rad)同步前 −1.0265 2.0688 −0.7234 2.8355 −0.5824 2.4639 同步后 −1.0265 2.0688 −0.7234 2.8355 −0.5824 2.4639 相位差 0 0 0 0 0 0 第2次
相位 (rad)同步前 −0.5541 0.4942 1.8092 −2.9202 −0.7091 0.4154 同步后 −0.2406 2.9902 −0.1096 −2.8706 1.3593 −1.8306 相位差 −0.3135 −2.4960 1.9188 −0.0496 −2.0683 2.2459 第3次
相位 (rad)同步前 0.1539 −0.9385 1.0209 0.4902 −0.3722 −0.9339 同步后 −2.3767 0.0884 −0.5297 2.6740 −1.0746 2.2049 相位差 2.5306 −1.0269 1.5506 −2.1838 0.7024 −3.1388
下载: 导出CSV
表 4 参考传声器不同频率下同步前后的相位及其变化值
频率 (Hz) 1010 1210 1410 1610 1810 2010 第1次
相位 (rad)同步前 −2.5066 −0.1400 1.4457 −3.1069 −0.9973 0.5982 同步后 −2.5066 −0.1400 1.4457 −3.1069 −0.9973 0.5982 相位差 0 0 0 0 0 0 第2次
相位 (rad)同步前 −2.8201 −2.6360 −2.9187 3.1267 −3.0656 2.8441 同步后 −2.5066 −0.1400 1.4457 −3.1069 −0.9973 0.5982 相位差 −0.3135 −2.4960 −4.3644 6.2336 −2.0683 2.2459 第3次
相位 (rad)同步前 0.0239 −1.1669 2.9963 0.9926 −0.2949 −2.5406 同步后 −2.5066 −0.1400 1.4457 −3.1069 −0.9973 0.5982 相位差 2.5305 −1.0269 1.5506 4.0995 0.7024 −3.1388
下载: 导出CSV
-
[1] Feng J, Liu Y, Li K M. Efficient frequency domain method for holographic localization of rotating sound sources with arbitrary array configurations. J. Sound Vib., 2023; 544: 117374 DOI: 10.1016/j.jsv.2022.117374
[2] Bi C X, Zhang F M, Zhang X Z, et al. Sound field reconstruction using block sparse Bayesian learning equivalent source method. J. Acoust. Soc. Am., 2022; 151(4): 2378−2390 DOI: 10.1121/10.0010103
[3] Zhang X Z, Bi C X, Zhang Y B, et al. Real-time nearfield acoustic holography for reconstructing the instantaneous surface normal velocity. Mech. Syst. Signal Process., 2015; 52-53: 663−671 DOI: 10.1016/j.ymssp.2014.06.009
[4] Xiao Y, Gao W, Chu C, et al. Identification of panel acoustic contribution based on patch nearfield acoustic holography. J. Phys. Conf. Ser., 2021; 1965(1): 012138 DOI: 10.1088/1742-6596/1965/1/012138
[5] Williams E G, Maynard J D, Skudrzyk E. Sound source reconstructions using a microphone array. J. Acoust. Soc. Am., 1980; 68(1): 340−344 DOI: 10.1121/1.384602
[6] Veronesi W A, Maynard J D. Digital holographic reconstruction of sources with arbitrarily shaped surfaces. J. Acoust. Soc. Am., 1989; 85(2): 588−598 DOI: 10.1121/1.397583
[7] Koopmann G H, Song L, Fahnline J B. A method for computing acoustic fields based on the principle of wave superposition. J. Acoust. Soc. Am., 1989; 86(6): 2433−2438 DOI: 10.1121/1.398450
[8] Valdivia N P, Williams E G. Study of the comparison of the methods of equivalent sources and boundary element methods for near-field acoustic holography. J. Acoust. Soc. Am., 2006; 120(6): 3694−3705 DOI: 10.1121/1.2359284
[9] Maynard J D, Williams E G, Lee Y. Nearfield acoustic holography: I. Theory of generalized holography and the development of NAH. J. Acoust. Soc. Am., 1985; 78(4): 1395−1413 DOI: 10.1121/1.392911
[10] Pinho M E V, Arruda J R F. On the use of the equivalent source method for nearfield acoustic holography. ABCM Symposium Series in Mechatronics, 2004: 590−599
[11] Sarkissian A. Extension of measurement surface in near-field acoustic holography. J. Acoust. Soc. Am., 2004; 115(4): 1593−1596 DOI: 10.1121/1.1645609
[12] 徐亮, 毕传兴, 陈剑, 等. 基于波叠加法的patch近场声全息及其实验研究. 物理学报, 2007; 56(5): 2776−2783 DOI: 10.3321/j.issn:1000-3290.2007.05.051 [13] 张小正, 毕传兴, 徐亮, 等. 基于波叠加法的近场声全息空间分辨率增强方法. 物理学报, 2010; 59(8): 5564−5571 DOI: 10.7498/aps.59.5564 [14] Wang J X, Zhang Z F, Huang Y Z, et al. A 3D convolutional neural network based near-field acoustical holography method with sparse sampling rate on measuring surface. Measurement, 2021; 177: 109297 DOI: 10.1016/j.measurement.2021.109297
[15] Wang J, Zhang Z, Li Z, et al. Research on joint training strategy for 3D convolutional neural network based near-field acoustical holography with optimized hyperparameters. Measurement, 2022; 202: 111790 DOI: 10.1016/j.measurement.2022.111790
[16] 胡定玉, 刘馨悦, 李曷冰, 等. 采用稀疏测量的高分辨率相干声场分离. 声学学报, 2020; 45(4): 563−570 DOI: 10.15949/j.cnki.0371-0025.2020.04.013 [17] Hu D Y, Li H B, Hu Y, et al. Sound field reconstruction with sparse sampling and the equivalent source method. Mech. Syst. Signal Process., 2018; 108: 317−325 DOI: 10.1016/j.ymssp.2018.02.031
[18] Bi C X, Liu Y, Xu L, et al. Sound field reconstruction using compressed modal equivalent point source method. J. Acoust. Soc. Am., 2017; 141(1): 73−79 DOI: 10.1121/1.4973567
[19] Zhang F M, Zhang X Z, Zhang Y B, et al. Sound field reconstruction using sparse Bayesian learning equivalent source method with hyperparametric-coupled prior. Appl. Acoust., 2023; 211: 109496 DOI: 10.1016/j.apacoust.2023.109496
[20] 杜向华, 朱海潮, 毛荣富, 等. 基于支持向量回归的patch近场声全息研究. 声学学报, 2012; 37(3): 286−293 DOI: 10.15949/j.cnki.0371-0025.2012.03.004 [21] Olivieri M, Pezzoli M, Malvermi R, et al. Near-field acoustic holography analysis with convolutional neural networks. INTER-NOISE and NOISE-CON Congress and Conference Proceedings, 2020: 5607−5618
[22] Yu L, Antoni J, Leclere Q. Spectral matrix completion by cyclic projection and application to sound source reconstruction from non-synchronous measurements. J. Sound Vib., 2016; 372: 31−49 DOI: 10.1016/j.jsv.2016.02.031
[23] Yu L, Antoni J, Leclere Q, et al. Acoustical source reconstruction from non-synchronous sequential measurements by fast iterative shrinkage thresholding algorithm. J. Sound Vib., 2017; 408: 351−367 DOI: 10.1016/j.jsv.2017.07.036
[24] Liu Q, Chu N, Yu L, et al. Efficient localization of low-frequency sound source with non-synchronous measurement at coprime positions by alternating direction method of multipliers. IEEE Trans. Instrum. Meas., 2022; 71: 1−12 DOI: 10.1109/TIM.2022.3156992
[25] Yu S, Zhang C, Hu D, et al. Sound sources reconstruction based on the non-synchronous measurements of microphone arrays and equivalent sources modeling. 5th International Conference on Information Communication and Signal Processing, IEEE, Shenzhen, China, 2022: 789−793
[26] Leclère Q. Multi-channel spectral analysis of multi-pass acquisition measurements. Mech. Syst. Signal Process., 2009; 23(5): 1415−1422 DOI: 10.1016/j.ymssp.2008.12.002
[27] 张永斌, 毕传兴, 陈剑, 等. 基于等效源法的平面近场声全息及其实验研究. 声学学报, 2007; 32(6): 489−496 DOI: 10.15949/j.cnki.0371-0025.2007.06.002 [28] 于飞, 陈心昭, 李卫兵, 等. 空间声场全息重建的波叠加方法研究. 物理学报, 2004; 53(8): 2607−2613 DOI: 10.3321/j.issn:1000-3290.2004.08.035 [29] Yoon S H, Nelson P A. A method for the efficient construction of acoustic pressure cross-spectral matrices. J. Sound Vib., 2000; 233(5): 897−920 DOI: 10.1006/jsvi.1999.2888
[30] Nelson P A, Yoon S H. Estimation of acoustic source strength by inverse methods: Part I, Conditioning of the inverse problem. J. Sound Vib., 2000; 233(4): 639−664 DOI: 10.1006/jsvi.1999.2837
计量
- 文章访问数: 59
- HTML全文浏览量: 4
- PDF下载量: 24
