本发明涉及机械振动信号处理和减振降噪领域，特别涉及基于landweber迭代法的运行工况传递路径分析方法。
背景技术：
水下兵器的声隐身性能声隐身性能是衡量其安全性和作战能力的重要指标，因此，壳体结构装备振动噪声的有效监测与控制对于提高装备性能具有重要工程意义。振动或噪声传递路径的识别是问题的关键，目前常采用传递路径分析(tpa)或者运行工况传递路径分析(otpa)的方法来识别振动或噪声传递路径，而传统的传递路径分析方法由于其复杂的频响函数测试以及载荷识别过程导致在工程实际中难以快速有效的识别传递路径，运行工况传递路径分析方法是近年来出现的一种传递路径快速分析方法，该方法利用试验工况数据识别传递利率函数矩阵，并将其用于实际工况数据，从而得到传递路径贡献量结果，该方法简单快捷，被广泛运用于工程实际中。
现有otpa方法均存在传递率函数矩阵拟合程度低、稳定性不足等缺陷，即解反映实际情况能力、抵抗噪声能力无法满足实际工程需求。因此，对于otpa，需要采用更加有效的数学手段以确保矩阵良态，在高效的前提下提高精度。
技术实现要素：
本发明的目的在于提供基于landweber迭代法的运行工况传递路径分析方法，以解决上述问题。
为实现上述问题，本发明采用以下技术方案：
基于landweber迭代法的运行工况传递路径分析方法，包括如下步骤：
步骤1，选取所分析机械系统激励点的临近位置作为参考点，选取待分析位置为目标点；
步骤2，根据所分析机械系统设计若干种工况试验，使机械系统在每种试验工况下运行，采集每种试验工况下参考点和目标点响应信号，得到参考点响应信号矩阵和目标点响应信号矩阵，并对其进行快速傅里叶变换得到相应的频响函数矩阵x和y；
步骤3，设置迭代初值x0＝αx*，计算迭代系数矩阵xk＝xk-1+x0(i-xxk-1)，在传递率函数矩阵tk＝xky满足morozov误差原则前持续迭代系数矩阵xk，直至迭代得到的传递率函数矩阵tk满足morozov误差原则，停止迭代；
步骤4，根据所分析机械系统进行实际运行工况实验，得到一组实际运行工况数据，根据该组实际运行工况数据和迭代得到的传递率函数矩阵计算得到估计的目标点总贡献量
步骤5，针对各条传递路径对各路径贡献量进行排序后计算出各传递路径的贡献量占比，完成运行工况传递路径分析。
进一步的，步骤2、步骤4中均采用振动加速度传感器测量响应信号。
进一步的，步骤2中试验工况数量大于等于参考点数量，步骤4中实际运行工况数量大于等于参考点数量。
进一步的，步骤2中x(ω)＝[x1(ω),x2(ω),...,xr(ω)]为全部试验工况下所有参考点响应信号，r表示试验工况数量，
进一步的，步骤2中y(ω)＝[y1(ω),y2(ω),...,yr(ω)]为全部试验工况下所有目标点响应信号，r表示试验工况数量，
进一步的，步骤3中迭代初值x0＝αx*，α为松弛因子，且0<α<1/||x||2，x*表示参考点频率响应函数x的共轭转置。
进一步的，步骤3中计算迭代系数矩阵xk＝xk-1+x0(i-xxk-1)，i为单位矩阵，k表示迭代次数。
进一步的，步骤3中morozov误差原则为：当迭代得到的传递率函数矩阵满足||x(ω)tk(ω)-y(ω)||≤τδ时停止迭代，常数τ＝2，误差水平δ＝10-3。
进一步的，步骤4中
与现有技术相比，本发明有以下技术效果：
本发明根据所分析机械系统设计试验工况，同时采集试验工况下参考点响应信号和目标点响应信号，针对传递率函数矩阵求解的不适定反问题特点，利用landweber迭代法对传递率函数矩阵的求解进行迭代控制，设置morozov误差原则为迭代终止准则，极大地提高了传递率函数矩阵求解的稳定性的同时有效地提高了其拟合程度，从而提高其抵抗噪声干扰能力和反映实际情况能力。然后采集一组实际运行工况下待分析机械系统的参考点响应信号和目标点响应信号，将迭代得到的传递率函数矩阵与实际工况下测得的参考点响应信号相乘，得到传递路径贡献量识别结果，对不同路径贡献量进行排序，得到各个传递路径的贡献量占比，完成运行工况传递路径分析。本发明克服了以往运行工况传递路径分析方法中传递率函数矩阵求解的拟合程度低和稳定性不足带来的识别误差问题，提高了传递路径贡献量识别精度。
附图说明
图1为本发明优选实例中所述的球状辐射声源仿真系统；
图2为本发明优选实例中所述的声源a源贡献量；
图3为本发明优选实例中所述的声源b源贡献量；
图4为本发明优选实例中所述的目标点总贡献量；
图5为本发明优选实例中所述的各贡献量识别精度柱状图；
图6为本发明优选实例中所述的声源a对目标点的源贡献量相对误差；
图7为本发明优选实例中所述的声源b对目标点的源贡献量相对误差；
图8为本发明优选实例中所述的特征频率处目标点总贡献量相对误差；
图9为本发明优选实例中所述的特征频率处目标点总贡献量相对误差；
图10为本发明优选实例中所述的求解稳定性比较柱状图；
具体实施方式
下面结合具体的实施例对本发明做进一步的详细说明，所述是对本发明的解释而不是限定。
基于landweber迭代法的运行工况传递路径分析方法，包括如下步骤：
步骤1，选取所分析机械系统激励点的临近位置作为参考点，选取待分析位置为目标点；
步骤2，根据所分析机械系统设计若干种工况试验，使机械系统在每种试验工况下运行，采集每种试验工况下参考点和目标点响应信号，得到参考点响应信号矩阵和目标点响应信号矩阵，并对其进行快速傅里叶变换得到相应的频响函数矩阵x和y；
步骤3，设置迭代初值x0＝αx*，计算迭代系数矩阵xk＝xk-1+x0(i-xxk-1)，在传递率函数矩阵tk＝xky满足morozov误差原则前持续迭代系数矩阵xk，直至迭代得到的传递率函数矩阵tk满足morozov误差原则，停止迭代；
步骤4，根据所分析机械系统进行实际运行工况实验，得到一组实际运行工况数据，根据该组实际运行工况数据和迭代得到的传递率函数矩阵计算得到估计的目标点总贡献量
步骤5，针对各条传递路径对各路径贡献量进行排序后计算出各传递路径的贡献量占比，完成运行工况传递路径分析。
步骤2、步骤4中均采用振动加速度传感器测量响应信号。
步骤2中试验工况数量大于等于参考点数量，步骤4中实际运行工况数量大于等于参考点数量。
步骤2中x(ω)＝[x1(ω),x2(ω),...,xr(ω)]为全部试验工况下所有参考点响应信号，r表示试验工况数量，
步骤2中y(ω)＝[y1(ω),y2(ω),...,yr(ω)]为全部试验工况下所有目标点响应信号，r表示试验工况数量，
步骤3中迭代初值x0＝αx*，α为松弛因子，且0<α<1/||x||2，x*表示参考点频率响应函数x的共轭转置。
步骤3中计算迭代系数矩阵xk＝xk-1+x0(i-xxk-1)，i为单位矩阵，k表示迭代次数。
步骤3中morozov误差原则为：当迭代得到的传递率函数矩阵满足||x(ω)tk(ω)-y(ω)||≤τδ时停止迭代，常数τ＝2，误差水平δ＝10-3。
步骤4中
本发明根据所分析机械系统设计试验工况，同时采集试验工况下参考点响应信号和目标点响应信号，针对传递率函数矩阵求解的不适定反问题特点，利用landweber迭代法对传递率函数矩阵的求解进行迭代控制，设置morozov误差原则为迭代终止准则，极大地提高了传递率函数矩阵求解的稳定性的同时有效地提高了其拟合程度，从而提高其抵抗噪声干扰能力和反映实际情况能力。本发明在设计试验工况过程中，选取所分析机械系统临近激励点位置作为参考点，待分析位置作为目标点，参考点个数与激励源个数相同，并且试验工况数量大于等于参考点数量。将迭代得到的传递率函数矩阵和实际工况下测得的参考点响应信号相乘，得到传递路径贡献量识别结果，对不同路径贡献量进行排序，得到各个传递路径的贡献量占比，完成运行工况传递路径分析。本发明克服了以往运行工况传递路径分析方法中传递率函数矩阵求解的拟合程度低和稳定性不足带来的识别误差问题，提高了传递路径贡献量识别精度。
优选实例，为说明基于landweber迭代法的otpa相对于以往otpa方法的优越性，对比基于截断奇异值分解的传统otpa和基于tikhonov正则化的otpa方法，设计球状辐射声源仿真系统，参见图1所示。
设声源1与声源2的半径均为a＝b＝0.34m；
声源表面振动速度均为v1＝v2＝60m/s，且同相位；
声源1和声源2到参考点距离均为r1ref＝r2ref＝0.85m；
声源1和声源2到目标点距离均为r1＝r2＝10.25；
声源1到声源2的距离r12＝13.6m。
参考点和目标点处的声压是声源1和声源2贡献量的向量和，可以表达为：
式中：noise(ω)——随机噪声。
设置分析频率为0～100hz，频率分辨率为δf＝1hz；
声源1和声源2分别在25hz和50hz输出信号；
随机噪声在0～100hz输出信号。
截断奇异值分解(tsvd)参数设置：
截断奇异值分解方法的关键是选择合适的截断系数k。截断系数k取值过大，会导致传递率函数矩阵严重病态(不稳定)；截断系数过小，会导致截断误差过大，传递率函数矩阵的拟合程度不足。
传统otpa根据奇异值贡献率：
来确定截断系数k，在本次仿真过程中取奇异值贡献率为3％。
tikhonov正则化参数设置：
tikhonov正则化参数λ的取值非常关键，λ＞0且λ∈(σn,σ1)。
最常用的正则化参数选择方法有广义交叉检验方法和l曲线方法，本次仿真选用l曲线法选取正则化参数，即求解[log(||axλ-b||)，log(||xλ||)]曲率极大值确定最优正则化参数λ。
landweber迭代法终止准则设置：
如前所述，利用landweber迭代法求解反问题的关键是要寻找一个合适的终止准则使迭代近似解在半收敛附近停止迭代。常用的迭代准则是morozov误差原则：
为使迭代近似解在半收敛附近停止迭代，设置常数τ为τ＝2，误差水平δ取δ＝10-3。
仿真步骤设定：
1)基于10组试验工况数据，利用上述三种方法分别求解传递率函数矩阵；
2)再仿真1组实际工况下的数据作为校验标准，包括单独开启声源a时目标点的响应、单独开启声源b时目标点的响应以及同时开启两个声源时目标点的响应；
3)利用步骤1)中求解的传递率函数矩阵和实际工况下的参考点响应分别求出声源a对目标点的源贡献量、声源b对目标点的源贡献量以及目标点的总贡献量；
4)对比步骤2)和步骤3)的结果差异，并分别计算出相对误差；
5)按整数倍增大噪声倍数至10倍，并分别对不同噪声倍数重复上述过程。
图2-图4所示为所选取的随机噪声增大至6倍时三种方法的贡献量对比图，其中practical代表实际工况下测得的目标点响应，图2所示为声源a源贡献量对比图，图3所示为声源b源贡献量对比图，图4所示为目标点总贡献量。
图5所示为随机噪声增大至6倍时三种方法的各贡献量识别精度的精确比较柱状图，表明了本方法的贡献量识别精度明显好于另外两种对比方法。
图6-图9所示为逐步增大噪声倍数至10倍时源贡献量相对误差和目标点总贡献量相对误差变化趋势折线图，图6所示为声源a对目标点的源贡献量相对误差对比折线图，图7所示为声源b对目标点的源贡献量相对误差对比折线图，图8所示为目标点总贡献量在特征频率25hz处相对误差对比折线图，图9所示为目标点总贡献量在特征频率50hz处相对误差对比折线图，表明了本方法的对源贡献量的识别精度和总贡献量的识别精度在整体上要明显优于另外两种对比方法，显示识别的传递率函数矩阵具有更高的拟合程度。
为了评价本方法识别传递率函数矩阵的稳定性，分别计算三种方法的各组相对误差相邻数值间的方差并求和，即：
图10所示为三种方法的贡献量求解稳定性的对比图，表明本方法计算得到的各贡献量识别相对误差相邻数值间的方差和远小于另外两种对比方法，即本方法的求解稳定性具有明显优势，显示了本方法传递率函数矩阵的抵抗噪声干扰能力明显优于另外两种对比方法，在实际工程应用中具有更好的应用价值。