[发明专利]复合材料格栅结构低速冲击定位方法

摘要

本发明公开一种复合材料格栅结构低速冲击定位方法，其属于复合材料结构类工程结构健康监测领域。本发明将基于牛顿法的无约束优化算法应用到复合材料格栅结构的低速冲击识别问题中，实现了结构低速冲击的快速识别。利用传感阵列技术，通过对结构施加冲击诱导产生Lamb波信号，并采集复合材料格栅结构的冲击响应信号；以其中某个传感器为基准，采用小波分析手段，提取Lamb波到达其他各个传感器的时间差；利用本发明基于牛顿法的无约束优化算法，求解定位非线性方程组，从而获得冲击源的位置。本发明能快速有效地实现复合材料格栅结构的冲击源识别，识别精度高，具有较好的工程应用价值。

主权项

一种复合材料格栅结构低速冲击定位方法，其特征在于：包括如下步骤步骤一：冲击定位根据时差定位法，在构件表面有规则地布设4个压电传感器PZT1，PZT2，PZT3，PZT4，并对传感信号进行适当的处理，测量出传感器接收到的信号的时间差并换算得到声发射源，即冲击位置，假设冲击源I的坐标为（x_I，y_I），由冲击引发的Lamb波传播的波速为Vg，到达各传感器的时间为t_i，那么通过求解下列方程组就可以得到冲击源的位置：(x₁‑x_I)²+(x‑y_I)²‑(t₁V_g)²＝0(x₂‑x_I)²+(x‑y_I)²‑[(t₁+Δt₁₂)V_g]²＝0   （1）(x₃‑x_I)²+(x‑y_I)²‑[(t₁+Δt₁₃)V_g]²＝0(x₄‑x_I)²+(x‑y_I)²‑[(t₁+Δt₁₄)V_g]²＝0其中，t₂=t₁+Δt₁₂，t₃=t₁+Δt₁₃，t₄=t₁+Δt₁₄，Δt_1j(j=2,3,4)为相对于PZT1的时间差，由公式（1）可知，假设Δt_1j已知，那么就可以通过求解自变量为x＝[x_I y_I t₁ V_g]^T的非线性方程组来获取冲击源的位置及Lamb波的传播速度；步骤二：连续小波变换通过小波变换分析声发射信号，对任意给定的f(t)，其连续小波变换（CWT）定义为：$\mathrm{CWT}(a,b)=\&lt;f\left(t\right)|{\mathrm{\&psi;}}_{a,b}\left(t\right)>=\frac{1}{\sqrt{a}}{\&Integral;}_{-\&infin;}^{+\&infin;}f\left(t\right){\mathrm{\&psi;}}^{*}\left(\frac{t-b}{a}\right)\mathrm{dt}---\left(2\right)$其中a为尺度因子，b为平移因子，ψ^*(t)是母小波ψ(t)的复共轭，ψ_a,b(t)为小波基函数，是由母小波ψ(t)经时间轴的平移、伸缩得到的：${\mathrm{\&psi;}}_{a,b}\left(t\right)=\frac{1}{\sqrt{a}}\mathrm{\&psi;}\left(\frac{t-b}{a}\right)---\left(3\right)$采用复数Morlet小波作为分析小波，其表达式为：${\mathrm{\&psi;}}_{\mathrm{morlet}}\left(t\right)=\frac{1}{\sqrt{\mathrm{\&pi;}{f}_b}}{e}^{i{\mathrm{\&omega;}}_{c}t}{e}^{-\frac{t^2}{f_b}}=\frac{1}{\sqrt{\mathrm{\&pi;}{f}_e}}{e}^{-\frac{t^2}{f_b}}[\cos \left({\mathrm{\&omega;}}_{c}t\right)+i\sin \left({\mathrm{\&omega;}}_{c}t\right)]---\left(4\right)$式中：f_b为小波带宽参数，f_c为小波中心频率参数，对其作傅里叶变换有：${\mathrm{\&psi;}}_{\mathrm{morlet}}\left(\mathrm{\&omega;}\right)={\&Integral;}_{-\&infin;}^{+\&infin;}{\mathrm{\&psi;}}_{\mathrm{morlet}}\left(t\right)e^{-\mathrm{i\&omega;t}}\mathrm{dt}={\&Integral;}_{-\&infin;}^{+\&infin;}\frac{1}{\sqrt{\mathrm{\&pi;}{f}_b}}{e}^{{\mathrm{j\&omega;}}_{c}t}{e}^{-\frac{t^2}{f_b}}{e}^{-\mathrm{i\&omega;t}}\mathrm{dt}---\left(5\right)$令ω_c＝2πf_c,ω_b＝2πf_b，则上式可改写为：${\mathrm{\&psi;}}_{\mathrm{morlet}}\left(\mathrm{\&omega;}\right)=\sqrt{\frac{2}{{\mathrm{\&omega;}}_b}}{\&Integral;}_{-\&infin;}^{+\&infin;}{e}^{-\frac{{2\mathrm{\&pi;t}}^2}{{\mathrm{\&omega;}}_b}+{\mathrm{i\&omega;}}_{c}t}{e}^{-\mathrm{i\&omega;t}}\mathrm{dt}=e^{-\frac{1}{2}\frac{{\mathrm{\&omega;}}_b}{4\mathrm{\&pi;}}{(\mathrm{\&omega;}-{\mathrm{\&omega;}}_c)}^2}---\left(6\right)$同理有，Morlet连续复数小波函数的傅立叶变换：${\mathrm{\&psi;}}_{\mathrm{morlet}-\mathrm{ab}}\left(\mathrm{\&omega;}\right)=\sqrt{a}{e}^{-\mathrm{i\&omega;b}}{e}^{-\frac{1}{2}\frac{{\mathrm{\&omega;}}_b}{4\mathrm{\&pi;}}{(\mathrm{a\&omega;}-{\mathrm{\&omega;}}_c)}^2}---\left(7\right)$由式（4）至（7）得出：Morlet连续复数母小波函数看作是中心在t=0，傅立叶变换的中心在ω＝ω_c的函数；Morlet连续复数小波ψ_morlet‑ab(t)＝ψ_morlet((t‑b)/a)的中心在t=b，其傅立叶变换的中心在ω＝ω_c/a；因此，Morlet连续复数小波变换表示的是f(t)在t=b，ω＝ω_c/a左右的时频成分；CWT系数的实部用于决定尺度的大小，而系数的模平方，亦称尺度谱，表征信号在任意时刻各个尺度下的能量密度，其表达式为：|CWT(a,b)|²＝CWT(a,b)·CWT^*(a,b)   （8）尺度谱系数越大，携带的能量也越高，能量最大处，在尺度谱上表现为能量波峰，所对应的瞬态频率即为被分析信号的主频，即分析频率；而且能量波峰在时域上的映射与应力波的到达时间一致，利用小波变换获取波达时间，从而计算各路传感信号相对于PZT1的时间差Δt_1j，分析频率与尺度因子之间有如下关系：$f=\frac{f_c}{a\&CenterDot;T}---\left(9\right)$其中，f为分析频率，f_c为小波变换的中心频率，T为采样周期；步骤三：基于牛顿法的无约束优化算法考虑非线性方程组（1），通过小波变换获取三个时延Δt₁₂、Δt₁₃、Δt₁₄，还需通过确定Lamb波的群速度进而获得声发射源的位置，其包括：3.1.牛顿法求解非线性方程组假设是二次利普希茨（Lipschitz）连续可微函数，为n维欧几里德空间，则非线性方程组可表示为：F(x)=0   （10）其中F是函数F_i（i=1，2，…）的向量，x是自变量x_j（j=1，2，…）的向量，同理，若x^*满足F(x^*)＝0，且当初始点x⁰充分接近x^*时，按照下述迭代公式产生的序列{xⁿ}收敛于x^*xⁿ⁺¹＝xⁿ+δxⁿ＝xⁿ‑J(xⁿ)^‑1·F(xⁿ)   （11）其中，δx=[‑J(x)^‑1·F(x)]，为牛顿步长；J(x)为雅可比（Jacobian）矩阵，包含了目标函数F(x)对变量的一阶偏导数，公式（11）中F＝[F₁ F₂ F₃ F₄]^T，x＝[x_I y_I t₁ V_g]^T，其雅可比矩阵表达式为：$J\left(x\right)=\frac{\&PartialD;F\left(x\right)}{\&PartialD;x}\left[\begin{array}{cccc}\frac{{\&PartialD;F}_1}{{\&PartialD;x}_I}& \frac{{\&PartialD;F}_1\left(x\right)}{{\&PartialD;y}_I}& \frac{{\&PartialD;F}_1\left(x\right)}{{\&PartialD;t}_1}& \frac{{\&PartialD;F}_1\left(x\right)}{{\&PartialD;V}_g}\\ \frac{{\&PartialD;F}_2\left(x\right)}{{\&PartialD;x}_I}& \frac{{\&PartialD;F}_2\left(x\right)}{{\&PartialD;y}_I}& \frac{{\&PartialD;F}_2\left(x\right)}{{\&PartialD;t}_1}& \frac{{\&PartialD;F}_2\left(x\right)}{{\&PartialD;V}_g}\\ \frac{{\&PartialD;F}_3\left(x\right)}{{\&PartialD;x}_I}& \frac{{\&PartialD;F}_3\left(x\right)}{{\&PartialD;y}_I}& \frac{{\&PartialD;F}_3\left(x\right)}{{\&PartialD;t}_1}& \frac{{\&PartialD;F}_3\left(x\right)}{{\&PartialD;V}_g}\\ \frac{{\&PartialD;F}_4\left(x\right)}{{\&PartialD;x}_I}& \frac{{\&PartialD;F}_4\left(x\right)}{{\&PartialD;y}_I}& \frac{{\&PartialD;F}_4\left(x\right)}{{\&PartialD;t}_1}& \frac{{\&PartialD;F}_4\left(x\right)}{{\&PartialD;V}_g}\end{array}\right]---\left(12\right)$3.2.无约束优化算法采用牛顿法和无约束优化技术相结合的算法来极小化目标函数（亦称性能函数）F：在无约束优化算法中，用F的模平方来表征目标函数：$f\left(x\right)=\frac{1}{2}{\left|\right|F\left(x\right)\left|\right|}^2=\frac{1}{2}F\left(x\right)\&CenterDot;F\left(x\right)---\left(14\right)$其中系数因子1/2是为了计算方便而引入的常数，可以证明，f的根均满足等式f(x^*)＝0，阻尼牛顿法是增加沿牛顿方向的线搜索，其迭代公式为：xⁿ⁺¹＝xⁿ+λ_nδxⁿ  0＜λ≤1   (15)其中λ为由一维搜索得到的最佳步长，其初始值为1，线搜索有精确线搜索和非精确线搜索之分；其中精确线搜索，是指求λ_n使目标函数f沿牛顿方向达到极小，即满足$f(x^n+{\mathrm{\&lambda;}}_n{\mathrm{\&delta;x}}^n)=\underset{\mathrm{\&lambda;}>0}{\min }f(x^n+{\mathrm{\&lambda;}}_n{\mathrm{\&delta;x}}^n)---\left(16\right)$采用Armijo‑Goldstein准则进行非精确线搜索设f(x)可微，取Armijo‑Goldstein准则可表示为：$f(x^n+{\mathrm{\&lambda;}}_n{\mathrm{\&delta;x}}^n)\&le;f\left(x^n\right)+\mathrm{\&alpha;}{\mathrm{\&lambda;}}_n\&dtri;f{\left(x^n\right)}^T\&CenterDot;{\mathrm{\&delta;x}}^n---\left(17\right)$取α=10^‑4，可以证明，此时算法的搜索方向是下降方向，因为：$\begin{array}{c}\&dtri;f{\left(x^n\right)}^T\&CenterDot;{\mathrm{\&delta;x}}^n=\frac{1}{2}\&dtri;[F\left(x^n\right)\&CenterDot;F\left(x^n\right)]\&CenterDot;{\mathrm{\&delta;x}}^n\\ =[F\left(x^n\right)\&CenterDot;J\left(x^n\right)]\&CenterDot;[-J^{-1}\left(x^n\right)\&CenterDot;F\left(x^n\right)]\\ =-F\left(x^n\right)\&CenterDot;F\left(x^n\right)\&lt;0\end{array}---\left(18\right)$采用后退线搜索法，即通过最小化多项式模型来求得λ的值：g(λⁿ)＝f(xⁿ+λ_nδxⁿ)   （19）因此，对任意一个牛顿下降方向δxⁿ，式（19）需满足式（17）和式（18），所以有：$g^{\&prime;}\left({\mathrm{\&lambda;}}^n\right)=\&dtri;f{\left(x^n\right)}^T\&CenterDot;{\mathrm{\&delta;x}}^n---\left(20\right)$一开始，假设模型g给定并且线性，即$g\left(0\right)=f\left(x^n\right),g^{\&prime;}\left(0\right)=\&dtri;f{\left(x^n\right)}^T\&CenterDot;{\mathrm{\&delta;x}}^n---\left(21\right)$令λ₀=0，如果多项式模型满足g(1)＝f(xⁿ+δxⁿ)＞g(0)+αg′(0)   （22）则终止搜索，否则，将通过插入之前求得的g(0)、g(1)和g′(0)三个已知量来构造g(λ)的二次函数模型：g_q(λ)≈[g(1)‑g(0)‑g′(0)]λ²+g′(0)λ+g(0)   （23）通过求上式的最小值，可解出λ₁${\mathrm{\&lambda;}}_1=-\frac{g^{\&prime;}\left(0\right){\mathrm{\&lambda;}}_0^2}{2[g\left(1\right)-g\left(0\right)-g^{\&prime;}\left(0\right)]}---\left(24\right)$如果λ₁太小，则上述二项式的建模就不精确，因此，λ₁<0.1，令λ₁=0.1，若此时g(λ_k)＝f(x^k+λ_kδx^k)仍不满足式（17），则需要进一步后退，考虑三次模型：g_c(λ)＝aλ³+bλ²+g′(0)λ+g(0)   （25）利用之前得到的前两个λ，即λ₀和λ₁，可求得上式的系数$\begin{array}{c}\left[\begin{array}{c}a\\ b\end{array}\right]=\frac{1}{{\mathrm{\&lambda;}}_0^2{\mathrm{\&lambda;}}_1^2({\mathrm{\&lambda;}}_1-{\mathrm{\&lambda;}}_0)}\left[\begin{array}{cc}{\mathrm{\&lambda;}}_0^2& -{\mathrm{\&lambda;}}_1^2\\ -{\mathrm{\&lambda;}}_0^3& {\mathrm{\&lambda;}}_1^3\end{array}\right]\\ \&times;\left[\begin{array}{c}g\left({\mathrm{\&lambda;}}_1\right)-g\left(0\right)-g^{\&prime;}\left(0\right){\mathrm{\&lambda;}}_1\\ g\left({\mathrm{\&lambda;}}_2\right)-g\left(0\right)-g^{\&prime;}\left(0\right){\mathrm{\&lambda;}}_2\end{array}\right]\end{array}---\left(26\right)$则g_c(λ)的极小值点为${\mathrm{\&lambda;}}_2=\frac{-b+\sqrt{b^2-3a{g}^{\&prime;}\left(0\right)}}{3a}---\left(27\right)$在后退过程中，若λ_k>0.5λ_k‑1，取λ_k=0.5λ_k‑1；若λ_k<0.1λ_k‑1，取λ_k=0.1λ_k‑1，这样能保证算法的稳定，多项式后退法的基本过程为：①令λ=1，α=1；②如果满足式（17），则令λ_k=λ停止搜索，输出λ_k；否则，转③；③建立如式（22）模型，并根据（23）求得λ₁，若λ₁<0.1，令λ₁=0.1；若满足式（17），则令λ_k=λ停止搜索，输出λ_k；否则，转④；④建立如式（25）多项式函数，根据式（26）和（27）求得λ_n，若λ_n<0.1，令λ_n=0.1，λ_n>0.5，令λ_n=0.5；若满足式（17），则令λ_k=λ停止搜索，输出λ_k；否则，转④。