CN113589232B - 基于重构正则化矩阵的奇异值分解微振动振源定位方法 - Google Patents

Abstract

本发明公开了一种基于重构正则化矩阵的奇异值分解微振动振源定位方法，在空间中放置N个加速度传感器，根据空间位置关系得到定位模型，选取其中一个传感器作为参考传感器，通过作差将非线性方程组转换为线性方程组，写成矩阵形式得到矩阵方程；令迭代索引j＝0，设置迭代门限值ε，采用最小二乘法求解矩阵方程，将方程得到的解作为迭代初始值然后对矩阵方程的系数矩阵进行奇异值分解，求正则化参数，构造正则化矩阵，并结合正则化参数计算得到正则化估计值；计算正则化估计值的均方误差，计算出均方误差取最小值时所对应的新的正则化矩阵，进而计算新的正则化估计值若则此时的为优化结果，确定环境微振动振源位置。

Description

基于重构正则化矩阵的奇异值分解微振动振源定位方法

技术领域

本发明涉及电子工业厂房环境微振动响应分析技术领域，主要涉及一种基于重构正则化矩阵的奇异值分解微振动振源定位方法。

背景说明

随着科技的不断发展，各种新型电子工业不断涌现，直拉炉、薄膜设备、光刻设备以及过程测试设备等半导体生产测试设备均对环境微振动提出严格要求，环境微振动的测量及控制成为影响半导体制造领域生产检测环节的重要因素。环境微振动属于随机振动，具有振动持续时间长、振动强度小、振动频带宽的特点，在实际的生产生活中，引起环境微振动的振源构成十分复杂，包括自然条件、临近公路与轨道交通荷载、建筑施工荷载、临近工厂机械设备振动和爆破冲击荷载等外部振源；还有厂房内部的生产设备、人员走动以及动力设备等内部振源。环境微振动经过土体传递至建筑基础或地坪，再传递至精密设备及厂房本体，影响设备工作。所以对电子工业厂房环境微振动的响应分析十分重要，对环境微振动振源进行定位，掌握其位置和分布信息，从而为厂房隔振提供依据，为精密加工提供保障，减小经济损失。

目前振源定位应用于矿难灾害应急通信、文化遗产地保护、卫星通讯等领域。传统方法有：高斯牛顿法、chan算法、粒子群算法等。TDOA定位方法的本质是求解线性方程组，但是当求解线性方程组时，若其系数矩阵的条件数较大，为病态线性方程组，传统算法没有考虑病态情况，应用于病态线性方程组求解时难以保证数值解的精度，无法解决实际问题，无法很好地应用于电子工业厂房环境微振动响应分析技术领域。

病态线性方程组的高精度求解问题被应用于电磁反演算法、遥感影像、GNSS 等领域，还没有应用于环境微振动振源定位算法中。解决病态问题最主要的方法是Tikhonov正则化，影响正则化解算的两个关键因素是正则化参数和正则化矩阵，但是目前对正则化矩阵的构造方法研究较少，并且小奇异值会对矩阵的病态造成很大影响，选取小奇异值边界的方法研究也较少，现有正则化方法的估计结果与真值偏差较大。

由于已有的振源定位算法无法很好地应用于环境微振动定位领域，因此本领域需要一种高精度求解病态线性方程组的振源定位算法。

发明内容：

本发明的目的在于提出一种基于重构正则化矩阵的奇异值分解微振动振源定位方法，解决传统定位算法无法很好地求解病态方程组的技术问题，可以将模型转化为良性问题来提高求解精度，提高解的稳定性，应用性强。

本发明是通过以下技术方案实现的：

一种基于重构正则化矩阵的奇异值分解微振动振源定位方法，包括以下步骤：

S1：设P为环境微振动振源，在空间中放置N个加速度传感器，建立三维空间直角坐标系，环境微振动振源坐标为P(x，y，z)，加速度传感器坐标为 S_i(x_i，y_i，z_i)，其中i＝1，2，…，N，根据空间位置关系得到TDOA定位模型，选取其中一个传感器作为参考传感器，通过作差将非线性方程组转换为线性方程组，写成矩阵形式得到矩阵方程；

S2：在空间中获得环境微振动振源信号到达N个传感器的加速度时程信号 x_i(t)，其中参考传感器的加速度时程信号为x₁(t)(即选取第1个传感器为参考传感器)，振源到第i个传感器的传播时间为t_i，振源到参考传感器的传播时间为 t₁，将第i个传感器的加速度时程信号x_i(t)和参考传感器的加速度时程信号x₁(t) 进行广义互相关运算，得到振源到第i个传感器与振源到参考传感器之间的时延值τ_i1＝t_i-t₁；

S3：令迭代索引j＝0，设置迭代门限值ε，采用最小二乘法求解步骤S1的矩阵方程，将方程得到的解作为迭代初始值

S4：对矩阵方程的系数矩阵进行奇异值分解；

S5：利用广义交叉验证GCV方法求关于步骤S1中矩阵方程的正则化参数；

S6：基于奇异值分解构造正则化矩阵，并结合正则化参数计算得到正则化估计值；

S7：计算步骤S6中正则化估计值的均方误差，计算出均方误差取最小值时所对应的新的正则化矩阵，进而计算新的正则化估计值

S8：若则转至步骤S9，否则更新迭代索引j，将S7得到的新的正则化估计值替换上一次循环步骤S7中的正则化估计值，去计算其均方误差并计算出均方误差取最小值时所对应的新的正则化矩阵，进而计算新的正则化估计值，循环执行步骤S7-S8；

S9：此时的为优化结果，确定环境微振动振源位置。

在上述技术方案中，步骤S1中：

TDOA定位模型为：

式中，d_i为环境微振动振源到第i个传感器的距离，d_i＝vt_i，t_i为环境微振动振源到第i个传感器所需的时间；v为环境微振动信号传播速度；

将第一个传感器作为参考传感器，作差可将非线性方程组转换为线性方程组，如式(2)所示：

式中，x_i1＝x_i-x₁，y_i1＝y_i-y₁，z_i1＝z_i-z₁，d_i1＝d_i-d₁，i＝1，2，…，N；

S1-3：将式(2)写成矩阵形式，得到矩阵方程，如式(3)所示：

AX＝b(3)

式中，

在上述技术方案中，步骤S2具体包括：

S2-1：计算第i个传感器的加速度时程信号x_i(t)和参考传感器的加速度时程信号x₁(t)的互相关函数对互相关函数/>进行加窗平滑处理，如式 (4)所示：

式中，为接收到的两组信号的互相关谱，φ(f)为窗函数，/>为广义相关谱；

S2-2：由式(4)得到时延值τ_i1，如式(5)所示：

此时，式(2)中d_i1＝d_i-d₁＝vτ_i1。

在上述技术方案中，步骤S3具体包括：

令迭代索引j＝0，设置迭代门限值ε，采用最小二乘法计算矩阵方程的解，如式(6)所示：

将式(6)得到的解作为迭代初始值如式(7)所示：

在上述技术方案中，步骤S4具体包括：

S4-1：对矩阵A进行奇异值分解可得

式中，A是m×n阶矩阵，U＝[u₁，u₂，…，u_m]和V＝[v₁，v₂，…，v_n]是正交矩阵， u_i是m维列向量，v_i是n维列向量，满足是矩阵A的左右奇异值向量，I是单位矩阵，即/>λ₁≥λ₂≥…≥λ_n≥0，λ₁，λ₂，…，λ_n为系数矩阵的奇异值；

S4-2：将式(8)代入式(7)得：

式中，W是U的前n列，即W＝[u₁，u₂，…，u_n]；

S4-3：最小二乘估计量的观测值校正为：

在上述技术方案中，步骤S5具体包括：

广义交叉验证GCV函数为：

式中，trace表示矩阵的迹，即矩阵对角线元素之和，当函数取最小值时的α即为最优正则化参数。

在上述技术方案中，步骤S6具体包括：

S6-1：构造的正则化矩阵R_t为：

式中，V_i为V的第i列；

S6-2：正则化估计值为：

在上述技术方案中，步骤S7具体包括：

S7-1：正则化估计量的方差为：

式中，

S7-2：正则化估计量的偏差的迹为：

S7-3：正则化估计量的均方误差为：

S7-4：当t取t＝1，2，…，n时，计算不同t值对应的均方误差，当均方误差取最小值时，令k＝t，即确定小奇异值边界λ_k＝λ_t，此时正则化矩阵R_k为：

S7-5：计算正则化估计量

在上述技术方案中，步骤S9确定的环境微振动振源位置：

本发明还提供一种计算机可读存储介质，存储有计算机程序，所述计算机程序被执行时实现上述的方法。

本发明的优点和有益效果为：

本发明的环境微振动振源TDOA定位方法相对于传统定位方法，考虑了线性方程组的病态问题，提高了求解精度和解的稳定性；与此同时，相对于其他奇异值分解算法，本发明的利用新的正则化矩阵构造方法，减小了小奇异值对矩阵病态的影响，提高精度。

附图说明

图1是本发明流程框图。

图2是三维TDOA定位模型。

对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，可以根据以上附图获得其他的相关附图。

具体实施方式：

为了使本技术领域的人员更好地理解本发明方案，下面结合具体实施例进一步说明本发明的技术方案。

实施例一

一种基于重构正则化矩阵的奇异值分解微振动振源定位方法，具体包括如下步骤：

S1：设P为环境微振动振源，在空间中放置N个加速度传感器，建立三维空间直角坐标系，环境微振动振源坐标为P(x，y，z)，加速度传感器坐标为 S_i(x_i，y_i，z_i)，其中i＝1，2，…，N，根据空间位置关系得到TDOA定位模型，TDOA 定位模型为：

式中，d_i为环境微振动振源到第i个传感器的距离，d_i＝vt_i，t_i为环境微振动振源到第i个传感器所需的时间；v为环境微振动信号传播速度；

将第一个传感器作为参考传感器，作差可将非线性方程组转换为线性方程组，如式(2)所示：

式中，x_i1＝x_i-x₁，y_i1＝y_i-y₁，z_i1＝z_i-z₁，d_i1＝d_i-d₁，i＝1，2，…，N；

将式(2)写成矩阵形式，得到矩阵方程，如式(3)所示：

AX＝b(3)

式中，

S2：在空间中获得环境微振动振源信号到达N个传感器的加速度时程信号 x_i(t)，其中参考传感器的加速度时程信号为x₁(t)(即选取第1个传感器为参考传感器)，振源到第i个传感器的传播时间为t_i，振源到参考传感器的传播时间为 t₁，将第i个传感器的加速度时程信号x_i(t)和参考传感器的加速度时程信号x₁(t) 进行广义互相关运算，得到振源到第i个传感器与振源到参考传感器之间的时延值τ_i1＝t_i-t₁；时延值τ_i1的具体计算步骤如下：

计算第i个传感器的加速度时程信号x_i(t)和参考传感器的加速度时程信号 x₁(t)的互相关函数对互相关函数/>进行加窗平滑处理，如式(4) 所示：

式中，为接收到的两组信号的互相关谱，φ(f)为窗函数，/>为广义相关谱；

由式(4)得到时延值τ_i1，如式(5)所示：

此时，式(2)中d_i1＝d_i-d₁＝vτ_i1。

S3：令迭代索引j＝0，设置迭代门限值ε，采用最小二乘法求解步骤S1的矩阵方程，将方程得到的解作为迭代初始值

采用最小二乘法计算矩阵方程的解，如式(6)所示：

将式(6)得到的解作为迭代初始值如式(7)所示：

S4：对矩阵方程的系数矩阵进行奇异值分解，具体的讲：

对矩阵A进行奇异值分解可得

式中，A是m×n阶矩阵，U＝[u₁，u₂，…，u_m]和V＝[v₁，v₂，…，v_n]是正交矩阵， u_i是m维列向量，v_i是n维列向量，满足是矩阵A的左右奇异值向量，即/>λ₁≥λ₂≥…≥λ_n≥0，λ₁，λ₂，…，λ_n为系数矩阵的奇异值；

将式(8)代入式(7)得：

式中，W是U的前n列，即W＝[u₁，u₂，…，u_n]；

最小二乘估计量的观测值校正为：

S5：利用广义交叉验证GCV方法求关于步骤S1中矩阵方程的正则化参数；广义交叉验证GCV函数为：

式中，trace表示矩阵的迹，即矩阵对角线元素之和，当函数取最小值时的α即为最优正则化参数。

S6：基于奇异值分解构造正则化矩阵，并结合正则化参数计算得到正则化估计值。具体的讲：

构造的正则化矩阵R_t为：

式中，V_i为V的第i列；

正则化估计值为：

S7：计算步骤S6中正则化估计值的均方误差，计算出均方误差取最小值时所对应的新的正则化矩阵，进而计算新的正则化估计值具体的讲：

正则化估计量的方差为：

式中，

正则化估计量的偏差的迹为：

正则化估计量的均方误差为：

当t取t＝1，2，…，n时，计算不同t值对应的均方误差，当均方误差取最小值时，令k＝t，即确定小奇异值边界λ_k＝λ_t，此时新的正则化矩阵R_k为：

计算新的正则化估计量

S8：若则转至步骤S9，否则更新迭代索引j，将S7得到的新的正则化估计值替换上一次循环步骤S7中的正则化估计值，去计算其均方误差并计算出均方误差取最小值时所对应的新的正则化矩阵，进而计算新的正则化估计值，循环执行步骤S7-S8。

S9：此时的为优化结果，确定环境微振动振源位置，环境微振动振源位置

实施例二

本实施应用上述方法对环境微振动振源定位精度进行验证：

本实施例中，在地表面放置5个传感器，选取其中一个传感器作为参考传感器，环境微振动振源S分别传递信号到设置好的参考传感器和其余传感器，根据接收到的信号求解时间差。

本实施例中通过测试获得振动波传播速度为v＝275m/s，建立二维直角坐标系，获得传感器坐标分别为：A(0，0)、B(4.225，0.1306)、C(2.3288，1.2978)、 D(4.1684，2.428)、E(-0.0586，2.4103)，模拟振源坐标为S(-4.1505，0)。

利用实施例一的方法，经过计算求得振源位置为(-4.2973，-0.1511)，实际振源位置为S(-4.1505，0)。本发明有效地解决了传统TDOA方案中线性方程组的病态问题，提高了求解精度和解的稳定性；与此同时，相对于其他奇异值分解算法，本发明的利用新的正则化矩阵构造方法，减小了小奇异值对矩阵病态的影响，提高精度。

Claims (9)

1.一种基于重构正则化矩阵的奇异值分解微振动振源定位方法，其特征在于，包括以下步骤：

S1：设P为环境微振动振源，在空间中放置N个加速度传感器，建立三维空间直角坐标系，环境微振动振源坐标为P(x，y，z)，加速度传感器坐标为S_i(x_i，y_i，z_i)，其中i＝1，2，…，N，根据空间位置关系得到TDOA定位模型，选取其中一个传感器作为参考传感器，通过作差将非线性方程组转换为线性方程组，写成矩阵形式得到矩阵方程；

S2：在空间中获得环境微振动振源信号到达N个传感器的加速度时程信号x_i(t)，其中参考传感器的加速度时程信号为x₁(t)(即选取第1个传感器为参考传感器)，振源到第i个传感器的传播时间为t_i，振源到参考传感器的传播时间为t₁，将第i个传感器的加速度时程信号x_i(t)和参考传感器的加速度时程信号x₁(t)进行广义互相关运算，得到振源到第i个传感器与振源到参考传感器之间的时延值τ_i1＝t_i-t₁；

S3：令迭代索引j＝0，设置迭代门限值ε，采用最小二乘法求解步骤S1的矩阵方程，将方程得到的解作为迭代初始值

S4：对矩阵方程的系数矩阵进行奇异值分解；

S5：利用广义交叉验证GCV方法求关于步骤S1中矩阵方程的正则化参数；

S6：基于奇异值分解构造正则化矩阵，并结合正则化参数计算得到正则化估计值；

S7：计算步骤S6中正则化估计值的均方误差，计算出均方误差取最小值时所对应的新的正则化矩阵，进而计算新的正则化估计值

S8：若则转至步骤S9，否则更新迭代索引j，将S7得到的新的正则化估计值替换上一次循环步骤S7中的正则化估计值，去计算其均方误差并得到新的正则化估计值，循环执行步骤S7-S8；

S9：此时的为优化结果，确定环境微振动振源位置。

2.根据权利要求1所述的基于重构正则化矩阵的奇异值分解微振动振源定位方法，其特征在于，步骤S1中：

TDOA定位模型为：

式中，d_i为环境微振动振源到第i个传感器的距离，d_i＝vt_i，t_i为环境微振动振源到第i个传感器所需的时间；v为环境微振动信号传播速度；

将第一个传感器作为参考传感器，作差可将非线性方程组转换为线性方程组，如式(2)所示：

式中，x_i1＝x_i-x₁，y_i1＝y_i-y₁，z_i1＝z_i-z₁，d_i1＝d_i-d₁，i＝1，2，…，N；

S1-3：将式(2)写成矩阵形式，得到矩阵方程，如式(3)所示：

AX＝b(3)

式中，

3.根据权利要求2所述的基于重构正则化矩阵的奇异值分解微振动振源定位方法，其特征在于，步骤S2具体包括：

S2-1：计算第i个传感器的加速度时程信号x_i(t)和参考传感器的加速度时程信号x₁(t)的互相关函数对互相关函数/>进行加窗平滑处理，如式(4)所示：

式中，为接收到的两组信号的互相关谱，φ(f)为窗函数，/>为广义相关谱；

S2-2：由式(4)得到时延值τ_i1，如式(5)所示：

此时，式(2)中d_i1＝d_i-d₁＝vτ_i1。

4.根据权利要求3所述的基于重构正则化矩阵的奇异值分解微振动振源定位方法，其特征在于，步骤S3具体包括：

采用最小二乘法计算矩阵方程的解，如式(6)所示：

将式(6)得到的解作为迭代初始值如式(7)所示：

5.根据权利要求4所述的基于重构正则化矩阵的奇异值分解微振动振源定位方法，其特征在于，步骤S4具体包括：

S4-1：对矩阵A进行奇异值分解可得

式中，A是m×n阶矩阵，U＝[u₁，u₂，…，u_m]和V＝[v₁，v₂，…，v_n]是正交矩阵，u_i是m维列向量，v_i是n维列向量，满足是矩阵A的左右奇异值向量，/>即/>λ₁≥λ₂≥…≥λ_n≥0，λ₁，λ₂，…，λ_n为系数矩阵的奇异值；

S4-2：将式(8)代入式(7)得：

式中，W是U的前n列，即W＝[u₁，u₂，…，u_n]；

S4-3：最小二乘估计量的观测值校正为：

6.根据权利要求5所述的基于重构正则化矩阵的奇异值分解微振动振源定位方法，其特征在于，步骤S5具体包括：

广义交叉验证GCV函数为：

式中，trace表示矩阵的迹，即矩阵对角线元素之和，当函数取最小值时的α即为最优正则化参数。

7.根据权利要求6所述的基于重构正则化矩阵的奇异值分解微振动振源定位方法，其特征在于，步骤S6具体包括：

S6-1：构造的正则化矩阵R_t为：

式中，V_i为V的第i列；

S6-2：正则化估计值为：

8.根据权利要求7所述的基于重构正则化矩阵的奇异值分解微振动振源定位方法，其特征在于，步骤S7具体包括：

S7-1：正则化估计量的方差为：

式中，

S7-2：正则化估计量的偏差的迹为：

S7-3：正则化估计量的均方误差为：

S7-4：当t取t＝1，2，…，n时，计算不同t值对应的均方误差，当均方误差取最小值时，令k＝t，即确定小奇异值边界λ_k＝λ_t，此时正则化矩阵R_k为：

S7-5：计算正则化估计量

9.一种计算机可读存储介质，其特征在于，存储有计算机程序，所述计算机程序被执行时实现如权利要求1至8中任一项所述的方法。