CN103913722B - 复合材料格栅结构低速冲击定位方法 - Google Patents

Abstract

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

Description

复合材料格栅结构低速冲击定位方法

技术领域：

本发明涉及一种复合材料格栅结构低速冲击定位方法，属于复合材料结构类工程结构健康监测领域。

背景技术：

格栅结构是由加强肋和蒙皮组成，加强肋呈正多边形网格分布，结构表现各向异性。格栅结构一经问世，就以其拓扑优化性、抗腐蚀能力好、防损伤扩散性能好以及可设计性好等优势受到航空航天领域的高度关注。为了满足航空航天器轻质化的要求，出现了复合材料格栅结构(AGS)。复合材料具有较高的比强度、比模量和抗腐蚀性，因此越来越多的学者将复合材料应用到格栅结构中，并逐渐在航空航天工业中推广应用。但是由于复合材料在外界冲击下，其内部结构易产生不可视的损伤，导致结构刚度下降，影响整体性能，因此很有必要将结构健康监测技术应用到复合材料格栅结构中，以实现结构的自诊断和自修复。

近年来，国内外很多学者对AGS的健康监测方面进行了深入而广泛的研究，而对AGS的损伤定位和损伤评估是结构健康监测领域的一个研究方向。日本东京大学的N.Takeda等人(TakeyaH,OzakiT,TakedaN.StructuralhealthmonitoringofadvancedgridstructureusingmultipointFBGsensors[C]//SmartStructuresandMaterials2005.ProcofSPIE:Vol.5762.2005:204-211.)(AmanoM,TakahashiI,OkabeY,etal.Identificationofdamagelocationinadvancedgridstructuresusingfiberbragggratingsensor[C]//SmartStructuresandMaterials2005.ProcofSPIE:Vol.5765,2005:644-655.)首先对飞机方向舵的格栅结构进行了健康检测的研究，通过对比损伤前后的应变差成功地检测到筋的开裂、筋与蒙皮的脱层；但是该方法只在结构损伤较为严重的分析是有效；南京航空航天大学的陈振英等人(《先进复合材料格栅加筋结构(AGS)的损伤定位研究》，硕士学位论文，南京航空航天大学2011年)应用模态振型曲率方法，成功地检测出复合材料格栅结构的肋板脱粘与断裂损伤，但该方法无法检测出格栅结构的蒙皮损伤；蓝友泽等人(《基于频响函数的复合材料格栅结构损伤检测研究》，硕士学位论文，南京航空航天大学2013年)针对传统模态应变能无法检测AGS蒙皮损伤的缺点，提出了一种基于结构频响函数虚部的损伤检测方法，但该方法极易受环境噪声的影响，因此很难在工程中广泛应用。

发明内容：

为了识别AGS在服役过程中遭受的外在冲击载荷的大小和位置，本发明提出了一种复合材料格栅结构低速冲击定位方法，基于牛顿法的无约束优化算法，通过迭代就能够获取声发射源的位置及Lamb波的群速度，本算法结合了线搜索和多项式后退技术，保证任意初始值开始迭代时都能快速收敛。

本发明采用如下技术方案：一种复合材料格栅结构低速冲击定位方法，其包括如下步骤：

步骤一：冲击定位

根据时差定位法，在构件表面有规则地布设4个压电传感器(分别为PZT1，PZT2，PZT3，PZT4)，并对传感信号进行适当的处理，测量出传感器接收到的信号的时间差并换算得到声发射源，即冲击位置，假设冲击源I的坐标为(x_I，y_I)，由冲击引发的Lamb波传播的波速为Vg，到达各传感器的时间为t_i，那么通过求解下列方程组就可以得到冲击源的位置：

$\begin{array}{c}{\left(x_1-x_I\right)}^2+{\left(y_1-y_I\right)}^2-{\left(t_1{V}_g\right)}^2=0\\ {\left(x_2-x_I\right)}^2+{\left(y_2-y_I\right)}^2-{\[\left(t_1+{\mathrm{\Δt}}_{12}\right)V_g\]}^2=0\\ {\left(x_3-x_I\right)}^2+{\left(y_3-y_I\right)}^2-{\[\left(t_1+{\mathrm{\Δt}}_{13}\right)V_g\]}^2=0\\ {\left(x_4-x_I\right)}^2+{\left(y_4-y_I\right)}^2-{\[\left(t_1+{\mathrm{\Δt}}_{14}\right)V_g\]}^2=0\end{array}---\left(1\right)$

其中，(x₁,y₁)，(x₂,y₂)，(x₃,y₃)，(x₄,y₄)为别为4个压电传感器的坐标，

t₂＝t₁+Δt₁₂，t₃＝t₁+Δt₁₃，t₄＝t₁+Δt₁₄，Δt_1j(j＝2,3,4)为相对于PZT1的时间差，

由公式(1)可知，假设Δt_1j已知，那么就可以通过求解自变量为x＝[x_Iy_It₁V_g]^T的非线性方程组来获取冲击源的位置及Lamb波的传播速度；

步骤二：连续小波变换

通过小波变换分析声发射信号，对任意给定的f(t)，其连续小波变换(CWT)定义为：

$CWT\left(a,b\right)=<f\left(t\right)|{\mathrm{\&psi;}}_{a,b}\left(t\right)>=\frac{1}{\sqrt{a}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}f\left(t\right){\mathrm{\&psi;}}^{*}\left(\frac{t-b}{a}\right)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;}}_{morlet}\left(t\right)=\frac{1}{\sqrt{{\mathrm{\&pi;f}}_b}}{e}^{{\mathrm{i\&omega;}}_{c}t}{e}^{-\frac{t^2}{f_b}}=\frac{1}{\sqrt{{\mathrm{\&pi;f}}_b}}{e}^{-\frac{t^2}{f_b}}\&lsqb;\cos \left({\mathrm{\&omega;}}_{c}t\right)+i\sin \left({\mathrm{\&omega;}}_{c}t\right)\&rsqb;---\left(4\right)$

式中：f_b为小波带宽参数，f_c为小波中心频率参数，对其作傅里叶变换有：

${\mathrm{\&psi;}}_{morlet}\left(\mathrm{\&omega;}\right)={\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}{\mathrm{\&psi;}}_{morlet}\left(t\right)e^{-i\mathrm{\&omega;}t}dt={\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}\frac{1}{\sqrt{{\mathrm{\&pi;f}}_b}}{e}^{{\mathrm{j\&omega;}}_{c}t}{e}^{-\frac{t^2}{f_b}}{e}^{-i\mathrm{\&omega;}t}dt---\left(5\right)$

令ω_c＝2πf_c,ω_b＝2πf_b，则上式可改写为：

${\mathrm{\&psi;}}_{morlet}\left(\mathrm{\&omega;}\right)=\sqrt{\frac{2}{{\mathrm{\&omega;}}_b}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}{e}^{-\frac{2{\mathrm{\&pi;t}}^2}{{\mathrm{\&omega;}}_b}+{\mathrm{i\&omega;}}_{c}t}{e}^{-i\mathrm{\&omega;}t}dt=e^{-\frac{1}{2}\frac{{\mathrm{\&omega;}}_b}{4\mathrm{\&pi;}}{\left(\mathrm{\&omega;}-{\mathrm{\&omega;}}_c\right)}^2}---\left(6\right)$

同理有，Morlet连续复数小波函数的傅立叶变换：

${\mathrm{\&psi;}}_{morlet-ab}\left(\mathrm{\&omega;}\right)=\sqrt{a}{e}^{-i\mathrm{\&omega;}b}{e}^{-\frac{1}{2}\frac{{\mathrm{\&omega;}}_b}{4\mathrm{\&pi;}}{\left(a\mathrm{\&omega;}-{\mathrm{\&omega;}}_c\right)}^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.牛顿法求解非线性方程组

假设F:是二次利普希茨(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)对变量的一阶偏导数；其中F＝[F₁F₂F₃F₄]^T，x＝[x_Iy_It₁V_g]^T，其雅可比矩阵表达式为：

$J\left(x\right)=\frac{\&part;F\left(x\right)}{\&part;x}=\left[\begin{array}{cccc}\frac{\&part;F_1\left(x\right)}{\&part;x_I}& \frac{\&part;F_1\left(x\right)}{\&part;y_I}& \frac{\&part;F_1\left(x\right)}{\&part;t_1}& \frac{\&part;F_1\left(x\right)}{\&part;V_g}\\ \frac{\&part;F_2\left(x\right)}{\&part;x_I}& \frac{\&part;F_2\left(x\right)}{\&part;y_I}& \frac{\&part;F_2\left(x\right)}{\&part;t_1}& \frac{\&part;F_2\left(x\right)}{\&part;V_g}\\ \frac{\&part;F_3\left(x\right)}{\&part;x_I}& \frac{\&part;F_3\left(x\right)}{\&part;y_I}& \frac{\&part;F_3\left(x\right)}{\&part;t_1}& \frac{\&part;F_3\left(x\right)}{\&part;V_g}\\ \frac{\&part;F_4\left(x\right)}{\&part;x_I}& \frac{\&part;F_4\left(x\right)}{\&part;y_I}& \frac{\&part;F_4\left(x\right)}{\&part;t_1}& \frac{\&part;F_4\left(x\right)}{\&part;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)|{|}^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\left(x^n+{\mathrm{\&lambda;}}_n{\mathrm{\&delta;x}}^n\right)=\underset{\mathrm{\&lambda;}>0}{\min }f\left(x^n+{\mathrm{\&lambda;}}_n{\mathrm{\&delta;x}}^n\right)---\left(16\right)$

采用Armijo-Goldstein准则进行非精确线搜索

设f(x)可微，取Armijo-Goldstein准则可表示为：

$f\left(x^n+{\mathrm{\&lambda;}}_n{\mathrm{\&delta;x}}^n\right)\&le;f\left(x^n\right)+{\mathrm{\&alpha;\&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;\&lsqb;F\left(x^n\right)\&CenterDot;F\left(x^n\right)\&rsqb;\&CenterDot;{\mathrm{\&delta;x}}^n\\ =\&lsqb;F\left(x^n\right)\&CenterDot;J\left(x^n\right)\&rsqb;\&CenterDot;\&lsqb;-J^{-1}\left(x^n\right)\&CenterDot;F\left(x^n\right)\&rsqb;\\ =-F\left(x^n\right)\&CenterDot;F\left(x^n\right)<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给定并且线性，即

$\begin{array}{cc}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\end{array}---\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\&lsqb;g\left(1\right)-g\left(0\right)-g^{\&prime;}\left(0\right)\&rsqb;}---\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\left({\mathrm{\&lambda;}}_1-{\mathrm{\&lambda;}}_0\right)}\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{{}^2-3{\mathrm{ag}}^{\&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；否则，转④。

本发明具有如下有益效果：

(1).本发明方法是基于牛顿法的一种无约束优化算法，在迭代过程中采用多项式后退技术并结合Armijo-Goldstein准则来搜索最佳步长，从而获得非线性规划问题的最优解，该算法运行时间短，且从任意初始点开始迭代都能全局收敛并且稳定，克服了传统牛顿法局部收敛的不足；

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

附图说明：

图1为本发明冲击定位示意图。

图2为复合材料格栅结构和传感器位置示意图。

图3为PZT1的冲击响应信号及其频谱。

图4为冲击位置在A3时PZT1的波达时间获取示意图。

具体实施方式：

1.冲击定位

从广义上说，冲击过程属于声发射现象。根据时差定位法可知，如果在构件表面有规则地布设压电传感器(本发明采用4个传感器实现冲击源的识别)，并对传感信号进行适当的处理，测量出传感器接收到的信号的时间差就可以换算得到声发射源，即冲击位置，如图1所示。假设冲击源I的坐标为(x_I，y_I)，由冲击引发的Lamb波传播的波速为Vg，到达各传感器的时间为t_i，那么通过求解下列方程组就可以得到冲击源的位置：

$\begin{array}{c}{\left(x_1-x_I\right)}^2+{\left(y_1-y_I\right)}^2-{\left(t_1{V}_g\right)}^2=0\\ {\left(x_2-x_I\right)}^2+{\left(y_2-y_I\right)}^2-{\&lsqb;\left(t_1+{\mathrm{\&Delta;t}}_{12}\right)V_g\&rsqb;}^2=0\\ {\left(x_3-x_I\right)}^2+{\left(y_3-y_I\right)}^2-{\&lsqb;\left(t_1+{\mathrm{\&Delta;t}}_{13}\right)V_g\&rsqb;}^2=0\\ {\left(x_4-x_I\right)}^2+{\left(y_4-y_I\right)}^2-{\&lsqb;\left(t_1+{\mathrm{\&Delta;t}}_{14}\right)V_g\&rsqb;}^2=0\end{array}---\left(1\right)$

其中，(x₁,y₁)，(x₂,y₂)，(x₃,y₃)，(x₄,y₄)为别为4个压电传感器的坐标，

t₂＝t₁+Δt₁₂，t₃＝t₁+Δt₁₃，t₄＝t₁+Δt₁₄，Δt_1j(j＝2,3,4)为相对于PZT1的时间差，

由公式(1)可知，假设Δt_1j已知，那么就可以通过求解自变量为x＝[x_Iy_It₁V_g]^T的非线性方程组来获取冲击源的位置及Lamb波的传播速度。

2.连续小波变换

由冲击诱导产生的声发射信号是一种非平稳的宽带随机信号且模式多变，传统的特征参数分析法无法满足需求。小波变换是一种多分辨率分析方法，对低频分量采用较宽的时窗，具有较高的频率分辨率；对高频信号采用较窄的时窗，具有较高的时间分辨率，并且遵循Heisenberg不等式，是分析声发射信号最有效的工具之一。

对任意给定的f(t)，其连续小波变换(CWT)定义为：

$CWT\left(a,b\right)=<f\left(t\right)|{\mathrm{\&psi;}}_{a,b}\left(t\right)>=\frac{1}{\sqrt{a}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}f\left(t\right){\mathrm{\&psi;}}^{*}\left(\frac{t-b}{a}\right)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小波具有能够分离幅值和相位、测量瞬时频率、时频窗口可调及能够提供最佳的时频分辨率等优点，采用复数Morlet小波作为分析小波，其表达式为：

${\mathrm{\&psi;}}_{morlet}\left(t\right)=\frac{1}{\sqrt{{\mathrm{\&pi;f}}_b}}{e}^{{\mathrm{i\&omega;}}_{c}t}{e}^{-\frac{t^2}{f_b}}=\frac{1}{\sqrt{{\mathrm{\&pi;f}}_b}}{e}^{-\frac{t^2}{f_b}}\&lsqb;\cos \left({\mathrm{\&omega;}}_{c}t\right)+i\sin \left({\mathrm{\&omega;}}_{c}t\right)\&rsqb;---\left(4\right)$

式中：f_b为小波带宽参数，f_c为小波中心频率参数。对其作傅里叶变换有：

${\mathrm{\&psi;}}_{morlet}\left(\mathrm{\&omega;}\right)={\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}{\mathrm{\&psi;}}_{morlet}\left(t\right)e^{-i\mathrm{\&omega;}t}dt={\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}\frac{1}{\sqrt{{\mathrm{\&pi;f}}_b}}{e}^{{\mathrm{j\&omega;}}_{c}t}{e}^{-\frac{t^2}{f_b}}{e}^{-i\mathrm{\&omega;}t}dt---\left(5\right)$

令ω_c＝2πf_c,ω_b＝2πf_b，则上式可改写为：

${\mathrm{\&psi;}}_{morlet}\left(\mathrm{\&omega;}\right)=\sqrt{\frac{2}{{\mathrm{\&omega;}}_b}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}{e}^{-\frac{2{\mathrm{\&pi;t}}^2}{{\mathrm{\&omega;}}_b}+{\mathrm{i\&omega;}}_{c}t}{e}^{-i\mathrm{\&omega;}t}dt=e^{-\frac{1}{2}\frac{{\mathrm{\&omega;}}_b}{4\mathrm{\&pi;}}{\left(\mathrm{\&omega;}-{\mathrm{\&omega;}}_c\right)}^2}---\left(6\right)$

同理有，Morlet连续复数小波函数的傅立叶变换：

${\mathrm{\&psi;}}_{morlet-ab}\left(\mathrm{\&omega;}\right)=\sqrt{a}{e}^{-i\mathrm{\&omega;}b}{e}^{-\frac{1}{2}\frac{{\mathrm{\&omega;}}_b}{4\mathrm{\&pi;}}{\left(a\mathrm{\&omega;}-{\mathrm{\&omega;}}_c\right)}^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为采样周期。

3.基于牛顿法的无约束优化算法

考虑非线性方程组(1)，三个时延Δt₁₂、Δt₁₃、Δt₁₄可以通过前面分析的小波变换获取，但是要想获得声发射源的位置，还需要确定Lamb波的群速度，而群速度的计算误差又会影响定位的精度。本发明采用一种基于牛顿法的无约束优化算法，该算法无需Lamb波群速度及结构的先验知识，且在任意初始值开始迭代都能保证收敛，具有很高的鲁棒性。

3.1.牛顿法求解非线性方程组的基本原理

假设F:是二次利普希茨(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)对变量的一阶偏导数；其中F＝[F₁F₂F₃F₄]^T，x＝[x_Iy_It₁V_g]^T，其雅可比矩阵表达式为：

$J\left(x\right)=\frac{\&part;F\left(x\right)}{\&part;x}=\left[\begin{array}{cccc}\frac{\&part;F_1\left(x\right)}{\&part;x_I}& \frac{\&part;F_1\left(x\right)}{\&part;y_I}& \frac{\&part;F_1\left(x\right)}{\&part;t_1}& \frac{\&part;F_1\left(x\right)}{\&part;V_g}\\ \frac{\&part;F_2\left(x\right)}{\&part;x_I}& \frac{\&part;F_2\left(x\right)}{\&part;y_I}& \frac{\&part;F_2\left(x\right)}{\&part;t_1}& \frac{\&part;F_2\left(x\right)}{\&part;V_g}\\ \frac{\&part;F_3\left(x\right)}{\&part;x_I}& \frac{\&part;F_3\left(x\right)}{\&part;y_I}& \frac{\&part;F_3\left(x\right)}{\&part;t_1}& \frac{\&part;F_3\left(x\right)}{\&part;V_g}\\ \frac{\&part;F_4\left(x\right)}{\&part;x_I}& \frac{\&part;F_4\left(x\right)}{\&part;y_I}& \frac{\&part;F_4\left(x\right)}{\&part;t_1}& \frac{\&part;F_4\left(x\right)}{\&part;V_g}\end{array}\right]---\left(12\right)$

3.2.无约束优化算法

牛顿法是局部收敛的，即当初始点选择不当时，可能会导致算法的振荡和发散，为了使算法在任意初始点xⁿ全局收敛，本文采用牛顿法和无约束优化技术相结合的算法来极小化目标函数(亦称性能函数)F：

在无约束优化算法中，经常用F的模平方来表征目标函数：

$f\left(x\right)=\frac{1}{2}\left|\right|F\left(x\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\left(x^n+{\mathrm{\&lambda;}}_n{\mathrm{\&delta;x}}^n\right)=\underset{\mathrm{\&lambda;}>0}{\min }f\left(x^n+{\mathrm{\&lambda;}}_n{\mathrm{\&delta;x}}^n\right)---\left(16\right)$

其基本思想是：首先确定包含问题最优解的搜索空间，然后采用某种插值或分割技术缩小这个空间，进行搜索求解；而非精确线搜索，是指选取λ_n使目标函数f得到可接受的下降量，即Δf_n＝f(xⁿ)-f(xⁿ+λ_nδxⁿ)＞0是可接受的。但计算量大、耗费资源多时精确搜索的主要不足，特别是当迭代点远离最优点时，精确线搜索通常不能达到预期效果。此时，非精确线搜索应运而生，其基本思想是：在当前迭代点xⁿ确定了下降方向δxⁿ之后，只需后继迭代点xⁿ⁺¹＝xⁿ+λ_nδxⁿ使得目标函数有满意的下降量即可。Armijo-Goldstein和Wolf-Powell准则是衡量该下降量常用的准则。本发明采用Armijo-Goldstein准则进行非精确线搜索。

设f(x)可微，取Armijo-Goldstein准则可表示为：

$f\left(x^n+{\mathrm{\&lambda;}}_n{\mathrm{\&delta;x}}^n\right)\&le;f\left(x^n\right)+{\mathrm{\&alpha;\&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;\&lsqb;F\left(x^n\right)\&CenterDot;F\left(x^n\right)\&rsqb;\&CenterDot;{\mathrm{\&delta;x}}^n\\ =\&lsqb;F\left(x^n\right)\&CenterDot;J\left(x^n\right)\&rsqb;\&CenterDot;\&lsqb;-J^{-1}\left(x^n\right)\&CenterDot;F\left(x^n\right)\&rsqb;\\ =-F\left(x^n\right)\&CenterDot;F\left(x^n\right)<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给定并且线性，即

$\begin{array}{cc}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\end{array}---\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\&lsqb;g\left(1\right)-g\left(0\right)-g^{\&prime;}\left(0\right)\&rsqb;}---\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\left({\mathrm{\&lambda;}}_1-{\mathrm{\&lambda;}}_0\right)}\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-3{\mathrm{ag}}^{\&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；否则，转④。

4.功能验证实验

4.1实验设置

实验对象为玻璃/环氧复合材料格栅结构，如图2所示。材料选用纤维含量45％的正交各向异性板结构铺层形式为玻璃/环氧复合材料，其中蒙皮使用[-45°,90°,45°,0°,-45°,90°,45°,0°]_s16层角对称铺设层合板，肋板采用[-45°,90°,45°,0°]_s8层角对称铺设单层板，各单层厚度为0.2mm，具体力学参数见表1。

表1玻璃/环氧复合材料力学参数

复合材料格栅结构板尺寸为长×宽×高＝330mm×156mm×30mm，肋板与蒙皮厚度均为t＝2mm，蒙皮闭室为边长a＝60mm的等边三角形。四边通过C型夹具固定在支架上。传感器为PZT-51型压电陶瓷片，尺寸为Ф8mm×0.5mm(本发明均采用此类型传感器，后面不再赘述)。以板的左下角为原点建立直角坐标系，压电片坐标如图2所示。实验中使用PCB力锤模拟冲击，并将冲击响应信号送入电荷放大器，放大后输入PXI控制器，在PXI中设置采样模式为外部触发模式，当信号幅值大于一定阈值时，数据采集系统被触发获取冲击信号，最后通过在PXI控制器中编写相关信号处理及定位算法的程序，实现冲击源的定位。

4.2典型冲击信号

在进行信号采集时，将1～4号压电片分别接至数采卡的1～4通道，力锤接至0通道，并将该通道作为触发通道，采样频率为1MHz。图3为对结构原点进行冲击时，PZT1的冲击响应信号及其频谱。从图中可以看出，信号的能量主要集中在0～20KHz这个范围内，高频信号的能量很小。低频部分的信号其波长越长，时间分辨率也越低，且信号能量较大，容易产生波包混叠现象，不利于波达时间的判断；高频部分的信号时间分辨率教高，但能量太小且模式复杂，亦不利于窄带信号的提取。综合考虑信号能量、时频分辨率和信号模式这三个方面的因素，在小波变换时，取其中心频率为20KHz，后面不再赘述。

本实验针对复合材料格栅结构上的四个不同位置A1(150,90)mm、A2(150,60)mm、A3(90,60)mm、A4(120,0)mm，利用PCB力锤进行冲击诱导产生Lamb波，用小波变换提取中心频率为20KHz的窄带信号，各传感器的波达时间见表2。图4为冲击位置在A3时PZT1的波达时间获取示意图。各传感器的波达时间差见表3。

表2监测区域内几个典型位置冲击后各传感器波达时间

表3冲击位置在A1-A4时各传感器的波达时间差

将采集得到的数据保存并导入到MATLAB中，运行基于牛顿法的优化算法程序，得到的冲击源位置及波速结果见表4。初始迭代点为x＝[0001200]，精度误差为10^-6。从表中可以看出，该算法的最大定位误差不超过1cm，这一误差水平与目前已有的精度较高的冲击定位方法的误差水平相当。而且，该算法的目标定位时间在1s以内，这也意味着该算法可以运用于冲击载荷定位的实时操作系统中。

表4复合材料格栅结构无约束优化算法定位结果

由实验结果可知，基于牛顿法的无约束优化算法是一种全局收敛的算法，对任意一个初始迭代点，都能保证全局收敛。该算法不需要冲击Lamb波的波速以及结构属性等的先验知识，从而也就避免了由波速误差带来的定位误差。实验结果表明，该算法对复合材料格栅结构和某飞行器机翼后缘模型上的冲击载荷的定位误差均在1cm以内，目标定位时间不超过1s，适用于冲击载荷的实时监测系统；并且，该算法同样能识别对发生在监测区域外的冲击。

本发明基于牛顿法和无约束优化技术实现了复合材料格栅结构的冲击载荷定位。通过布设在复合材料格栅结构上的4个压电传感器来采集冲击响应信号，并以PZT1为基准，采用连续小波变换(CWT)计算各个传感器信号相对于PZT1的时延，最后将计算得到的结果输入牛顿法中，结合线搜索和多项式后退技术实现冲击源的定位。实验结果表明，该算法能够快速、准确地识别出冲击位置。因此本发明方法在航空航天飞行器的冲击监测和结构健康监测中有较好地应用前景。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下还可以作出若干改进，这些改进也应视为本发明的保护范围。

Claims (1)

1.一种复合材料格栅结构低速冲击定位方法，其特征在于：包括如下步骤

步骤一：冲击定位

根据时差定位法，在构件表面有规则地布设4个压电传感器PZT1，PZT2，PZT3，PZT4，并对传感信号进行适当的处理，测量出传感器接收到的信号的时间差并换算得到声发射源，即冲击位置，假设冲击源I的坐标为(x_I，y_I)，由冲击引发的Lamb波传播的波速为Vg，到达各传感器的时间为t_i，那么通过求解下列方程组就可以得到冲击源的位置：

$\begin{array}{c}{(x_1-x_I)}^2+{(y_1-y_I)}^2-{\left(t_1{V}_g\right)}^2=0\\ {(x_2-x_I)}^2+{(y_2-y_I)}^2-{\&lsqb;(t_1+{\mathrm{\&Delta;t}}_{12})V_g\&rsqb;}^2=0\\ {(x_3-x_I)}^2+{(y_3-y_I)}^2-{\&lsqb;(t_1+{\mathrm{\&Delta;t}}_{13})V_g\&rsqb;}^2=0\\ {(x_4-x_I)}^2+{(y_4-y_I)}^2-{\&lsqb;(t_1+{\mathrm{\&Delta;t}}_{14})V_g\&rsqb;}^2=0\end{array}---\left(1\right)$

其中，(x₁,y₁)，(x₂,y₂)，(x₃,y₃)，(x₄,y₄)为别为4个压电传感器的坐标，

t₂＝t₁+Δt₁₂，t₃＝t₁+Δt₁₃，t₄＝t₁+Δt₁₄，Δt_1j(j＝2,3,4)为相对于PZT1的时间差，

由公式(1)可知，假设Δt_1j已知，那么就可以通过求解自变量为x＝[x_Iy_It₁V_g]^T的非线性方程组来获取冲击源的位置及Lamb波的传播速度；

步骤二：连续小波变换

通过小波变换分析声发射信号，对任意给定的f(t)，其连续小波变换(CWT)定义为：

$CWT(a,b)=<f\left(t\right)|{\mathrm{\&psi;}}_{a,b}\left(t\right)>=\frac{1}{\sqrt{a}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}f\left(t\right){\mathrm{\&psi;}}^{*}\left(\frac{t-b}{a}\right)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;}}_{morlet}\left(t\right)=\frac{1}{\sqrt{{\mathrm{\&pi;f}}_b}}{e}^{{\mathrm{i\&omega;}}_{c}t}{e}^{-\frac{t^2}{f_b}}=\frac{1}{\sqrt{{\mathrm{\&pi;f}}_b}}{e}^{-\frac{t^2}{f_b}}\&lsqb;cos\left({\mathrm{\&omega;}}_{c}t\right)+isin\left({\mathrm{\&omega;}}_{c}t\right)\&rsqb;---\left(4\right)$

式中：f_b为小波带宽参数，f_c为小波中心频率参数，对其作傅里叶变换有：

${\mathrm{\&psi;}}_{morlet}\left(\mathrm{\&omega;}\right)={\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}{\mathrm{\&psi;}}_{morlet}\left(t\right)e^{-i\mathrm{\&omega;}t}dt={\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}\frac{1}{\sqrt{{\mathrm{\&pi;f}}_b}}{e}^{{\mathrm{j\&omega;}}_{c}t}{e}^{-\frac{t^2}{f_b}}{e}^{-i\mathrm{\&omega;}t}dt---\left(5\right)$

令ω_c＝2πf_c,ω_b＝2πf_b，则上式可改写为：

${\mathrm{\&psi;}}_{morlet}\left(\mathrm{\&omega;}\right)=\sqrt{\frac{2}{{\mathrm{\&omega;}}_b}}{\&Integral;}_{-\mathrm{\&infin;}}^{+\mathrm{\&infin;}}{e}^{-\frac{2{\mathrm{\&pi;t}}^2}{{\mathrm{\&omega;}}_b}+{\mathrm{i\&omega;}}_{c}t}{e}^{-i\mathrm{\&omega;}t}dt=e^{-\frac{1}{2}\frac{{\mathrm{\&omega;}}_b}{4\mathrm{\&pi;}}{(\mathrm{\&omega;}-{\mathrm{\&omega;}}_c)}^2}---\left(6\right)$

同理有，Morlet连续复数小波函数的傅立叶变换：

${\mathrm{\&psi;}}_{morlet-ab}\left(\mathrm{\&omega;}\right)=\sqrt{a}{e}^{-i\mathrm{\&omega;}b}{e}^{-\frac{1}{2}\frac{{\mathrm{\&omega;}}_b}{4\mathrm{\&pi;}}{(a\mathrm{\&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.牛顿法求解非线性方程组

假设F:是二次利普希茨(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_Iy_It₁V_g]^T，其雅可比矩阵表达式为：

$J\left(x\right)=\frac{\&part;F\left(x\right)}{\&part;x}=\left[\begin{array}{cccc}\frac{\&part;F_1\left(x\right)}{\&part;x_I}& \frac{\&part;F_1\left(x\right)}{\&part;y_I}& \frac{\&part;F_1\left(x\right)}{\&part;t_1}& \frac{\&part;F_1\left(x\right)}{\&part;V_g}\\ \frac{\&part;F_2\left(x\right)}{\&part;x_I}& \frac{\&part;F_2\left(x\right)}{\&part;y_I}& \frac{\&part;F_2\left(x\right)}{\&part;t_1}& \frac{\&part;F_2\left(x\right)}{\&part;V_g}\\ \frac{\&part;F_3\left(x\right)}{\&part;x_I}& \frac{\&part;F_3\left(x\right)}{\&part;y_I}& \frac{\&part;F_3\left(x\right)}{\&part;t_1}& \frac{\&part;F_3\left(x\right)}{\&part;V_g}\\ \frac{\&part;F_4\left(x\right)}{\&part;x_I}& \frac{\&part;F_4\left(x\right)}{\&part;y_I}& \frac{\&part;F_4\left(x\right)}{\&part;t_1}& \frac{\&part;F_4\left(x\right)}{\&part;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)|{|}^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;\&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;\&lsqb;F\left(x^n\right)\&CenterDot;F\left(x^n\right)\&rsqb;\&CenterDot;{\mathrm{\&delta;x}}^n\\ =\&lsqb;F\left(x^n\right)\&CenterDot;J\left(x^n\right)\&rsqb;\&CenterDot;\&lsqb;-J^{-1}\left(x^n\right)\&CenterDot;F\left(x^n\right)\&rsqb;\\ =-F\left(x^n\right)\&CenterDot;F\left(x^n\right)<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(0)＝f(xⁿ) $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\&lsqb;g\left(1\right)-g\left(0\right)-g^{\&prime;}\left(0\right)\&rsqb;}---\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-3{\mathrm{ag}}^{\&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；否则，转④。