CN103559340A - 一种基于comi-pso算法的不相关多源频域载荷识别方法 - Google Patents

Abstract

本发明一种基于COMI-PSO算法的不相关多源频域载荷识别方法，首先根据多个载荷源不相关的性质，对载荷识别的动力学方程进行了解耦；其次，利用基于一元线性回归模型和最小二乘解的传递函数进行识别；再次，寻找一个可行载荷解，使得所测得的响应的最大相对误差最小；最后，采用COMI-PSO算法来搜索求解该单目标最优化问题，识别多个不相关载荷源；本发明根据多个测点的响应信号，能同时识别多个不相关频域载荷源，能彻底解决传递函数获取困难，在共振频率处矩阵求逆所出现的病态问题，通过本发明识别的多元载荷激励的精度和对测量噪声的敏感程度均优于传统方法，基本能满足3db的工程精度要求。

Description

一种基于COMI-PSO算法的不相关多源频域载荷识别方法

技术领域

本发明涉及一种基于COMI-PSO算法的不相关多源频域载荷识别方法。

背景技术

载荷识别是通过测量结构动态响应和系统特征来求结构所受激励的方法，属于振动问题中的第二类反问题。激励源是引起系统结构产生振动的主要原因，但在实际工程中，如导弹在空中飞行、火车在轨道上行驶、海洋平台和桥梁等大型建筑物受风浪等随机激励作用等情况下，很难对作用于结构的外载荷进行直接测量或计算，甚至有时因载荷作用点不可到达，使这种动态载荷不可测。载荷识别技术为那些无法直接测量载荷的结构或者系统提供了一种识别动态载荷的有效方法。而准确地确定载荷、科学地制定相应的载荷谱是可靠性试验、力学控制、铁路交通和桥梁设计等重大工程设计时面临的迫切问题。

载荷识别作为一个动力学反问题，存在不适定性。从结构响应数据中反求载荷是目前研究的热点和难点，其存在性、唯一性和确定性都缺乏严格的理论证明。载荷识别主要分为频域法和时域法两类。其中频域法提出较早，主要利用激励和响应之间的频响函数求逆来实现，但是矩阵求逆法在应用时通常需要求解广义逆，且经常会遇到系数矩阵的病态问题和奇异值分解问题。为克服频响函数求逆时在固有频率附近秩缺或病态的问题，Ojalvo和张令弥等采用了小量分解法，后张德文发展了改进小量分解法，李东升等又进一步提出了广义小量分解法。张丽霞等将神经网络应用于频域载荷识别，胡迪科、毛文涛等将支持向量机应用于频域载荷识别，胡杰等人利用优化技术将载荷识别的反问题转化为正问题来处理。段瑞玲、王慧儒等根据正问题的激励和响应数据建立逆系统的ARMA模型，从而将载荷识别问题变为参数辨识问题。

神经网络、支持向量机、优化技术、参数辨识等方法参数较多，物理意义和概念不明确，且在应用上比较复杂。当系统和激励复杂时，寻找最优激励的过程非常长，效率很低，且存在过拟合和过学习的情况，精度也无法保证。

发明内容

本发明的目的在于提出一种基于COMI-PSO算法的不相关多源频域载荷识别方法，克服了传统载荷识别方法中出现的阵求逆载传递函数获取困难、直接求逆会出现矩阵奇异性和病态、寻优过程复杂、时间效率低、识别精度无法保证等问题，有效提高载荷识别精度。

一种基于COMI-PSO算法的多源频域载荷识别方法，适用于响应测点的个数大于等于激励的个数、线性不变的系统、平稳随机的载荷激励、各个载荷源近似不相关或独立，包括如下步骤：

步骤1、设系统有m个载荷激励输入f_i(i=1,…,m)，n个测点输出y_j(j=1,…,n)，对应于每一个输出y_j，有m个脉冲响应函数h_ji(t),i=1,…,m，j=1,…n,测得多输入多输出系统在m个不相关载荷同时作用下n个响应测点的时域响应y_j(t)j=1,2,…,n，计算响应测点间的互功率谱密度矩阵S_yy(ω)，其中ω为频率；

步骤2、对系统依次分立对各个载荷点施加单个激励输入f_i(i=1,…,m)，计算其自功率谱为

测得系统在该单独激励输入下各响应测点的响应输出并计算其自功率谱

按照(15)式的一元线性回归模型识别出的传递函数模的平方

$\left\{\begin{array}{c}{\widehat{\mathrm{\β}}}_0={\stackrel{\‾}{S}}_{{\mathrm{yy}}_{j,i}}\left(\mathrm{\ω}\right)-{\left|{\hat{H}}_{j,i}\left(\mathrm{\ω}\right)\right|}^2{\stackrel{\‾}{S}}_{f_i}\left(\mathrm{\ω}\right)\\ {\left|{\hat{H}}_{j,i}\left(\mathrm{\ω}\right)\right|}^2=\frac{\underset{l=1}{\overset{k}{\mathrm{\Σ}}}[S_{{\mathrm{yy}}_{j,i}^l}\left(\mathrm{\ω}\right)-{\stackrel{\‾}{S}}_{{\mathrm{yy}}_{j,i}}\left(\mathrm{\ω}\right)][S_{f_i^l}\left(\mathrm{\ω}\right)-{\stackrel{\‾}{S}}_{f_i}\left(\mathrm{\ω}\right)]}{\underset{l=1}{\overset{k}{\mathrm{\Σ}}}{[S_{f_i^l}\left(\mathrm{\ω}\right)-{\stackrel{\‾}{S}}_{f_i}\left(\mathrm{\ω}\right)]}^2}\end{array}\right.---\left(15\right)$

(15)式中，为载荷点i所施加的k次载荷的平均值，为响应点j在载荷点i所施加的k次响应的平均值；

步骤3、利用载荷源的不相关特性对载荷识别方程(1)式进行解耦得到（4）式：

$\begin{array}{c}{S}_{\mathrm{yy}}\left(\mathrm{\ω}\right)=\frac{1}{2\mathrm{\π}}{\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}h\left(u\right)C_{\mathrm{ff}}(\mathrm{\τ}+u-v)\·h^T\left(v\right)e^{-\mathrm{j\ω\τ}}\mathrm{dudvu\τ}\\ =\stackrel{\‾}{H}\left(\mathrm{\ω}\right)S_{\mathrm{ff}}\left(\mathrm{\ω}\right)H^T\left(\mathrm{\ω}\right)\end{array}---\left(1\right)$

(1)式中h(u)是系统的单位脉冲响应矩阵，h^T(u)是系统的单位脉冲响应矩阵的转置，

是输入各激励之间的协方差函数矩阵，是系统频率特性矩阵，H^T(ω)是系统频率特性矩阵的转置，是系统频率特性矩阵的共轭，S_ff(ω)是输入各激励之间的互功率谱密度，

在m个输入载荷激励都是零均值的平稳随机过程，且互不相关的情况下，输入协方差函数矩阵

为对角阵，即：输入功率谱矩阵S_ff(ω)也为对角阵此时，输出功率谱中主对角线上的任意一元素 $S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\ω}\right)(j=\mathrm{1,2},\·\·\·,n)$ 满足：

$\begin{array}{c}{S}_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\ω}\right)=\left[\begin{array}{ccccc}{\stackrel{\‾}{H}}_{j1}\left(\mathrm{\ω}\right)& \·\·\·& {\stackrel{\‾}{H}}_{\mathrm{ji}}\left(\mathrm{\ω}\right)& \·\·\·& {\stackrel{\‾}{H}}_{\mathrm{jn}}\end{array}\left(\mathrm{\ω}\right)\right]\·\mathrm{diag}\left[S_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\ω}\right)\right]\·{\left[\begin{array}{ccccc}{H}_{j1}\left(\mathrm{\ω}\right)& \·\·\·& H_{\mathrm{ji}}\left(\mathrm{\ω}\right)& \·\·\·& H_{\mathrm{jn}}\left(\mathrm{\ω}\right)\end{array}\right]}^T\\ =\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\stackrel{\‾}{H}}_{\mathrm{ji}}\left(\mathrm{\ω}\right)S_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\ω}\right)H_{\mathrm{ji}}^T\left(\mathrm{\ω}\right)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\left|H_{\mathrm{ji}}\left(\mathrm{\ω}\right)\right|}^2{S}_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\ω}\right)\end{array}---\left(3\right)$

(3)式写成矩阵后的形式为：

其中，|H_j,i(ω)|²是输入f_i对响应y_j的传递函数模的平方，

是待识别的载荷源f_i的自功率谱，

是响应y_j的自功率谱；

记 ${\stackrel{\→}{S}}_Y\left(\mathrm{\ω}\right)\stackrel{\mathrm{\Δ}}{=}{\left[\begin{array}{ccc}{S}_{{\mathrm{yy}}_{11}}\left(\mathrm{\ω}\right)& \·\·\·& S_{{\mathrm{yy}}_{\mathrm{nn}}}\left(\mathrm{\ω}\right)\end{array}\right]}^T$ ，

则(4)式可简写为： ${\stackrel{\→}{S}}_Y\left(\mathrm{\ω}\right)=B\left(\mathrm{\ω}\right){\stackrel{\→}{S}}_F\left(\mathrm{\ω}\right);$

步骤4、对(4)式识别出的S_yy(ω)和

采用(16)式传统最小二乘广义逆计算(4)式载荷值：

${\stackrel{\→}{S}}_F^{\′}\left(\mathrm{\ω}\right)={[B{\left(\mathrm{\ω}\right)}^{T}B\left(\mathrm{\ω}\right)]}^{-1}B{\left(\mathrm{\ω}\right)}^T{\stackrel{\→}{S}}_Y\left(\mathrm{\ω}\right)---\left(16\right)$

步骤5、利用(27)式条件数的定义计算（4）式的条件数，若条件数小于阈值，直接使用（16）式的载荷值作为最终的识别载荷值,载荷识别结束；如果条件数大于阈值，则使用（16）式的计算结果控制粒子群的初始化位置，并随机初始化粒子的初试速度；

在求解方程Ax=b的过程中，引入的舍入误差总会导致其数值解x或多或少地不等于其理论解x'=A^-1B，

定义误差为： $\mathrm{\δx}\stackrel{\mathrm{\Δ}}{=}x-x^{\′}---\left(25\right)$

定义剩余为： $\mathrm{\δb}\stackrel{\mathrm{\Δ}}{=}\mathrm{Ax}-b---\left(26\right)$

定义条件数为： $\mathrm{cond}\left(A\right)\stackrel{\mathrm{\Δ}}{=}\left|\right|A^{-1}\left|\right|\left|\right|A\left|\right|---\left(27\right);$

步骤6、以(20)式作为适应值函数，计算每个粒子的适应值：

$\min f\left(S_{\mathrm{ff}}^{\′}\left(\mathrm{\ω}\right)\right)=\min \left(\underset{j}{\max }\right\{\frac{|S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\ω}\right)-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\left|H_{j,i}\left(\mathrm{\ω}\right)\right|}^2{S}_{{\mathrm{ff}}_{\mathrm{ii}}^{\′}}\left(\mathrm{\ω}\right)|}{\left|S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\ω}\right)\right|}\left\}\right)---\left(20\right)$

步骤7、找到每个粒子的历史最优位置和最优适应值：

对每个粒子而言，如果它的当前适应度比其历史最优适应度还要小，则使用当前的适应度替代该粒子的最优适应度，并保存当前位置为该粒子历史最优位置；

步骤8、找到整个群体的最优位置和最优适应值；

对每个粒子而言，如果它的当前适应度比整个群体最优适应度还要小，则使用当前的适应度替代粒子群最优适应度，并保存当前位置为整个群体最优位置；

步骤9、根据COMI-PSO算法，使用(21)式、(22)式计算各个粒子最新速度和最新位置

$v_{\mathrm{id}}^{t+1}=w*v_{\mathrm{id}}^t+c_1{r}_1(p_{\mathrm{id}}-x_{\mathrm{id}}^t)+c_2{r}_2(p_{\mathrm{gd}}-x_{\mathrm{id}}^t)---\left(21\right)$

$x_{\mathrm{id}}^{t+1}=x_{\mathrm{id}}^t+v_{\mathrm{id}}^{t+1}---\left(22\right)$

其中惯性参数w=max_W-(max_W-min_W)*ln(1+(e-1)*CT/TS)，max_W为惯性参数的上限，min_W为惯性参数的下限，e为自然对数的底数，上式能保证w的取值范围是[min_W,max_W]；r₁,r₂为(0,1)之间的随机数，p_id表示当前群体中第i个个体历史最优位置在第d维上数值，p_gd表示当前群体最优位置在第d维上的数值，CT为当前的迭代次数，TS为预设的总迭代次数，加速因子c₁,c₂分别为：

c₁=4*(w-min_W)²/(max_W-min_W)²

c₂=4*(max_W-w)²/(max_W-min_W)²

通过上述操作完成迭代过程中第i个个体从第t代的速度和位置

更新为第(t+1)代新的速度

和位置

步骤10、通过(20)式计算新位置下，每个粒子的适应度，根据遗传算法中的选择算子，保留种群中适应度较好的部分粒子直接进入下一次迭代，对另一部分适应度较小的粒子采用遗传算法中的交叉算子和突变算子进行预处理，计算预处理之后的子代的适应值，将适应度比父代好的部分相同数目的粒子取代原粒子群中的父辈进入下一次迭代；

步骤11、根据是否达到最大迭代次数或达到最好的适应值判断迭代终止条件，如果达到转入步骤12，否则转入步骤6；

步骤12、输出全局最优解，即群体最优位置，作为最终识别的多个不相关载荷源值,载荷识别算法结束。

本发明一种基于COMI-PSO算法的不相关多源频域载荷识别方法，首先根据多个载荷源不相关的性质，对载荷识别的动力学方程进行了解耦，从而简化了方程，减少了传递函数获取的数量；其次，为了消除测量噪声和系统弱非线性的影响，使获得的传递函数更为准确，提出了一种基于一元线性回归模型和最小二乘解的传递函数识别方法；再次，在深入了解载荷识别的工作原理的基础上建立相应的数学模型，将不相关多源载荷识别问题转化为一个单目标最优化正问题，目标是寻找一个可行载荷解，使得所测得的传感器响应的最大相对误差最小；然后，针对标准粒子群优化算法(Particle Swarm optimization,PSO)在寻优过程中出现的早熟、容易陷入局部最优等问题，提出了一种在粒子初始化、参数选择以及粒子群迭代三方面进行控制的综合改进粒子群优化算法(Comprehensive Improved Particle Swarm optimization,COMI-PSO)；最后，采用COMI-PSO算法来搜索求解该单目标最优化问题，识别多个不相关载荷源。

本发明根据多个测点的响应信号，能同时识别多个不相关频域载荷源，能彻底解决传递函数获取困难，在共振频率处矩阵求逆所出现的病态问题，通过本发明识别的多元载荷激励的精度和对测量噪声的敏感程度均优于传统方法，基本能满足3db的工程精度要求。

附图说明

图1为发明多输入多输出系统的示意图；

图2为本发明利用一元线性回归模型进行拟合识别传递函数示意图；

图3为COMI-PSO算法流程示意图；

图4为本发明悬臂圆柱薄壳Nastran不相关多源激励模型；

图5为本发明载荷点1到6个响应测点的传递函数；

图6为本发明载荷点2到6个响应测点的传递函数；

图7为本发明联合施加的2个不相关平稳随机集中力载荷；

图8为本发明6个响应输出点在两个集中力载荷联合作用下的功率谱响应示意图；

图9为广义逆方法识别载荷点1和载荷点2结果与真实载荷的误差比较。

图10为COMI-PSO方法识别载荷点1和载荷点2结果与真实载荷的误差比较。

以下结合具体实施例和附图对本发明作进一步详述。

具体实施方式

如图1所示，设系统有m个载荷激励输入f_i(i=1,…,m)，n个测点输出y_j(j=1,…,n)，根据叠加原理，线性系统的每一个输出都可以由各个分立输入所引起的响应叠加而成，因此，对应于每一个输出y_j，有m个脉冲响应函数h_ji(t),i=1,…,m，j=1,…,n,计算响应测点间的互功率谱密度矩阵S_yy(ω)，其中ω为频率。

$\begin{array}{c}{S}_{\mathrm{yy}}\left(\mathrm{\ω}\right)=\frac{1}{2\mathrm{\π}}{\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}h\left(u\right)C_{\mathrm{ff}}(\mathrm{\τ}+u-v)\·h^T\left(v\right)e^{-\mathrm{j\ω\τ}}\mathrm{dudvu\τ}\\ =\stackrel{\‾}{H}\left(\mathrm{\ω}\right)S_{\mathrm{ff}}\left(\mathrm{\ω}\right)H^T\left(\mathrm{\ω}\right)\end{array}---\left(1\right)$

(1)式中h(u)是系统的单位脉冲响应矩阵，h^T(u)是系统的单位脉冲响应矩阵的转置，是输入各激励之间的协方差函数矩阵，

是系统频率特性矩阵，H^T(ω)是系统频率特性矩阵的转置，

是系统频率特性矩阵的共轭，S_ff(ω)是输入各激励之间的互功率谱密度；

(1)式给出了多输入/多输出情形下输出功率谱矩阵与输入功率谱矩阵之间的关系式，它显示了输入与输出功率谱关系的简明特点，正是频域分析法的优点所在。

在实际情况中，m与n不相等，因此要求取载荷谱矩阵，须对频响函数矩阵求广义逆，则在频域中的载荷识别公式可表示为:

$S_{\mathrm{ff}}\left(\mathrm{\ω}\right)={\left[H^T\left(\mathrm{\ω}\right)\stackrel{\‾}{H}\left(\mathrm{\ω}\right)\right]}^{-1}{H}^T\left(\mathrm{\ω}\right)S_{\mathrm{yy}}\left(\mathrm{\ω}\right)\stackrel{\‾}{H}\left(\mathrm{\ω}\right){\left[H^T\left(\mathrm{\ω}\right)\stackrel{\‾}{H}\left(\mathrm{\ω}\right)\right]}^{-1}---\left(2\right)$

(2)式中

代表H^T(ω)矩阵的逆，(2)式的主要问题是用试验获得系统的复频响应函数矩阵H(ω)的工作量太大，而用有限元法来获得H(ω)又存在仿真建模与试验的误差问题。

在m个输入载荷激励都是零均值的平稳随机过程，且互不相关的情况下，输入协方差函数矩阵

为对角阵，即：

输入功率谱矩阵S_ff(ω)也为对角阵

此时，输出功率谱中主对角线上的任意一元素 $S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\ω}\right)(j=\mathrm{1,2},\·\·\·,n)$ 满足：

$\begin{array}{c}{S}_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\ω}\right)=\left[\begin{array}{ccccc}{\stackrel{\‾}{H}}_{j1}\left(\mathrm{\ω}\right)& \·\·\·& {\stackrel{\‾}{H}}_{\mathrm{ji}}\left(\mathrm{\ω}\right)& \·\·\·& {\stackrel{\‾}{H}}_{\mathrm{jn}}\end{array}\left(\mathrm{\ω}\right)\right]\·\mathrm{diag}\left[S_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\ω}\right)\right]\·{\left[\begin{array}{ccccc}{H}_{j1}\left(\mathrm{\ω}\right)& \·\·\·& H_{\mathrm{ji}}\left(\mathrm{\ω}\right)& \·\·\·& H_{\mathrm{jn}}\left(\mathrm{\ω}\right)\end{array}\right]}^T\\ =\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\stackrel{\‾}{H}}_{\mathrm{ji}}\left(\mathrm{\ω}\right)S_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\ω}\right)H_{\mathrm{ji}}^T\left(\mathrm{\ω}\right)=\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\left|H_{\mathrm{ji}}\left(\mathrm{\ω}\right)\right|}^2{S}_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\ω}\right)\end{array}---\left(3\right)$

(3)式写成矩阵后的形式为：

其中，|H_j,i(ω)|²是输入f_i对响应y_j的传递函数模的平方，

是待识别的载荷源f_i的自功率谱，

是响应y_j的自功率谱；

记 ${\stackrel{\→}{S}}_Y\left(\mathrm{\ω}\right)\stackrel{\mathrm{\Δ}}{=}{\left[\begin{array}{ccc}{S}_{{\mathrm{yy}}_{11}}\left(\mathrm{\ω}\right)& \·\·\·& S_{{\mathrm{yy}}_{\mathrm{nn}}}\left(\mathrm{\ω}\right)\end{array}\right]}^T$ ，

则(4)式可简写为： ${\stackrel{\→}{S}}_Y\left(\mathrm{\ω}\right)=B\left(\mathrm{\ω}\right){\stackrel{\→}{S}}_F\left(\mathrm{\ω}\right)$

1)当n<m，(4)式为欠定方程，对应的满足(4)式的解有无穷组；

2)当n=m，(4)式为正定方程，对应的满足(4)式的解唯一；

3)当n>m，(4)式为过定方程，无对应的满足(4)式的解。

为保证反演出载荷激励的精度，(4)式中应满足n>m，并将该问题转化为一个优化问题，目标是找一组载荷激励

使得系统的n个测点的响应能达到

为验证该方法的正确性和精度，识别出来的激励可以与实际加载的激励

进行比较。但是该问题本身是一个多目标优化问题，需要转化成单目标优化问题进行求解。

在具体的频率下，基于一元线性回归模型获取系统载荷作用点到各响应测点的传递函数：对系统依次分立施加单个输入f_i，计算其自功率谱为

测得系统在该激励f_i下的输出

并计算其自功率谱

则输入f_i对输出

的传递函数H_j,i(ω)满足：

${\left|H_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2=\frac{S_{{\mathrm{yy}}_{\mathrm{jj}}^i}\left(\mathrm{\&omega;}\right)}{S_{f_i}\left(\mathrm{\&omega;}\right)}---\left(5\right)$

如果系统是线性的，且不存在测量噪声，由(5)式，单个输入f_i的自功率谱

和输出

的自功率谱

之间存在比例关系，其比值即传递函数模的平方|H_j,i(ω)|²。但实验中载荷源和响应测点都存在测量噪声，使得在不同量级不同波形激励下每次识别出的传递函数略有不同。

设系统真实的频率特性为H(ω)，令F与Y表示系统真正的输入与输出，F′与Y′表示测量得到的输入与输出，n₁与n₂分别表示输入与输出的测量噪声，假设n₁与n₂是统计独立的零均值平稳过程，且测量噪声与系统真正的输入F（或输出Y）都是不相关或统计独立的，记输入信噪比为：

$S_F\left(\mathrm{\&omega;}\right)/S_{n_1}\left(\mathrm{\&omega;}\right)=1/{\mathrm{\&alpha;}}_1\left(\mathrm{\&omega;}\right)---\left(6\right)$

输出信噪比为：

$S_Y\left(\mathrm{\&omega;}\right)/S_{n_2}\left(\mathrm{\&omega;}\right)=1/{\mathrm{\&alpha;}}_2\left(\mathrm{\&omega;}\right)---\left(7\right)$

则在载荷和响应均含噪声的情况下，使用不同的方法求得的传递函数H₀(ω)、H₁(ω)、H₂(ω)与真实的传递函数H(ω)的关系为：

$\begin{array}{c}{\left|H_0\left(\mathrm{\&omega;}\right)\right|}^2=S_{Y^{\&prime;}{Y}^{\&prime;}}\left(\mathrm{\&omega;}\right)/S_{F^{\&prime;}{F}^{\&prime;}}\left(\mathrm{\&omega;}\right)\\ =\frac{1+{\mathrm{\&alpha;}}_2\left(\mathrm{\&omega;}\right)}{1+{\mathrm{\&alpha;}}_1\left(\mathrm{\&omega;}\right)}\&CenterDot;\frac{S_{\mathrm{YY}}\left(\mathrm{\&omega;}\right)}{S_{\mathrm{FF}}\left(\mathrm{\&omega;}\right)}=\frac{1+{\mathrm{\&alpha;}}_2\left(\mathrm{\&omega;}\right)}{1+{\mathrm{\&alpha;}}_1\left(\mathrm{\&omega;}\right)}\&CenterDot;{\left|H\left(\mathrm{\&omega;}\right)\right|}^2\end{array}---\left(8\right)$

$H_1\left(\mathrm{\&omega;}\right)=\frac{S_{F^{\&prime;}{Y}^{\&prime;}}\left(\mathrm{\&omega;}\right)}{S_{F^{\&prime;}{F}^{\&prime;}}\left(\mathrm{\&omega;}\right)}=H\left(\mathrm{\&omega;}\right)\&CenterDot;\frac{1}{1+{\mathrm{\&alpha;}}_1\left(\mathrm{\&omega;}\right)}---\left(9\right)$

$H_2\left(\mathrm{\&omega;}\right)=\frac{S_{Y^{\&prime;}{Y}^{\&prime;}}\left(\mathrm{\&omega;}\right)}{S_{Y^{\&prime;}{F}^{\&prime;}}\left(\mathrm{\&omega;}\right)}=H\left(\mathrm{\&omega;}\right)\&CenterDot;(1+{\mathrm{\&alpha;}}_2\left(\mathrm{\&omega;}\right))---\left(10\right)$

H₀(ω)、H₁(ω)、H₂(ω)与H(ω)的大小关系为：

|H₁(ω)|≤|H(ω)|≤|H₂(ω)|    (11)

|H₁(ω)|≤|H₀(ω)|≤|H₂(ω)|    (12)

而|H₀(ω)|与|H(ω)|的大小关系取决于输入信噪比和输出信噪比，无法事先确定。

载荷和响应的谱相干函数定义为：

$\begin{array}{c}{\mathrm{\&gamma;}}_{F^{\&prime;}{Y}^{\&prime;}}^2\left(\mathrm{\&omega;}\right)=\frac{{\left|S_{F^{\&prime;}{Y}^{\&prime;}}\left(\mathrm{\&omega;}\right)\right|}^2}{S_{F^{\&prime;}}\left(\mathrm{\&omega;}\right)S_{Y^{\&prime;}}\left(\mathrm{\&omega;}\right)}\\ =\frac{{\left|S_{\mathrm{FY}}\left(\mathrm{\&omega;}\right)\right|}^2}{S_F\left(\mathrm{\&omega;}\right)S_Y\left(\mathrm{\&omega;}\right)(1+{\mathrm{\&alpha;}}_1)(1+{\mathrm{\&alpha;}}_2)}=\frac{1}{(1+{\mathrm{\&alpha;}}_1)(1+{\mathrm{\&alpha;}}_2)}\&le;1\end{array}---\left(13\right)$

由(13)式，无论是输入测量噪声，还是输出测量噪声，都将使输入输出谱相干函数小于1，而且会使传递函数和载荷识别带来误差。

为了消除测量噪声和系统弱非线性的影响，使识别出的传递函数更准确，可对载荷点i多次分别施加不同谱形、不同量级的k次单源输入

其自功率谱分别为

分别测得系统的响应点j的单点输出响应

其自功率谱分别为然后，利用一元回归模型(18)式对这些点进行线性拟合，如图2所示，使用最小二乘法求解(14)式，得到的斜率即传递函数模的平方|H_j,i(ω)|²的一元线性回归模型的最小二乘估计

如式(15)所示。

$S_{{\mathrm{yy}}_{j,i}^l}\left(\mathrm{\&omega;}\right)={\mathrm{\&beta;}}_0+{\left|H_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2{S}_{f_i^l}\left(\mathrm{\&omega;}\right)+{\mathrm{\&epsiv;}}_l{\mathrm{\&epsiv;}}_l~\mathrm{NID}(0,{\mathrm{\&sigma;}}^2)---\left(14\right)$

$\left\{\begin{array}{c}{\widehat{\mathrm{\&beta;}}}_0={\stackrel{\&OverBar;}{S}}_{{\mathrm{yy}}_{j,i}}\left(\mathrm{\&omega;}\right)-{\left|{\hat{H}}_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2{\stackrel{\&OverBar;}{S}}_{f_i}\left(\mathrm{\&omega;}\right)\\ {\left|{\hat{H}}_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2=\frac{\underset{l=1}{\overset{k}{\mathrm{\&Sigma;}}}[S_{{\mathrm{yy}}_{j,i}^l}\left(\mathrm{\&omega;}\right)-{\stackrel{\&OverBar;}{S}}_{{\mathrm{yy}}_{j,i}}\left(\mathrm{\&omega;}\right)][S_{f_i^l}\left(\mathrm{\&omega;}\right)-{\stackrel{\&OverBar;}{S}}_{f_i}\left(\mathrm{\&omega;}\right)]}{\underset{l=1}{\overset{k}{\mathrm{\&Sigma;}}}{[S_{f_i^l}\left(\mathrm{\&omega;}\right)-{\stackrel{\&OverBar;}{S}}_{f_i}\left(\mathrm{\&omega;}\right)]}^2}\end{array}\right.---\left(15\right)$

(14)式中ε_l～NID(0,σ²)，l=1,2,…,k为均值为0，方差为σ的高斯随机白噪声，β₀为系统测量偏差项，应该很小，接近0。

(15)式中，为载荷点i所施加的k次载荷的平均值，

为响应点j在载荷点i所施加的k次响应的平均值。

为验证一元线性回归模型及其最小二乘解消除测量噪声和系统弱非线性影响的效果，将(15)式识别出的传递函数模的平方

与有限元仿真计算或传递实验法的结果进行比较。

利用(27)式计算(4)式的条件数，如果条件数小于阈值，直接利用响应测点间的互功率谱密度矩阵S_yy(ω)，采用一元线性回归模型最小二乘解获得的传递函数模的平方

利用(16)式传统最小二乘广义逆的方法，计算得到的载荷值作为最终的识别载荷值,载荷识别结束；如果条件数大于阈值，则将(4)式变成单目标最优化搜索的正问题，寻找一个可行载荷解 $S_{\mathrm{ff}}^{\&prime;}\left(\mathrm{\&omega;}\right)={\left[\begin{array}{ccccc}{S}_{{\mathrm{ff}}_{11}}^{\&prime;}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& S_{{\mathrm{ff}}_{\mathrm{ii}}}^{\&prime;}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& S_{{\mathrm{ff}}_{\mathrm{mm}}}^{\&prime;}\left(\mathrm{\&omega;}\right)\end{array}\right]}^T,$ 使得所测得的响应值的最大相对误差最小，最后，采用COMI-PSO算法来搜索求解(20)式的单目标最优化问题，识别多个不相关载荷源。

现在考虑方程Ax=b。在求解该方程的过程中，引入的舍入误差总会导致其数值解x或多或少的不等于其理论解x'=A^-1B。

定义误差为： $\mathrm{\&delta;x}\stackrel{\mathrm{\&Delta;}}{=}x-x^{\&prime;}---\left(25\right)$

定义剩余为： $\mathrm{\&delta;b}\stackrel{\mathrm{\&Delta;}}{=}\mathrm{Ax}-b---\left(26\right)$

在现实背景下建立的线性方程组系统中，系数矩阵和右端项往往不是精确给出的，以载荷识别为例，左端的传递函数和右端的响应必然会存在一定的误差，计算过程中也会存在一些舍入误差。因此，在实际工程中需要考虑在线性方程组系数加入一个数值扰动，计算结果改变幅度的问题。

定义条件数如下：

$\mathrm{cond}\left(A\right)\stackrel{\mathrm{\&Delta;}}{=}\left|\right|A^{-1}\left|\right|\left|\right|A\left|\right|---\left(27\right)$

如果m≠0,则有：

$\frac{\left|\right|\mathrm{\&delta;x}\left|\right|}{\left|\right|x\left|\right|}\&le;\mathrm{cond}\left(A\right)\frac{\left|\right|\mathrm{\&delta;b}\left|\right|}{\left|\right|b\left|\right|}---\left(28\right)$

解的相对误差是右端项b的相对误差的cond(A)倍，如果条件数很大，则解的误差将成倍增长。条件数很大的矩阵称之为病态矩阵，病态矩阵对应的方程组为病态方程组；反之，则称之为良态矩阵。

由以上分析可知，在条件数比较大的情况下，不相关多源频域载荷识别中直接采用矩阵求逆的方式求解会产生较大误差，无法保证精度。由此，本发明将不相关多源载荷识别问题转化为单目标的优化问题，进一步求解。

定理1：对于(4)式，当n≥m时，在响应误差平方和最小的单目标优化准则下的解S′_F(ω)为：

${\stackrel{\&RightArrow;}{S}}_F^{\&prime;}\left(\mathrm{\&omega;}\right)={[B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right)]}^{-1}B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)---\left(16\right)$

证明：在(4)式中，响应误差平方和的一半

为：

$\begin{array}{c}\frac{1}{2}{(B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))}^T(B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))\\ =\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(b_i\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-S_{{\mathrm{yy}}_{\mathrm{ii}}}\left(\mathrm{\&omega;}\right))}^2\stackrel{\mathrm{\&Delta;}}{=}J\left({\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)\right)\end{array}---\left(17\right)$

为了使J最小化，以

为参数，求J的梯度，可以得到(18)式：

$\begin{array}{c}{\&dtri;}_{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}J\left({\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)\right)={\&dtri;}_{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}\left(\frac{1}{2}{(B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))}^T(B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))\right)\\ =\frac{1}{2}{\&dtri;}_{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}({\stackrel{\&RightArrow;}{S}}_F{\left(\mathrm{\&omega;}\right)}^{T}B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}^{T}B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)+{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))\\ =\frac{1}{2}{\&dtri;}_{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}\mathrm{tr}({\stackrel{\&RightArrow;}{S}}_F{\left(\mathrm{\&omega;}\right)}^T{B\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^T{B\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)+{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))\\ =\frac{1}{2}{\&dtri;}_{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}(\mathrm{tr}{\stackrel{\&RightArrow;}{S}}_F{\left(\mathrm{\&omega;}\right)}^{T}B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-2\mathrm{tr}{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right))\\ =\frac{1}{2}({B\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)+B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-2B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))\\ ={B\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)\end{array}---\left(18\right)$

为了使J最小化，使(18)式最后结果等于零，从而得到以下等式：

$B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)=B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)---\left(19\right)$

对(19)式等式左右两边同时左乘B(ω)^TB(ω)矩阵的逆[B(ω)^TB(ω)]^-1，最后得到(16)式的结果，证毕。

(16)式又叫(4)式的Moore-Penrose逆或最小二乘广义逆。

在上述不相关多源频率载荷识别的数学模型中，为避免传统的矩阵求逆方式求解过程存在病态问题，本发明将不相关多源载荷识别问题转化为一个单目标最优化搜索的正问题，寻找一个可行载荷解 $S_{\mathrm{ff}}^{\&prime;}\left(\mathrm{\&omega;}\right)={\left[\begin{array}{ccccc}{S}_{{\mathrm{ff}}_{11}}^{\&prime;}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& S_{{\mathrm{ff}}_{\mathrm{ii}}}^{\&prime;}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& S_{{\mathrm{ff}}_{\mathrm{mm}}}^{\&prime;}\left(\mathrm{\&omega;}\right)\end{array}\right]}^T,$ 使得所测得的响应值的最大相对误差最小：

$\min f\left(S_{\mathrm{ff}}^{\&prime;}\left(\mathrm{\&omega;}\right)\right)=\min \left(\underset{j}{\max }\right\{\frac{|S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\&omega;}\right)-\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left|H_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2{S}_{{\mathrm{ff}}_{\mathrm{ii}}^{\&prime;}}\left(\mathrm{\&omega;}\right)|}{\left|S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\&omega;}\right)\right|}\left\}\right)---\left(20\right)$

(20)式的物理意义是：求得识别的载荷S′_ff(ω)后，代入（4）式，计算得到的响应会和真实的响应

之间有一定的误差，目标函数(20)式是使n个响应点中的最大相对误差达到最小。

为验证该方法的正确性和精度，将识别出来的激励S′_ff(ω)与实际加载的激励进行比较。

粒子群优化算法(PSO)是解决上述优化问题的常用方法之一，但考虑到标准的PSO算法存在早熟、容易陷入局部最优等问题，本发明采用改进后的COMI-PSO算法处理上述最优化问题。COMI-PSO算法是在标准PSO算法基础上对粒子初始化、参数选择以及粒子群迭代进行综合控制，主要从以下三个方面提高PSO算法的寻优效率：

(I）在标准PSO算法中，对粒子群个体进行随机初始化，存在一定的盲目性，为有效缩短了粒子群收敛时间并综合传统识别结果，COMI-PSO采用最小二乘法计算得到的多源载荷识别结果控制粒子群个体初始化；

(II）为了有效控制PSO算法中粒子群局部搜索和全局搜索，进一步提高算法搜索时效，COMP-PSO对标准PSO算法中的惯性参数、加速因子进行优化，其中惯性参数采用对数非线性递减策略，加速因子采用非线性渐变策略；

（III）在标准PSO算法处理过程中，随着种群多样性的降低，种群中大部分个体集中在一个较小的范围内，致使整个群体失去大范围搜索能力，算法陷入早熟。为使种群多样性不致丧失过早或者使其始终保持在一定水平，COMI-PSO算法在粒子群整体迭代之前，采用GA算法(遗传算法)对粒子群进行预处理，包括使用选择算子选取适应度较小的部分粒子以及使用交叉算子和变异算子处理选择的粒子。

(IV)使用(4)式的计算结果控制粒子群初始化位置，并随机初始化粒子的初试速度；

(V)以(20)式作为适应度函数，计算每个粒子的适应度；

(VI)对每个粒子而言，如果它的当前适应度比其历史最优适应度还要小，则使用当前的适应度替代该粒子的最优适应度，并保存当前位置为该粒子历史最优位置；

(VII)对每个粒子而言，如果它的当前适应度比整个群体最优适应度还要小，则使用当前的适应度替代粒子群最优适应度，并保存当前位置为整个群体最优位置；

(VIII)根据COMI-PSO算法，使用(21)式、(22)式计算各个粒子最新速度

和最新位置

$v_{\mathrm{id}}^{t+1}=w*v_{\mathrm{id}}^t+c_1{r}_1(p_{\mathrm{id}}-x_{\mathrm{id}}^t)+c_2{r}_2(p_{\mathrm{gd}}-x_{\mathrm{id}}^t)---\left(21\right)$

$x_{\mathrm{id}}^{t+1}=x_{\mathrm{id}}^t+v_{\mathrm{id}}^{t+1}---\left(22\right)$

其中惯性参数w=max_W-(max_W-min_W)*ln(1+(e-1)*CT/TS)，max_W为惯性参数的上限，min_W为惯性参数的下限，e为自然对数的底数，上式能保证w的取值范围是[min_W,max_W]；r₁,r₂为(0,1)之间的随机数，p_id表示当前群体中第i个个体历史最优位置在第d维上数值，p_gd表示当前群体最优位置在第d维上的数值，CT为当前的迭代次数，TS为预设的总迭代次数，加速因子c₁,c₂分别为：

c₁=4*(w-min_W)²/(max_W-min_W)²     (23)

c₂=4*(max_W-w)²/(max_W-min_W)²     (24)

通过上述操作完成迭代过程中第i个个体从第t代的速度

和位置

更新为第(t+1)代新的速度

和位置

(Ⅸ)通过(20)式计算新位置下，每个粒子的适应度，根据遗传算法中的选择算子，保留种群中适应度较好的部分粒子直接进入下一次迭代，对另一部分适应度较小的粒子采用遗传算法中的交叉算子和突变算子进行预处理，计算预处理之后的子代的适应值，将适应度比父代好的部分相同数目的粒子取代原粒子群中的父辈进入下一次迭代；

(X)根据是否达到最大迭代次数或达到最好的适应值判断迭代终止条件，如果达到转入步骤(XI)，否则转入步骤(V)；

(XI)输出全局最优解，即群体最优位置，从而识别出多个不相关载荷源，算法结束。

实施例：以悬臂圆柱薄壳的有限元模型为对象，利用Hypermesh建模，划分网格，设定材料，施加2个不相关的平稳随机集中力载荷，添加6个响应输出点的传感器，其位置如图4所示。选用有限元仿真计算软件Nastran，计算出2个载荷点到6个响应点的传递函数模型，如图5和6所示，采用有限元功率谱线性叠加理论，同时施加载荷谱如图7所示的2个不相关的平稳随机集中力载荷作用，计算得到6个响应输出点在2个集中力载荷联合作用下的响应，如图8所示。鉴于模型单元数较小，激励的施加频率和响应测点的计算频率范围在10～500Hz。以图5和6所示的传递函数和图8所示的响应数据，分别利用(16)式使用传统的最小二乘广义逆方式和本专利提到的COMI-PSO方式求解，并将识别的载荷和图7所联合施加的2个真实载荷进行比较，采用传统的最小二乘广义逆方式对联合施加的两个不相关的平稳随机集中力载荷进行识别的结果如图9所示，采用COMI-PSO算法识对联合施加的两个不相关的平稳随机集中力载荷进行识别的结果如图10所示。

如图3所示，一种基于COMI-PSO算法的多源频域载荷识别方法，适用于响应测点的个数大于等于激励的个数、线性不变的系统、平稳随机的载荷激励、各个载荷源近似不相关或独立，具体包括如下步骤：

步骤1、测得图1所示的多输入多输出系统在m个不相关载荷同时作用下n个响应测点的时域响应y_j(t)j=1,2,…,n，计算响应测点间的互功率谱密度矩阵S_yy(ω)；

步骤2、对系统依次分立对各个载荷点施加单个输入f_i(i=1,…,m)，计算其自功率谱为

测得系统在该单独激励下各响应测点的响应输出

并计算其自功率谱

按照(15)式的一元线性回归模型识别出的传递函数模的平方

对系统依次分立施加单个输入f_i，计算其自功率谱

测得系统在该激励f_i下的输出

并计算其自功率谱

则输入f_i对输出的传递函数H_j,i(ω)满足：

${\left|H_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2=\frac{S_{{\mathrm{yy}}_{\mathrm{jj}}^i}\left(\mathrm{\&omega;}\right)}{S_{f_i}\left(\mathrm{\&omega;}\right)}---\left(5\right)$

如果系统是线性的，且不存在测量噪声，由(5)式，单个输入f_i的自功率谱

和输出

的自功率谱

之间存在比例关系，其比值即传递函数模的平方|H_j,i(ω)|²。但实验中载荷源和响应测点都存在测量噪声，使得在不同量级不同波形激励下每次识别出的传递函数略有不同。

设系统真实的频率特性为H(ω)，令F与Y表示系统真正的输入与输出，F′与Y′表示测量得到的输入与输出，n₁与n₂分别表示输入与输出的测量噪声，假设n₁与n₂是统计独立的零均值平稳过程，且测量噪声与系统真正的输入F（或输出Y）都是不相关或统计独立的，记输入信噪比为：

$S_F\left(\mathrm{\&omega;}\right)/S_{n_1}\left(\mathrm{\&omega;}\right)=1/{\mathrm{\&alpha;}}_1\left(\mathrm{\&omega;}\right)---\left(6\right)$

输出信噪比为：

$S_Y\left(\mathrm{\&omega;}\right)/S_{n_2}\left(\mathrm{\&omega;}\right)=1/{\mathrm{\&alpha;}}_2\left(\mathrm{\&omega;}\right)---\left(7\right)$

则在载荷和响应均含噪声的情况下，使用不同的方法求得的传递函数H₀(ω)、H₁(ω)、H₂(ω)与真实的传递函数H(ω)的关系为：

$\begin{array}{c}{\left|H_0\left(\mathrm{\&omega;}\right)\right|}^2=S_{Y^{\&prime;}{Y}^{\&prime;}}\left(\mathrm{\&omega;}\right)/S_{F^{\&prime;}{F}^{\&prime;}}\left(\mathrm{\&omega;}\right)\\ =\frac{1+{\mathrm{\&alpha;}}_2\left(\mathrm{\&omega;}\right)}{1+{\mathrm{\&alpha;}}_1\left(\mathrm{\&omega;}\right)}\&CenterDot;\frac{S_{\mathrm{YY}}\left(\mathrm{\&omega;}\right)}{S_{\mathrm{FF}}\left(\mathrm{\&omega;}\right)}=\frac{1+{\mathrm{\&alpha;}}_2\left(\mathrm{\&omega;}\right)}{1+{\mathrm{\&alpha;}}_1\left(\mathrm{\&omega;}\right)}\&CenterDot;{\left|H\left(\mathrm{\&omega;}\right)\right|}^2\end{array}---\left(8\right)$

$H_1\left(\mathrm{\&omega;}\right)=\frac{S_{F^{\&prime;}{Y}^{\&prime;}}\left(\mathrm{\&omega;}\right)}{S_{F^{\&prime;}{F}^{\&prime;}}\left(\mathrm{\&omega;}\right)}=H\left(\mathrm{\&omega;}\right)\&CenterDot;\frac{1}{1+{\mathrm{\&alpha;}}_1\left(\mathrm{\&omega;}\right)}---\left(9\right)$

$H_2\left(\mathrm{\&omega;}\right)=\frac{S_{Y^{\&prime;}{Y}^{\&prime;}}\left(\mathrm{\&omega;}\right)}{S_{Y^{\&prime;}{F}^{\&prime;}}\left(\mathrm{\&omega;}\right)}=H\left(\mathrm{\&omega;}\right)\&CenterDot;(1+{\mathrm{\&alpha;}}_2\left(\mathrm{\&omega;}\right))---\left(10\right)$

H₀(ω)、H₁(ω)、H₂(ω)与H(ω)的大小关系为：

|H₁(ω)|≤|H(ω)|≤|H₂(ω)|    (11)

|H₁(ω)|≤|H₀(ω)|≤|H₂(ω)|    (12)

而|H₀(ω)|与|H(ω)|的大小关系取决于输入信噪比和输出信噪比，无法事先确定。

载荷和响应的谱相干函数定义为：

$\begin{array}{c}{\mathrm{\&gamma;}}_{F^{\&prime;}{Y}^{\&prime;}}^2\left(\mathrm{\&omega;}\right)=\frac{{\left|S_{F^{\&prime;}{Y}^{\&prime;}}\left(\mathrm{\&omega;}\right)\right|}^2}{S_{F^{\&prime;}}\left(\mathrm{\&omega;}\right)S_{Y^{\&prime;}}\left(\mathrm{\&omega;}\right)}\\ =\frac{{\left|S_{\mathrm{FY}}\left(\mathrm{\&omega;}\right)\right|}^2}{S_F\left(\mathrm{\&omega;}\right)S_Y\left(\mathrm{\&omega;}\right)(1+{\mathrm{\&alpha;}}_1)(1+{\mathrm{\&alpha;}}_2)}=\frac{1}{(1+{\mathrm{\&alpha;}}_1)(1+{\mathrm{\&alpha;}}_2)}\&le;1\end{array}---\left(13\right)$

由(13)式，无论是输入测量噪声，还是输出测量噪声，都将使输入输出谱相干函数小于1，而且会使传递函数和载荷识别带来误差。

为了消除测量噪声和系统弱非线性的影响，使识别出的传递函数更准确，可对载荷点i多次分别施加不同谱形、不同量级的k次单源输入

其自功率谱分别为

分别测得系统的响应点j的单点输出响应

其自功率谱分别为然后，利用一元回归模型(18)式对这些点进行线性拟合，如图2所示，使用最小二乘法求解(14)式，得到的斜率即传递函数模的平方|H_j,i(ω)|²的一元线性回归模型的最小二乘估计如式(15)所示。

$S_{{\mathrm{yy}}_{j,i}^l}\left(\mathrm{\&omega;}\right)={\mathrm{\&beta;}}_0+{\left|H_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2{S}_{f_i^l}\left(\mathrm{\&omega;}\right)+{\mathrm{\&epsiv;}}_l{\mathrm{\&epsiv;}}_l~\mathrm{NID}(0,{\mathrm{\&sigma;}}^2)---\left(14\right)$

$\left\{\begin{array}{c}{\widehat{\mathrm{\&beta;}}}_0={\stackrel{\&OverBar;}{S}}_{{\mathrm{yy}}_{j,i}}\left(\mathrm{\&omega;}\right)-{\left|{\hat{H}}_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2{\stackrel{\&OverBar;}{S}}_{f_i}\left(\mathrm{\&omega;}\right)\\ {\left|{\hat{H}}_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2=\frac{\underset{l=1}{\overset{k}{\mathrm{\&Sigma;}}}[S_{{\mathrm{yy}}_{j,i}^l}\left(\mathrm{\&omega;}\right)-{\stackrel{\&OverBar;}{S}}_{{\mathrm{yy}}_{j,i}}\left(\mathrm{\&omega;}\right)][S_{f_i^l}\left(\mathrm{\&omega;}\right)-{\stackrel{\&OverBar;}{S}}_{f_i}\left(\mathrm{\&omega;}\right)]}{\underset{l=1}{\overset{k}{\mathrm{\&Sigma;}}}{[S_{f_i^l}\left(\mathrm{\&omega;}\right)-{\stackrel{\&OverBar;}{S}}_{f_i}\left(\mathrm{\&omega;}\right)]}^2}\end{array}\right.---\left(15\right)$

(14)式中ε_l～NID(0,σ²)，l=1,2,…,k为均值为0，方差为σ的高斯随机白噪声，β₀为系统测量偏差项，应该很小，接近0。

(15)式中，

为载荷点i所施加的k次载荷的平均值，

为响应点j在载荷点i所施加的k次响应的平均值。

为验证一元线性回归模型及其最小二乘解消除测量噪声和系统弱非线性影响的效果，将(15)式识别出的传递函数模的平方

与有限元仿真计算或传递实验法的结果进行比较。

步骤3、利用载荷源的不相关特性对载荷识别方程(1)式进行解耦得到（4）式：

设系统有m个载荷激励输入f_i(i=1,…,m)，n个测点输出y_j(j=1,…,n)，对应于每一个输出y_j，有m个脉冲响应函数h_ji(t),i=1,…,m，j=1,…n,计算响应测点间的互功率谱密度矩阵S_yy(ω)，其中ω为频率。

$\begin{array}{c}{S}_{\mathrm{yy}}\left(\mathrm{\&omega;}\right)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}h\left(u\right)C_{\mathrm{ff}}(\mathrm{\&tau;}+u-v)\&CenterDot;h^T\left(v\right)e^{-\mathrm{j\&omega;\&tau;}}\mathrm{dudvu\&tau;}\\ =\stackrel{\&OverBar;}{H}\left(\mathrm{\&omega;}\right)S_{\mathrm{ff}}\left(\mathrm{\&omega;}\right)H^T\left(\mathrm{\&omega;}\right)\end{array}---\left(1\right)$

(1)式中h(u)是系统的单位脉冲响应矩阵，h^T(u)是系统的单位脉冲响应矩阵的转置，

是输入各激励之间的协方差函数矩阵，

是系统频率特性矩阵，H^T(ω)是系统频率特性矩阵的转置，

是系统频率特性矩阵的共轭，S_ff(ω)是输入各激励之间的互功率谱密度；

在m个输入载荷激励都是零均值的平稳随机过程，且互不相关的情况下，输入协方差函数矩阵

为对角阵，即：输入功率谱矩阵S_ff(ω)也为对角阵此时，输出功率谱中主对角线上的任意一元素 $S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\&omega;}\right)(j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n)$ 满足：

$\begin{array}{c}{S}_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\&omega;}\right)=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{H}}_{j1}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& {\stackrel{\&OverBar;}{H}}_{\mathrm{ji}}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& {\stackrel{\&OverBar;}{H}}_{\mathrm{jn}}\end{array}\left(\mathrm{\&omega;}\right)\right]\&CenterDot;\mathrm{diag}\left[S_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\&omega;}\right)\right]\&CenterDot;{\left[\begin{array}{ccccc}{H}_{j1}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& H_{\mathrm{ji}}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& H_{\mathrm{jn}}\left(\mathrm{\&omega;}\right)\end{array}\right]}^T\\ =\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\stackrel{\&OverBar;}{H}}_{\mathrm{ji}}\left(\mathrm{\&omega;}\right)S_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\&omega;}\right)H_{\mathrm{ji}}^T\left(\mathrm{\&omega;}\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left|H_{\mathrm{ji}}\left(\mathrm{\&omega;}\right)\right|}^2{S}_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\&omega;}\right)\end{array}---\left(3\right)$

(3)式写成矩阵后的形式为：

其中，|H_j,i(ω)|²是输入f_i对响应y_j的传递函数模的平方，

是待识别的载荷源f_i的自功率谱，

是响应y_j的自功率谱；

记 ${\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)\stackrel{\mathrm{\&Delta;}}{=}{\left[\begin{array}{ccc}{S}_{{\mathrm{yy}}_{11}}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& S_{{\mathrm{yy}}_{\mathrm{nn}}}\left(\mathrm{\&omega;}\right)\end{array}\right]}^T$ ，

则(4)式可简写为： ${\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)=B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)$

1)当n<m，(4)式为欠定方程，对应的满足(4)式的解有无穷组；

2)当n=m，(4)式为正定方程，对应的满足(4)式的解唯一；

3)当n>m，(4)式为过定方程，无对应的满足(4)式的解。

为保证反演出载荷激励的精度，(4)式中应满足n>m，并将该问题转化为一个优化问题，目标是找一组载荷激励

使得系统的n个测点的响应能达到为验证该方法的正确性和精度，识别出来的激励

可以与实际加载的激励

进行比较。但是该问题本身是一个多目标优化问题，需要转化成单目标优化问题进行求解。

步骤4、采用(16)式传统最小二乘广义逆计算(4)式载荷值：定理1：对于(4)式，当n≥m时，在响应误差平方和最小的单目标优化准则下的解S′_F(ω)为：

${\stackrel{\&RightArrow;}{S}}_F^{\&prime;}\left(\mathrm{\&omega;}\right)={[B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right)]}^{-1}B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)---\left(16\right)$

证明：在(4)式中，响应误差平方和的一半

为：

$\begin{array}{c}\frac{1}{2}{(B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))}^T(B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))\\ =\frac{1}{2}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(b_i\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-S_{{\mathrm{yy}}_{\mathrm{ii}}}\left(\mathrm{\&omega;}\right))}^2\stackrel{\mathrm{\&Delta;}}{=}J\left({\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)\right)\end{array}---\left(17\right)$

为了使J最小化，以

为参数，求J的梯度，可以得到(18)式：

$\begin{array}{c}{\&dtri;}_{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}J\left({\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)\right)={\&dtri;}_{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}\left(\frac{1}{2}{(B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))}^T(B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))\right)\\ =\frac{1}{2}{\&dtri;}_{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}({\stackrel{\&RightArrow;}{S}}_F{\left(\mathrm{\&omega;}\right)}^{T}B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^{T}B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)+{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))\\ =\frac{1}{2}{\&dtri;}_{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}\mathrm{tr}({\stackrel{\&RightArrow;}{S}}_F{\left(\mathrm{\&omega;}\right)}^T{B\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^T{B\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)-{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)+{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))\\ =\frac{1}{2}{\&dtri;}_{{\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)}(\mathrm{tr}{\stackrel{\&RightArrow;}{S}}_F{\left(\mathrm{\&omega;}\right)}^{T}B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-2\mathrm{tr}{\stackrel{\&RightArrow;}{S}}_Y{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right))\\ =\frac{1}{2}({B\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)+B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-2B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right))\\ ={B\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)-B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)\end{array}---\left(18\right)$

为了使J最小化，使(18)式最后结果等于零，从而得到以下等式：

$B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right)=B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)---\left(19\right)$

对(19)式等式左右两边同时左乘B(ω)^TB(ω)矩阵的逆[B(ω)^TB(ω)]^-1，最后得到(16)式的结果，证毕。

(16)式又叫(4)式的Moore-Penrose逆或最小二乘广义逆。

步骤5、利用(27)式计算（4）式的条件数，若条件数小于阈值，直接使用（16）式的载荷值作为最终的识别载荷值,载荷识别结束；如果条件数大于阈值，则使用（16）式的计算结果控制粒子群的初始化位置，并随机初始化粒子的初试速度；

现在考虑方程Ax=b。在求解该方程的过程中，引入的舍入误差总会导致其数值解x或多或少的不等于其理论解x'=A^-1B。

定义误差为： $\mathrm{\&delta;x}\stackrel{\mathrm{\&Delta;}}{=}x-x^{\&prime;}---\left(25\right)$

定义剩余为： $\mathrm{\&delta;b}\stackrel{\mathrm{\&Delta;}}{=}\mathrm{Ax}-b---\left(26\right)$

在现实背景下建立的线性方程组系统中，系数矩阵和右端项往往不是精确给出的，以载荷识别为例，左端的传递函数和右端的响应必然会存在一定的误差，计算过程中也会存在一些舍入误差。因此，在实际工程中需要考虑在线性方程组系数加入一个数值扰动，计算结果改变幅度的问题。

定义条件数如下：

$\mathrm{cond}\left(A\right)\stackrel{\mathrm{\&Delta;}}{=}\left|\right|A^{-1}\left|\right|\left|\right|A\left|\right|---\left(27\right)$

如果m≠0,则有：

$\frac{\left|\right|\mathrm{\&delta;x}\left|\right|}{\left|\right|x\left|\right|}\&le;\mathrm{cond}\left(A\right)\frac{\left|\right|\mathrm{\&delta;b}\left|\right|}{\left|\right|b\left|\right|}---\left(28\right)$

解的相对误差是右端项b的相对误差的cond(A)倍，如果条件数很大，则解的误差将成倍增长。条件数很大的矩阵称之为病态矩阵，病态矩阵对应的方程组为病态方程组；反之，则称之为良态矩阵。

由以上分析可知，在条件数比较大的情况下，不相关多源频域载荷识别中直接采用矩阵求逆的方式求解会产生较大误差，无法保证精度。

步骤6、以(20)式作为适应值函数，计算每个粒子的适应值；

在上述不相关多源频率载荷识别的数学模型中，为避免传统的矩阵求逆方式求解过程存在病态问题，本发明将不相关多源载荷识别问题转化为一个单目标最优化搜索的正问题，寻找一个可行载荷解 $S_{\mathrm{ff}}^{\&prime;}\left(\mathrm{\&omega;}\right)={\left[\begin{array}{ccccc}{S}_{{\mathrm{ff}}_{11}}^{\&prime;}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& S_{{\mathrm{ff}}_{\mathrm{ii}}}^{\&prime;}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& S_{{\mathrm{ff}}_{\mathrm{mm}}}^{\&prime;}\left(\mathrm{\&omega;}\right)\end{array}\right]}^T,$ 使得所测得的响应值的最大相对误差最小：

$\min f\left(S_{\mathrm{ff}}^{\&prime;}\left(\mathrm{\&omega;}\right)\right)=\min \left(\underset{j}{\max }\right\{\frac{|S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\&omega;}\right)-\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left|H_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2{S}_{{\mathrm{ff}}_{\mathrm{ii}}^{\&prime;}}\left(\mathrm{\&omega;}\right)|}{\left|S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\&omega;}\right)\right|}\left\}\right)---\left(20\right)$

(20)式的物理意义是：求得识别的载荷S′_ff(ω)后，代入（4）式，计算得到的响应

会和真实的响应之间有一定的误差，目标函数(20)式是使n个响应点中的最大相对误差达到最小。

为验证该方法的正确性和精度，将识别出来的激励S′_ff(ω)与实际加载的激励

进行比较；

步骤7、找到每个粒子的历史最优位置和最优适应值；

对每个粒子而言，如果它的当前适应度比其历史最优适应度还要小，则使用当前的适应度替代该粒子的最优适应度，并保存当前位置为该粒子历史最优位置；

步骤8、找到整个群体的最优位置和最优适应值；

对每个粒子而言，如果它的当前适应度比整个群体最优适应度还要小，则使用当前的适应度替代粒子群最优适应度，并保存当前位置为整个群体最优位置；

步骤9、根据COMI-PSO算法，使用(21)式、(22)式计算各个粒子最新速度

和最新位置

$v_{\mathrm{id}}^{t+1}=w*v_{\mathrm{id}}^t+c_1{r}_1(p_{\mathrm{id}}-x_{\mathrm{id}}^t)+c_2{r}_2(p_{\mathrm{gd}}-x_{\mathrm{id}}^t)---\left(21\right)$

$x_{\mathrm{id}}^{t+1}=x_{\mathrm{id}}^t+v_{\mathrm{id}}^{t+1}---\left(22\right)$

其中惯性参数w=max_W-(max_W-min_W)*ln(1+(e-1)*CT/TS)，max_W为惯性参数的上限，min_W为惯性参数的下限，e为自然对数的底数，上式能保证w的取值范围是[min_W,max_W]；r₁,r₂为(0,1)之间的随机数，p_id表示当前群体中第i个个体历史最优位置在第d维上数值，p_gd表示当前群体最优位置在第d维上的数值，CT为当前的迭代次数，TS为预设的总迭代次数，加速因子c₁,c₂分别为：

c₁=4*(w-min_W)²/(max_W-min_W)²

c₂=4*(max_W-w)²/(max_W-min_W)²

通过上述操作完成迭代过程中第i个个体从第t代的速度和位置

更新为第(t+1)代新的速度

和位置

步骤10、通过(20)式计算新位置下，每个粒子的适应度，根据遗传算法中的选择算子，保留种群中适应度较好的部分粒子直接进入下一次迭代，对另一部分适应度较小的粒子采用遗传算法中的交叉算子和突变算子进行预处理，计算预处理之后的子代的适应值，将适应度比父代好的部分相同数目的粒子取代原粒子群中的父辈进入下一次迭代；

步骤11、根据是否达到最大迭代次数或达到最好的适应值判断迭代终止条件，如果达到转入步骤12，否则转入步骤6；

步骤12、输出全局最优解，即群体最优位置，作为最终识别的多个不相关载荷源值,载荷识别算法结束。

以上所述，仅是本发明较佳实施例而已，并非对本发明的技术范围作任何限制，故凡是依据本发明的技术实质对以上实施例所作的任何细微修改、等同变化与修饰，均仍属于本发明技术方案的范围内。

Claims (1)

1.一种基于COMI-PSO算法的多源频域载荷识别方法，适用于响应测点的个数大于等于激励的个数、线性不变的系统、平稳随机的载荷激励、各个载荷源近似不相关或独立，其特征在于包括如下步骤：

步骤1、设系统有m个载荷激励输入f_i(i=1,…,m)，n个测点输出y_j(j=1,…,n)，对应于每一个输出y_j，有m个脉冲响应函数h_ji(t),i=1,…,m，j=1,…n,测得多输入多输出系统在m个不相关载荷同时作用下n个响应测点的时域响应y_j(t)j=1,2,…,n，计算响应测点间的互功率谱密度矩阵S_yy(ω)，其中ω为频率；

步骤2、对系统依次分立对各个载荷点施加单个激励输入f_i(i=1,…,m)，计算其自功率谱为

测得系统在该单独激励输入下各响应测点的响应输出

并计算其自功率谱

按照(15)式的一元线性回归模型识别出的传递函数模的平方

$\left\{\begin{array}{c}{\widehat{\mathrm{\&beta;}}}_0={\stackrel{\&OverBar;}{S}}_{{\mathrm{yy}}_{j,i}}\left(\mathrm{\&omega;}\right)-{\left|{\hat{H}}_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2{\stackrel{\&OverBar;}{S}}_{f_i}\left(\mathrm{\&omega;}\right)\\ {\left|{\hat{H}}_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2=\frac{\underset{l=1}{\overset{k}{\mathrm{\&Sigma;}}}[S_{{\mathrm{yy}}_{j,i}^l}\left(\mathrm{\&omega;}\right)-{\stackrel{\&OverBar;}{S}}_{{\mathrm{yy}}_{j,i}}\left(\mathrm{\&omega;}\right)][S_{f_i^l}\left(\mathrm{\&omega;}\right)-{\stackrel{\&OverBar;}{S}}_{f_i}\left(\mathrm{\&omega;}\right)]}{\underset{l=1}{\overset{k}{\mathrm{\&Sigma;}}}{[S_{f_i^l}\left(\mathrm{\&omega;}\right)-{\stackrel{\&OverBar;}{S}}_{f_i}\left(\mathrm{\&omega;}\right)]}^2}\end{array}\right.---\left(15\right)$

(15)式中，

为载荷点i所施加的k次载荷的平均值，

为响应点j在载荷点i所施加的k次响应的平均值；

步骤3、利用载荷源的不相关特性对载荷识别方程(1)式进行解耦得到（4）式：

$\begin{array}{c}{S}_{\mathrm{yy}}\left(\mathrm{\&omega;}\right)=\frac{1}{2\mathrm{\&pi;}}{\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}h\left(u\right)C_{\mathrm{ff}}(\mathrm{\&tau;}+u-v)\&CenterDot;h^T\left(v\right)e^{-\mathrm{j\&omega;\&tau;}}\mathrm{dudvu\&tau;}\\ =\stackrel{\&OverBar;}{H}\left(\mathrm{\&omega;}\right)S_{\mathrm{ff}}\left(\mathrm{\&omega;}\right)H^T\left(\mathrm{\&omega;}\right)\end{array}---\left(1\right)$

(1)式中h(u)是系统的单位脉冲响应矩阵，h^T(u)是系统的单位脉冲响应矩阵的转置，

是输入各激励之间的协方差函数矩阵，

是系统频率特性矩阵，H^T(ω)是系统频率特性矩阵的转置，

是系统频率特性矩阵的共轭，S_ff(ω)是输入各激励之间的互功率谱密度，

在m个输入载荷激励都是零均值的平稳随机过程，且互不相关的情况下，输入协方差函数矩阵

为对角阵，即：

输入功率谱矩阵S_ff(ω)也为对角阵

此时，输出功率谱中主对角线上的任意一元素 $S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\&omega;}\right)(j=\mathrm{1,2},\&CenterDot;\&CenterDot;\&CenterDot;,n)$ 满足：

$\begin{array}{c}{S}_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\&omega;}\right)=\left[\begin{array}{ccccc}{\stackrel{\&OverBar;}{H}}_{j1}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& {\stackrel{\&OverBar;}{H}}_{\mathrm{ji}}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& {\stackrel{\&OverBar;}{H}}_{\mathrm{jn}}\end{array}\left(\mathrm{\&omega;}\right)\right]\&CenterDot;\mathrm{diag}\left[S_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\&omega;}\right)\right]\&CenterDot;{\left[\begin{array}{ccccc}{H}_{j1}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& H_{\mathrm{ji}}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& H_{\mathrm{jn}}\left(\mathrm{\&omega;}\right)\end{array}\right]}^T\\ =\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\stackrel{\&OverBar;}{H}}_{\mathrm{ji}}\left(\mathrm{\&omega;}\right)S_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\&omega;}\right)H_{\mathrm{ji}}^T\left(\mathrm{\&omega;}\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left|H_{\mathrm{ji}}\left(\mathrm{\&omega;}\right)\right|}^2{S}_{{\mathrm{ff}}_{\mathrm{ii}}}\left(\mathrm{\&omega;}\right)\end{array}---\left(3\right)$

(3)式写成矩阵后的形式为：

其中，|H_j,i(ω)|²是输入f_i对响应y_j的传递函数模的平方，

是待识别的载荷源f_i的自功率谱，是响应y_j的自功率谱；

记 ${\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)\stackrel{\mathrm{\&Delta;}}{=}{\left[\begin{array}{ccc}{S}_{{\mathrm{yy}}_{11}}\left(\mathrm{\&omega;}\right)& \&CenterDot;\&CenterDot;\&CenterDot;& S_{{\mathrm{yy}}_{\mathrm{nn}}}\left(\mathrm{\&omega;}\right)\end{array}\right]}^T$ ，

则(4)式可简写为： ${\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)=B\left(\mathrm{\&omega;}\right){\stackrel{\&RightArrow;}{S}}_F\left(\mathrm{\&omega;}\right);$

步骤4、对(4)式识别出的S_yy(ω)和

采用(16)式传统最小二乘广义逆计算(4)式载荷值：

${\stackrel{\&RightArrow;}{S}}_F^{\&prime;}\left(\mathrm{\&omega;}\right)={[B{\left(\mathrm{\&omega;}\right)}^{T}B\left(\mathrm{\&omega;}\right)]}^{-1}B{\left(\mathrm{\&omega;}\right)}^T{\stackrel{\&RightArrow;}{S}}_Y\left(\mathrm{\&omega;}\right)---\left(16\right)$

步骤5、利用(27)式条件数的定义计算（4）式的条件数，若条件数小于阈值，直接使用（16）式的载荷值作为最终的识别载荷值,载荷识别结束；如果条件数大于阈值，则使用（16）式的计算结果控制粒子群的初始化位置，并随机初始化粒子的初试速度；

在求解方程Ax=b的过程中，引入的舍入误差总会导致其数值解x或多或少地不等于其理论解x'=A^-1B，

定义误差为： $\mathrm{\&delta;x}\stackrel{\mathrm{\&Delta;}}{=}x-x^{\&prime;}---\left(25\right)$

定义剩余为： $\mathrm{\&delta;b}\stackrel{\mathrm{\&Delta;}}{=}\mathrm{Ax}-b---\left(26\right)$

定义条件数为： $\mathrm{cond}\left(A\right)\stackrel{\mathrm{\&Delta;}}{=}\left|\right|A^{-1}\left|\right|\left|\right|A\left|\right|---\left(27\right);$

步骤6、以(20)式作为适应值函数，计算每个粒子的适应值：

$\min f\left(S_{\mathrm{ff}}^{\&prime;}\left(\mathrm{\&omega;}\right)\right)=\min \left(\underset{j}{\max }\right\{\frac{|S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\&omega;}\right)-\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{\left|H_{j,i}\left(\mathrm{\&omega;}\right)\right|}^2{S}_{{\mathrm{ff}}_{\mathrm{ii}}^{\&prime;}}\left(\mathrm{\&omega;}\right)|}{\left|S_{{\mathrm{yy}}_{\mathrm{jj}}}\left(\mathrm{\&omega;}\right)\right|}\left\}\right)---\left(20\right)$

步骤7、找到每个粒子的历史最优位置和最优适应值：

对每个粒子而言，如果它的当前适应度比其历史最优适应度还要小，则使用当前的适应度替代该粒子的最优适应度，并保存当前位置为该粒子历史最优位置；

步骤8、找到整个群体的最优位置和最优适应值；

对每个粒子而言，如果它的当前适应度比整个群体最优适应度还要小，则使用当前的适应度替代粒子群最优适应度，并保存当前位置为整个群体最优位置；

步骤9、根据COMI-PSO算法，使用(21)式、(22)式计算各个粒子最新速度

和最新位置

$v_{\mathrm{id}}^{t+1}=w*v_{\mathrm{id}}^t+c_1{r}_1(p_{\mathrm{id}}-x_{\mathrm{id}}^t)+c_2{r}_2(p_{\mathrm{gd}}-x_{\mathrm{id}}^t)---\left(21\right)$

$x_{\mathrm{id}}^{t+1}=x_{\mathrm{id}}^t+v_{\mathrm{id}}^{t+1}---\left(22\right)$

其中惯性参数w=max_W-(max_W-min_W)*ln(1+(e-1)*CT/TS)，max_W为惯性参数的上限，min_W为惯性参数的下限，e为自然对数的底数，上式能保证w的取值范围是[min_W,max_W]；r₁,r₂为(0,1)之间的随机数，p_id表示当前群体中第i个个体历史最优位置在第d维上数值，p_gd表示当前群体最优位置在第d维上的数值，CT为当前的迭代次数，TS为预设的总迭代次数，加速因子c₁,c₂分别为：

c₁=4*(w-min_W)²/(max_W-min_W)²

c₂=4*(max_W-w)²/(max_W-min_W)²

通过上述操作完成迭代过程中第i个个体从第t代的速度

和位置

更新为第(t+1)代新的速度

和位置

步骤10、通过(20)式计算新位置下，每个粒子的适应度，根据遗传算法中的选择算子，保留种群中适应度较好的部分粒子直接进入下一次迭代，对另一部分适应度较小的粒子采用遗传算法中的交叉算子和突变算子进行预处理，计算预处理之后的子代的适应值，将适应度比父代好的部分相同数目的粒子取代原粒子群中的父辈进入下一次迭代；

步骤11、根据是否达到最大迭代次数或达到最好的适应值判断迭代终止条件，如果达到转入步骤12，否则转入步骤6；

步骤12、输出全局最优解，即群体最优位置，作为最终识别的多个不相关载荷源值,载荷识别算法结束。

Images