CN103267907A - 一种变压器绕组模态参数识别方法 - Google Patents

Abstract

本发明公开了一种变压器绕组模态参数的识别方法，包括下列步骤：采集变压器绕组各个测点的振动信号；对振动信号进行傅里叶变换，得到变压器绕组的振动频响曲线；做傅里叶反变换，得到变压器绕组的自由振动信号；使用Morlet小波对自由振动信号进行小波变换，做小波变换时频图；提取小波变换时频图中的小波脊线，各条小波脊线的纵坐标为变压器绕组的各阶固有频率；计算各阶固有频率的阻尼比；计算各阶固有频率的绕组振型。本发明能准确、高效地识别变压器绕组的模态参数，方便对绕组结构进行设计及制造、绕组变形检测和状态评估。

Description

一种变压器绕组模态参数识别方法

技术领域

本发明涉及一种变压器参数识别技术领域，尤其是涉及一种变压器绕组模态参数识别方法。 

背景技术

变压器是电力系统中最重要的设备之一，其运行的稳定性对电力系统安全影响重大。目前变电站设备及线路的运行环境始终不容乐观，因外部短路造成变压器绕组受冲击而引发的变形是变压器运行过程中较为常见的故障，对系统的安全运行造成了很大的威胁。 

变压器遭受突发短路后，巨大的短路电动力的作用使得绕组振动加剧，使得绕组可能首先发生松动或变形。大量的实验研究分析表明，变压器绕组变形具有累积效应，如果对于松动或变形不能及时发现和修复，那么在变压器的松动或变形累积到一定程度后会使变压器的抗短路能力大幅下降而在遭受较小的冲击电流下也会引发大的事故发生。绕组的变形一方面会导致机械抗短路电流冲击能力的下降，另一方面也会导致线圈内部局部绝缘距离发生变化，使局部出现绝缘薄弱点，当遇到过电压作用时，绕组有可能发生饼间或匝间短路导致变压器绝缘击穿事故，或者由于局部场强增大而引起局部放电，绝缘损伤部位会逐渐扩大，最终导致变压器发生绝缘击穿事故而引发进一步的事态扩大。此外，当短路电动力的激励频率与绕组的固有频率相近时，可能会产生共振现象，将导致绕组的严重变形、绝缘损坏乃至轴向失稳造成绕组坍塌。 

绕组结构件振动的主要描述参数是固有频率（特征值）和对应的振型（特征矢量），绕组的松动、变形或振动加剧等现象均意味着绕组机械动力性特性的变化，因此，准确识别变压器绕组的模态参数（固有频率、阻尼和振型），使其远离激励频率对变压器的设计与制造意义重大，这是保证变压器可靠运行的重要前提。对运行中的变压器来说，当变压器经历了外部短路事故后或 运行一段时间后的常规检修中，如何有效地检测出变压器绕组是否存在松动和变形，从而判断变压器是否需要检修处理显得十分重要，是保障变压器安全运行的一个重要手段。 

现有的模态参数识别方法如多参考最小二乘法（Multi-reference Lease Squares，LSCE）、多参考最小二乘复频域法（Poly-reference Least Squares Complex Frequency Domain Method，PolyMax）和多参考极大似然估计法（Poly-reference Maximum Likelihood Estimation，PMLE）（郑锦涛，车身试验模态分析方法对比研究，华南理工大学硕士学位论文，2012年5月）在识别变压器绕组这类结构复杂的非线性系统时精度有限，且对噪声较为敏感，需要寻求新的方法。 

发明内容

本发明的目的是克服上述现有技术的不足之处，提供一种变压器绕组模态参数识别方法，用于变压器绕组的固有频率、阻尼以及振型的识别。 

为了实现上述发明目的，本发明的技术解决方案如下： 

一种变压器绕组模态参数识别方法，其特点在于，该方法包括下列步骤： 

(1)在变压器绕组上放置激振器，向激振器输入信号V_i对变压器绕组进行激振，采用放置在变压器绕组表面的N个振动加速度传感器采集和记录各个测点的振动信号Voi(i=1,2,…,N)。 

(2)分别对输入的白噪声信号V_i和各个测点的振动信号V_oi进行傅里叶变换（傅里叶变换是本领域内常用的数学方法，因此发明人在此不再进行详细的描述），得到变压器绕组的振动频响曲线H(ω)： 

$H\left(\mathrm{\ω}\right)=\left(\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{V}_{\mathrm{oi}}\left(\mathrm{\ω}\right)\right)/V_i\left(\mathrm{\ω}\right)$

式中，V_i(ω)为输入的白噪声信号的傅里叶变换；V_oi(ω)为各个测点振动信号的傅里叶变换；N为振动加速度传感器的个数。 

(3)对振动频响曲线H(ω)做傅里叶反变换，得到变压器绕组的自由振动信号h(t)（傅里叶反变换是本领域内常用的数学方法，因此发明人在此不再进行详细的描述），此处自由振动信号h(t)可表示为 

$h\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{A}_i{e}^{{-\mathrm{\ϵ}}_i{\mathrm{\ω}}_{\mathrm{ni}}t}\cos ({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)$

式中，P为系统模态阶数；A_i为第i阶系统模态响应信号振幅；ω_ni为第i阶系统无阻尼固有频率；ε_i为第i阶系统阻尼比；θ_i为第i阶模态响应初始相位；ω_di为第i阶系统有阻尼固有频率，且有因此，当系统阻尼较小时，有ω_di＝ω_ni。 

(4)使用Morlet小波对自由振动信号h(t)进行小波变换，得到自由振动信号h(t)的小波变换系数矩阵WT(a,b)，其中 

Z(t)＝h(t)+jH(h(t)) 

$=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{A}_i{e}^{-{\mathrm{\ϵ}}_i{\mathrm{\ω}}_{\mathrm{ni}}t}\cos ({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)+j\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{A}_i{e}^{-{\mathrm{\ϵ}}_i{\mathrm{\ω}}_{\mathrm{ni}}t}\sin ({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)$

$=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{A}_i{e}^{-{\mathrm{\ϵ}}_i{\mathrm{\ω}}_{\mathrm{ni}}t}[\cos ({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)+j\sin ({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)]=\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{A}_i{e}^{-{\mathrm{\ϵ}}_i{\mathrm{\ω}}_{\mathrm{ni}}t}{e}^{j({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)}$

式中，H(h(t))为h(t)的Hilbert变换（Hilbert变换对于本领域内技术人员来说是熟知的，故本文不再对Hilbert变换进行详细描述）。 

$\mathrm{WT}(a,b)=\⟨h\left(t\right),{\mathrm{\ψ}}_{a,b}\left(t\right)\⟩=\frac{1}{2}\⟨Z\left(t\right),{\mathrm{\ψ}}_{a,b}^{*}\left(t\right)\⟩$

$=\frac{1}{2\sqrt{2\mathrm{a\π}}}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\underset{-\∞}{\overset{+\∞}{\∫}}{A}_i{e}^{-{\mathrm{\ϵ}}_i{\mathrm{\ω}}_{\mathrm{ni}}t}{e}^{j({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)}{e}^{(\frac{t-b}{a}{)}^2}{e}^{-j{\mathrm{\ω}}_0\left(\frac{t-b}{a}\right)}\mathrm{dt}$

$=\frac{1}{2\sqrt{2\mathrm{a\π}}}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\underset{-\∞}{\overset{+\∞}{\∫}}{B}_i{e}^{-{\mathrm{\ϵ}}_i{\mathrm{\ω}}_{\mathrm{ni}}t}{e}^{j({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)}{e}^{(\frac{t-b}{a}{)}^2}{e}^{-j{\mathrm{\ω}}_0\left(\frac{t-b}{a}\right)}\mathrm{dt}$

此处，令 $B_i\left(t\right)=A_i{e}^{-{\mathrm{\ϵ}}_i{\mathrm{\ω}}_{\mathrm{ni}}t}.$

在t＝b处对上式进行Taylor级数展开（Taylor级数对于本领域内技术人员来说是熟知的，故本文不再对Taylor级数进行详细描述），有 

$\mathrm{WT}(a,b)=\frac{1}{2\sqrt{2\mathrm{a\π}}}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\underset{-\∞}{\overset{+\∞}{\∫}}\{B_i\left(b\right)+O[B_i^{\′}\left(b\right)\left]\right\}{e}^{j({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)}{e}^{(\frac{t-b}{a}{)}^2}{e}^{-j{\mathrm{\ω}}_0\left(\frac{t-b}{a}\right)}\mathrm{dt}$

忽略上式中的高阶无穷小项，有 

$\mathrm{WT}(a,b)\≈\frac{1}{2\sqrt{2\mathrm{a\π}}}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}\underset{-\∞}{\overset{+\∞}{\∫}}{B}_i\left(b\right)e^{j({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)}{e}^{(\frac{t-b}{a}{)}^2}{e}^{-j{\mathrm{\ω}}_0\left(\frac{t-b}{a}\right)}\mathrm{dt}$

$\≈\frac{\sqrt{\mathrm{a\π}}}{2}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{B}_i\left(b\right)e^{j({\mathrm{\ω}}_{\mathrm{di}}t+{\mathrm{\θ}}_i)}{e}^{-({\mathrm{a\ω}}_{\mathrm{di}}{\mathrm{\ω}}_0{)}^{2/2}}\≈\frac{\sqrt{\mathrm{a\π}}}{2}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{A}_i{e}^{-{\mathrm{\ϵ}}_i{\mathrm{\ω}}_{\mathrm{ni}}b}{e}^{-(a{\mathrm{\ω}}_{\mathrm{di}}-{\mathrm{\ω}}_0{)}^{2/2}}{e}^{j({\mathrm{\ω}}_{\mathrm{di}}b+{\mathrm{\θ}}_i)}$

式中，WT(a,b)表示小波变换系数矩阵中第a行和第b列对应的元素，其中小波变换系数矩阵的总行数A对应频率，小波变换系数矩阵的总列数B对应 时间；e为自然对数的底数，是一常数，取值为2.71828；j为虚数单位，为

ω₀为Morlet小波的中心频率，即为Morlet小波窗口频率的中心，为一常数。 

上述Morlet小波是本领域内技术人员所熟知的，其作为小波基是一系列正交的小波函数，是在时域和频域都具有紧支集且均值为零的函数。若设ψ(t)为平方可积函数，且ψ(t)的傅里叶变换ψ(ω)满足下述容许条件： 

${\mathrm{\ψ}}_{(a,b)}=\frac{1}{a}\mathrm{\ψ}\left(\frac{t-b}{a}\right)b\∈R,a>0$

式中，a为频率尺度；b为时间。 

那么Morlet小波基的时域和频域表达式分别为 

式中，ω₀为Morlet小波的中心频率，为一常数。 

由于Morlet小波对于本领域内技术人员来说是熟知的，故本文不再对Morlet小波进行详细描述，上述内容只是对Morlet小波进行的简单介绍，并不作为对本技术方案的限制。 

(5)对小波变换系数矩阵WT(a,b)取模，计算小波变换系数矩阵WT(a,b)的模值|WT(a,b)|，根据小波变换系数矩阵的模值|WT(a,b)|做小波变换时频图：即以时间为横坐标轴，频率为纵坐标轴，以小波变换系数矩阵WT(a,b)的模值为显示结果作图，就可以得到小波变换时频图，其中 

$\left|\mathrm{WT}(a,b)\right|\≈\frac{\sqrt{\mathrm{a\π}}}{2}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{A}_i{e}^{-{\mathrm{\ϵ}}_i{\mathrm{\ω}}_{\mathrm{ni}}b}{e}^{-({\mathrm{a\ω}}_{\mathrm{di}}{\mathrm{\ω}}_0{)}^{2/2}}$

(6)提取小波变换时频图中的小波脊线（小波脊线为时频图中呈现出的类似地形图中山脊的形状的点集连线，其是本领域内技术人员知晓的技术术语），小波变换时频图中各条小波脊线的纵坐标为变压器绕组的各阶固有频率。 

提取小波变换时频图中的小波脊线可以采用如下疯狂爬山法： 

6a.定义并初始化度量密度矩阵D，初始值为0，度量密度矩阵D的行列数与小波变换系数矩阵WT(a,b)的行列数相同，在度量密度矩阵D所有元素中，随机选取n1个元素作为运动点，此处，n1为小波变换系数矩阵WT(a,b)元素总个数的四分之一；定义系统循环变量为T_t，其初值为 T₀＝max(WT)-min(WT)，其中，max(WT)和min(WT)分别为小波变换系数矩阵WT(a,b)的模极大值和模极小值（该模极大值和模极小值是从上述步骤（4）中计算得到的一系列小波变换系数矩阵WT(a,b)的模值中获得的，也就是说上述小波变换系数矩阵WT(a,b)的模值为一系列值，其中包括有一极大值和一极小值）； 

6b.在初始时刻t＝t₀，确定运动点X_k对应的位置为(i,j)，即第k个运动点在t₀时刻满足X_k(t₀)＝(i,j)，此处有0＜k≤n1，i为第k个运动点在小波变换时频图中的横坐标位置，j为第k个运动点在小波变换时频图中的纵坐标位置，且1≤i≤A和1≤j≤B；A为小波变换系数矩阵WT(a,b)的行数，B为小波变换系数矩阵WT(a,b)的列数； 

6c.在下一时刻即t+1时刻， 

在水平方向上， 

若满足1＜i＜A，则以概率P＝50%向左或向右移动运动点X_k(t₀)到位置(i,j′)，即有 

j′＝j+1或j′＝j-1 

若满足i＝1，则向右移动运动点X_k(t₀)到位置(i,j′)，j′＝j+1； 

若满足i＝B，则向左移动运动点X_k(t₀)到位置(i,j′)，j′＝j-1； 

在垂直方向上， 

若满足j′＝1，则向下移动初始运动点X_k(t₀)到位置(i′,j′)，i′＝i+1； 

若满足j′＝A，则向上移动初始运动点X_k(t₀)到位置(i′,j′)，i′＝i-1； 

若满足1＜j′＜B，先以概率P＝50%的概率向上或向下移动运动点X_k(t₀)到位置(i_p,j′)，即有 

i_p＝i+1或i_p＝i-1 

若满足WT(i_p,j′)＞WT(i,j′)，则运动点X_k(t₀)在当前时刻的位置为 

(i′,j′)，i′＝i_p，转至步骤6d； 

若满足WT(i_p,j′)＜WT(i,j′)，根据下式计算垂直方向上的移动概率p： 

$p=\exp \left[\frac{\mathrm{wt}(i_p,j^{\′})-\mathrm{wT}(i,j^{\′})}{T_t}\right]$

如上式所定义的，T_t为系统循环变量； 

将运动点X_k(t₀)从位置(i_p,j′)按照概率p保持不动，即i′＝i_p；按照概率 1-p从位置(i_p,j′)返回到位置(i,j′)，即i′＝i。 

6d.记录该运动点X_k在t+1时刻的位置，即为(i′,j′)，根据下式更新度量密度矩阵，即有 

D_t+1(i′,j′)＝D_t(i′,j′)+WT(i′,j′) 

6e.对所有n1个点重复步骤6a～6d； 

6f.计算系统当前循环参量为T_t＝T_t/(t+1)²，对所有n1个点重复步骤6a～6e，直到T_t＜T₀/1000； 

6g.确定最后的密度度量矩阵D′(i,j)，其中D′(i,j)中的各个元素为脊点： 

$D^{\&prime;}(i,j)=\left\{\begin{array}{cc}D(i,j)& D(i,j)\&GreaterEqual;T_h\\ 0& D(i,j)<T_h\end{array}\right.$

式中，T_h为设定的密度阈值； 

6h.根据脊点做出小波脊线，即以时间为横坐标轴，频率为纵坐标轴，最终密度度量矩阵D′(i,j)各个元素的模值为显示结果作图，得到小波脊线。 

(7)计算小波系数矩阵WT(a,b)中第a行各个元素的自然对数值，此处，a的取值为步骤(6)中小波脊线的纵坐标，计算小波系数矩阵WT(a,b)中第a行各个元素所在直线对时间的斜率λ，则可计算得到步骤(6)中各阶固有频率对应的阻尼比ε，其计算公式为 

ε_k＝-λ/(2πf_k) 

式中，f_k为第k阶固有频率；ε_k为第k阶固有频率对应的阻尼比。 

(8)根据小波变换系数矩阵WT(a,b)识别步骤(6)中各个固有频率对应的振型，计算过程为： 

(8a).选取小波变换系数矩阵WT(a,b)的模值最小的点作为参考点，记为r点； 

(8b).使用最小二乘法计算r点对应的第m阶小波变换系数WT_r(a,b)的最佳常数解，最小二乘法是本领域内技术人员来说是熟知的，故本文不再对最小二乘法进行详细描述。 

(8c).使用最小二乘法计算第p点对应的第m阶小波变换系数WT_p(a,b)的最佳常数解，此处，p点指除了r点之外的任一点，则可得第p个响应点对参考点r的最佳最小二乘解即振动信息，其计算公式为 

${\mathrm{\&phi;}}_k^j=\frac{{\mathrm{WT}}_r^T{\mathrm{WT}}_p}{{\mathrm{WT}}_r^T{\mathrm{WT}}_r}$

式中，为r点的小波系数变换矩阵的转置；WT_p为p点的小波系数变换矩阵；WT_r为r点的小波系数变换矩阵。 

依次计算各个响应点对参考点的最佳最小二乘解，即可得第k阶固有频率对应的振型。 

进一步地，在本发明所述的变压器绕组模态参数的识别方法中，向激振器输入的信号V_i为白噪声信号。 

与现有方法相比，本发明变压器绕组模态参数方法通过测试变压器绕组的振动频响曲线使用复小波变换法及疯狂爬山算法对变压器绕组的模态参数进行识别，从而得到较为准确的变压器绕组模态参数，以便于对变压器绕组进行设计制造，同时可作为变压器绕组状态监测的基准值，进而能够及时发现问题，对变压器进行及时检修。 

具体实施方式

下面将结合附图和实施例对本发明变压器绕组运行状态的诊断方法做进一步的详细说明。 

本实施例中，按照下述步骤对一台10kV的变压器绕组进行监测诊断： 

(1)在该变压器绕组上放置激振器，将经过功率放大器放大的20kHz的白噪声信号V_i输入激振器对变压器绕组进行激振，在该变压器绕组表面放置20个振动加速度传感器采集和记录各个测点的振动信号V_oi(i＝1,2,…,20)，采集时间为0.04s。 

(2)分别对输入的白噪声信号V_i和各个测点的振动信号V_oi进行傅里叶变换，得到变压器绕组的振动频响曲线H(ω)： 

$H\left(\mathrm{\&omega;}\right)=\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{V}_{\mathrm{oi}}\left(\mathrm{\&omega;}\right)\right)/V_i\left(\mathrm{\&omega;}\right)$

式中，V_i(ω)为输入的白噪声信号的傅里叶变换；V_oi(ω)为各个测点振动信号的傅里叶变换；N=20。 

(3)对振动频响曲线H(ω)做傅里叶反变换，得到变压器绕组的自由振动信号h(t)（傅里叶反变换是本领域内常用的数学方法，因此发明人在此不再 进行详细的描述），此处自由振动信号h(t)可表示为 

$h\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{A}_i{e}^{-{\mathrm{\&epsiv;}}_i{\mathrm{\&omega;}}_{\mathrm{ni}}t}\cos ({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)$

式中，P为系统模态阶数；A_i为第i阶系统模态响应信号振幅；ω_ni为第i阶系统无阻尼固有频率；ε_i为第i阶系统阻尼比；θ_i为第i阶模态响应初始相位；ω_di为第i阶系统有阻尼固有频率，且有

因此，当系统阻尼较小时，有ω_di＝ω_ni。 

(4)使用Morlet小波对自由振动信号H(t)进行小波变换，得到自由振动信号H(t)的小波变换系数矩阵WT(a,b)，其中 

Z(t)＝h(t)+jH(h(t)) 

$=\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{A}_i{e}^{-{\mathrm{\&epsiv;}}_i{\mathrm{\&omega;}}_{\mathrm{ni}}t}\cos ({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)+j\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{A}_i{e}^{-{\mathrm{\&epsiv;}}_i{\mathrm{\&omega;}}_{\mathrm{ni}}t}\sin ({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)$

$=\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{A}_i{e}^{-{\mathrm{\&epsiv;}}_i{\mathrm{\&omega;}}_{\mathrm{ni}}t}[\cos ({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)+j\sin ({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)]=\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{A}_i{e}^{-{\mathrm{\&epsiv;}}_i{\mathrm{\&omega;}}_{\mathrm{ni}}t}{e}^{j({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)}$

式中，H(h(t))为h(t)的Hilbert变换（Hilbert变换对于本领域内技术人员来说是熟知的，故本文不再对Hilbert变换进行详细描述）。 

$\mathrm{WT}(a,b)=\&lang;h\left(t\right),{\mathrm{\&psi;}}_{a,b}\left(t\right)\&rang;=\frac{1}{2}\&lang;Z\left(t\right),{\mathrm{\&psi;}}_{a,b}^{*}\left(t\right)\&rang;$

$=\frac{1}{2\sqrt{2\mathrm{a\&pi;}}}\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}{A}_i{e}^{-{\mathrm{\&epsiv;}}_i{\mathrm{\&omega;}}_{\mathrm{ni}}t}{e}^{j({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)}{e}^{(\frac{t-b}{a}{)}^2}{e}^{-j{\mathrm{\&omega;}}_0\left(\frac{t-b}{a}\right)}\mathrm{dt}$

$=\frac{1}{2\sqrt{2\mathrm{a\&pi;}}}\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}{B}_i\left(t\right)e^{({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)}{e}^{(\frac{t-b}{a}{)}^2}{e}^{-j{\mathrm{\&omega;}}_0\left(\frac{t-b}{a}\right)}\mathrm{dt}$

此处，令Bi(t)＝Aie-εiωnit。 

在t＝b处对上式进行Taylor级数展开（Taylor级数对于本领域内技术人员来说是熟知的，故本文不再对Taylor级数进行详细描述），有 

$\mathrm{WT}(a,b)=\frac{1}{2\sqrt{2\mathrm{a\&pi;}}}\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}\{B_i\left(b\right)+O[B_i^{\&prime;}\left(b\right)\left]\right\}{e}^{j({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)}{e}^{(\frac{t-b}{a}{)}^2}{e}^{-j{\mathrm{\&omega;}}_0\left(\frac{t-b}{a}\right)}\mathrm{dt}$

忽略上式中的高阶无穷小项，有 

$\mathrm{WT}(a,b)\&ap;\frac{1}{2\sqrt{2\mathrm{a\&pi;}}}\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}{B}_i\left(b\right)e^{j({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)}{e}^{(\frac{t-b}{a}{)}^2}{e}^{-j{\mathrm{\&omega;}}_0\left(\frac{t-b}{a}\right)}\mathrm{dt}$

$\&ap;\frac{\sqrt{\mathrm{a\&pi;}}}{2}\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{B}_i\left(b\right)e^{j({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)}{e}^{-({\mathrm{a\&omega;}}_{\mathrm{di}}{\mathrm{\&omega;}}_0{)}^{2/2}}\&ap;\frac{\sqrt{\mathrm{a\&pi;}}}{2}\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{A}_i{e}^{-{\mathrm{\&epsiv;}}_i{\mathrm{\&omega;}}_{\mathrm{ni}}b}{e}^{-(a{\mathrm{\&omega;}}_{\mathrm{di}}-{\mathrm{\&omega;}}_0{)}^{2/2}}{e}^{j({\mathrm{\&omega;}}_{\mathrm{di}}b+{\mathrm{\&theta;}}_i)}$

式中，WT(a,b)表示小波变换系数矩阵中第a行和第b列对应的元素，其中小波变换系数矩阵的总行数A对应频率，为512；小波变换系数矩阵的总列数B对应时间，为100；e为自然对数的底数，是一常数，取值为2.71828；j为虚数单位，为

ω₀为Morlet小波的中心频率，即为Morlet小波窗口频率的中心，为常数69.08。 

(5)对小波变换系数矩阵WT(a,b)取模，计算小波变换系数矩阵WT(a,b)的模值|WT(a,b)|，根据小波变换系数矩阵的模值|WT(a,b)|做小波变换时频图：即以时间为横坐标轴，频率为纵坐标轴，以小波变换系数矩阵WT(a,b)的模值为显示结果作图，就可以得到小波变换时频图，其中 

$\left|\mathrm{WT}(a,b)\right|\&ap;\frac{\sqrt{\mathrm{a\&pi;}}}{2}\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{A}_i{e}^{-{\mathrm{\&epsiv;}}_i{\mathrm{\&omega;}}_{\mathrm{ni}}b}{e}^{-({\mathrm{a\&omega;}}_{\mathrm{di}}{\mathrm{\&omega;}}_0{)}^{2/2}}$

(6)采用疯狂爬山法提取小波变换时频图中的小波脊线，小波变换时频图中各条小波脊线的纵坐标为变压器绕组的各阶固有频率： 

6a.定义并初始化度量密度矩阵D，初始值为0，度量密度矩阵D的行列数与小波变换系数矩阵WT(a,b)的行列数相同，即为512行和100列，在度量密度矩阵D所有元素中，随机选取3200个元素作为运动点（运动点的个数为小波变换系数矩阵元素数目总数的四分之一），用于记录相应位置的密度参量；定义系统循环变量为T_t，其初值为T₀＝max(WT)-min(WT，)其中，max(WT)和min(WT)分别为小波变换系数矩阵WT(a,b)的模极大值和模极小值，本实施例中分别为0.2123和4.9813×10^-7； 

6b.在初始时刻t＝t₀＝1，确定运动点X_k对应的位置为(i,j)，即第k个运动点在t₀时刻满足X_k(t₀)＝(i,j)，此处有0＜k≤3200，其中i和j分别为运动点在小波变换时频图中的位置坐标，且1≤i≤512和1≤j≤100； 

6c.在下一时刻即t+1时刻， 

在水平方向上， 

若满足1＜i＜A，则以概率P＝50%向左或向右移动运动点X_k(t₀)到位置(i,j′)，即有 

j′＝j+1或j′＝j-1 

若满足i＝1，则向右移动运动点X_k(t₀)到位置(i,j′)，j′＝j+1； 

若满足i＝B，则向左移动运动点X_k(t₀)到位置(i,j′)，j′＝j-1； 

在垂直方向上， 

若满足j′＝1，则向下移动初始运动点X_k(t₀)到位置(i′,j′)，i′＝i+1； 

若满足j′＝A，则向上移动初始运动点X_k(t₀)到位置(i′,j′)，i′＝i-1； 

若满足1＜j′＜B，先以概率P＝50%的概率向上或向下移动运动点X_k(t₀)到位置(i_p,j′)，即有 

i_p＝i+1或i_p＝i-1 

若满足WT(i_p,j′)＞WT(i,j′)，则运动点X_k(t₀)在当前时刻的位置为(i′,j′)，i′＝i_p，转至步骤6d； 

若满足WT(i_p,j′)＜WT(i,j′)，根据下式计算垂直方向上的移动概率p： 

$p=\exp \left[\frac{\mathrm{wt}(i_p,j^{\&prime;})-\mathrm{wT}(i,j^{\&prime;})}{T_t}\right]$

如上式所定义的，T_t为系统循环变量； 

将运动点X_k(t₀)从位置(i_p,j′)按照概率p保持不动，即i′＝i_p；按照概率1-p从位置(i_p,j′)返回到位置(i,j′)，即i′＝i。 

6d.记录该运动点X_k在t+1时刻的位置，即为(i′,j′)，根据下式更新度量密度矩阵，即有 

D_t+1(i′,j′)＝D_t(i′,j′)+WT(i′,j′) 

6e.对所有n1个点重复步骤6a～6d； 

6f.计算系统当前循环参量为T_t＝T_t/(t+1)²，对所有n1个点重复步骤6a～6e，直到T_t＜T₀/1000； 

6g.确定最后的密度度量矩阵D′(i,j)，其中D′(i,j)中的各个元素为脊点： 

$D^{\&prime;}(i,j)=\left\{\begin{array}{cc}D(i,j)& D(i,j)\&GreaterEqual;T_h\\ 0& D(i,j)<T_h\end{array}\right.$

式中，T_h为设定的密度阈值，本实施例中为0.05； 

6h.根据脊点做出小波脊线，即以时间为横坐标轴，频率为纵坐标轴，最终密度度量矩阵D′(i,j)各个元素的模值为显示结果作图，得到小波脊线图。 

(7)计算小波系数矩阵WT(a,b)中第a行各个元素的自然对数值，此处，a的取值为步骤(6)中小波脊线的纵坐标，计算小波系数矩阵WT(a,b)中第a行 各个元素所在直线对时间的斜率λ，则可计算得到步骤(6)中各阶固有频率对应的阻尼比ε，其计算公式为 

ε_k＝-λ/(2πf_k) 

式中，f_k为第k阶固有频率；ε_k为第k阶固有频率对应的阻尼比。此处，k的取值为4。该变压器绕组的各阶固有频率及其阻尼比，表中同时给出了使用PolyMax法得到的变压器绕组固有频率和阻尼比。 

(8)根据小波变换系数矩阵WT(a,b)识别步骤(6)中各个固有频率对应的振型，计算过程为： 

(8a).选取小波变换系数矩阵WT(a,b)的模值最小的点作为参考点，记为r点； 

(8b).使用最小二乘法计算r点对应的第m阶小波变换系数WT_r(a,b)的最佳常数解，此处，m＝1,2,3,4。最小二乘法是本领域内技术人员来说是熟知的，故本文不再对最小二乘法进行详细描述。 

(8c).使用最小二乘法计算第p点对应的第m阶小波变换系数WT_p(a,b)的最佳常数解，此处，p点指除了r点之外的任一点，为p＝1,…,19。则可得第p个响应点对参考点r的最佳最小二乘解即振动信息，其计算公式为 

${\mathrm{\&phi;}}_k^j=\frac{{\mathrm{WT}}_r^T{\mathrm{WT}}_p}{{\mathrm{WT}}_r^T{\mathrm{WT}}_r}$

式中，

为r点的小波系数变换矩阵的转置；WT_p为p点的小波系数变换矩阵；WT_r为r点的小波系数变换矩阵。 

依次计算各个相应点对参考点r的最佳最小二乘解，即可得第k阶固有频率对应的振型。 

要注意的是，以上列举的仅为本发明的具体实施例，显然本发明不限于以上实施例，随之有着许多的类似变化。本领域的技术人员如果从本发明公开的内容直接导出或联想到的所有内容，均应属于本发明的保护范围。 

Claims (4)

1.一种变压器绕组模态参数识别方法，其特征在于，该方法包括下列步骤：

（1）在变压器绕组表面设置N个测点，将N个振动加速度传感器分别对应放置在各测点，在变压器绕组上放置激振器，并向该激振器输入信号V_i对变压器绕组进行激振，所述的N个振动加速度传感器分别采集各个测点的振动信号V_oi，i=1，2，……，N；

（2）分别对输入信号V_i和各个测点的振动信号V_oi进行傅里叶变换，得到变压器绕组的振动频响曲线H(ω)：

$H\left(\mathrm{\&omega;}\right)=\left(\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{V}_{\mathrm{oi}}\left(\mathrm{\&omega;}\right)\right)/V_i\left(\mathrm{\&omega;}\right)$

式中，V_i(ω)为输入信号V_i的傅里叶变换；V_oi(ω)为各个测点振动信号V_oi的傅里叶变换；N为振动加速度传感器的个数；

（3）对所述的振动频响曲线H(ω)做傅里叶反变换，得到变压器绕组的自由振动信号h(t)：

$h\left(t\right)=\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{A}_i{e}^{-{\mathrm{\&epsiv;}}_i{\mathrm{\&omega;}}_{\mathrm{ni}}t}\cos ({\mathrm{\&omega;}}_{\mathrm{di}}t+{\mathrm{\&theta;}}_i)$

式中，P为系统模态阶数；A_i为第i阶系统模态响应信号振幅；ω_ni为第i阶系统无阻尼固有频率；ε_i为第i阶系统阻尼比；θ_i为第i阶模态响应初始相位；ω_di为第i阶系统有阻尼固有频率，且有

因此，当系统阻尼较小时，有ω_di＝ω_ni；

（4）利用Morlet小波对自由振动信号h(t)进行小波变换，得到自由振动信号h(t)的小波变换系数矩阵WT(a,b)：

$\mathrm{WT}(a,b)\&ap;\frac{\sqrt{\mathrm{a\&pi;}}}{2}\underset{i=1}{\overset{P}{\mathrm{\&Sigma;}}}{A}_i{e}^{-{\mathrm{\&epsiv;}}_i{\mathrm{\&omega;}}_{\mathrm{ni}}b}{e}^{-{(a{\mathrm{\&omega;}}_{\mathrm{di}}-{\mathrm{\&omega;}}_0)}^2/2}{e}^{j({\mathrm{\&omega;}}_{\mathrm{di}}b+{\mathrm{\&theta;}}_i)}$

式中，ω₀为Morlet小波的中心频率；a为频率尺度；b为时间。

（5）对小波变换系数矩阵WT(a,b)取模，计算小波变换系数矩阵WT(a,b)的模值|WT(a,b)，根据小波变换系数矩阵WT(a,b)的模值|WT(a,b)做小波变换时频图；

（6）提取小波变换时频图中的小波脊线，各条小波脊线的纵坐标为变压器绕组的各阶固有频率；

（7）计算小波系数模极大值max(WT)对时间t的导数λ，根据该导数λ，计算步骤(6)得到的各阶固有频率分别所对应的阻尼比ε_k，公式为：

ε_k＝-λ/(2πf_k)

式中，f_k为第k阶固有频率；ε_k为第k阶固有频率所对应的阻尼比；

（8）根据小波变换系数矩阵WT(a,b)，识别步骤(6)得到的各个固有频率分别所对应的振型。

2.如权利要求1所述的变压器绕组模态参数识别方法，其特征在于，所述的向激振器输入信号V_i为白噪声信号。

3.如权利要求1所述的变压器绕组模态参数识别方法，其特征在于，所述步骤（6）中提取小波变换时频图中的小波脊线，具体方法如下：

6a.定义并初始化度量密度矩阵D，初始值为0，度量密度矩阵D的行列数与小波变换系数矩阵WT(a,b)的行列数相同，在度量密度矩阵D所有元素中，随机选取n1个元素作为运动点，用于记录相应位置的密度参量，此处，n1为小波变换系数矩阵WT(a,b)元素总个数的四分之一；定义系统循环变量为T_t，其初值为T₀＝max(WT)-min(WT)，其中max(WT)和min(WT)分别为小波变换系数矩阵WT(a,b)的模极大值和模极小值；

6b.在初始时刻t＝t₀，确定运动点X_k对应的位置为(i,j)，即第k个运动点在t₀时刻满足X_k(t₀)＝(i,j)，此处有0＜k≤n1，i为第k个运动点在小波变换时频图中的横坐标位置，j为第k个运动点在小波变换时频图中的纵坐标位置，且1≤i≤A和1≤j≤B；A为小波变换系数矩阵WT(a,b)的行数，B为小波变换系数矩阵WT(a,b)的列数；

6c.在下一时刻即t+1时刻，

在水平方向上，

若满足1＜i＜A，则以概率P＝50%向左或向右移动运动点X_k(t₀)到位置(i,j')，即有

j'＝j+1或j'＝j-1

若满足i＝1，则向右移动运动点X_k(t₀)，该运动点在当前时刻t+1的位置为(i,j')，j'＝j+1；

若满足i＝B，则向左移动运动点X_k(t₀)，该运动点在当前时刻t+1的位置为位置(i,j')，j'＝j-1；

在垂直方向上，

若满足j'＝1，则向下移动运动点X_k(t₀)，该运动点在当前时刻t+1的位置为(i',j')，i'＝i+1；

若满足j'＝A，则向上移动运动点X_k(t₀)，该运动点在当前时刻t+1的位置为(i',j')，i'＝i-1；

若满足1＜j'＜B，先以概率P＝50%的概率向上或向下移动运动点X_k(t₀)，该运动点在当前时刻t+1的位置为(i_p,j')，即有

i_p＝i+1或i_p＝i-1

若满足WT(i_p,j')＞WT(i,j')，则运动点X_k(t₀)在当前时刻的位置为(i',j')，i'＝i_p，转至步骤6d；

若满足WT(i_p,j')＜WT(i,j')，根据下式计算垂直方向上的移动概率p：

$p=\exp \left[\frac{\mathrm{WT}(i_p,j^{\&prime;})-\mathrm{WT}(i,j^{\&prime;})}{T_t}\right]$

如上式所定义的，T_t为系统循环变量；

将运动点X_k(t₀)从位置(i_p,j')按照概率p保持不动，即i'＝i_p；按照概率1-p从位置(i_p,j')返回到位置(i,j')，即i'＝i。

6d.记录该运动点X_k在t+1时刻的位置，即为(i',j')，根据下式更新度量密度矩阵，即有

D_t+1(i',j')＝D_t(i',j')+WT(i',j')

6e.对所有n1个点重复步骤6a～6d；

6f.计算系统当前循环参量T_t＝T_t/(t+1)²，对所有n1个点重复步骤6a～6e，直到T_t＜T₀/1000；

6g.确定最后的密度度量矩阵D'(i,j)，其中D'(i,j)中的各个元素为脊点：

$D^{\&prime;}(i,j)=\left\{\begin{array}{cc}D(i,j)& D(i,j)\&GreaterEqual;T_h\\ 0& D(i,j)<T_h\end{array}\right.$

式中，T_h为设定的密度阈值；

6h.根据脊点做出小波脊线。

4.如权利要求1所述的变压器绕组模态参数识别方法，其特征在于，所述的步骤（8）根据小波变换系数矩阵WT(a,b)，识别步骤(6)得到的各个固有频率分别所对应的振型，具体方法如下：

8a.选取小波变换系数矩阵WT(a,b)的模值最小的点作为参考点，记为r点；

8b.利用最小二乘法计算r点对应的第m阶小波变换系数WT^r(a,b)的最佳常数解；

8c.利用最小二乘法计算第p个响应点对应的第m阶小波变换系数WT^p(a,b)的最佳常数解，计算得到第p个响应点对r点的最佳最小二乘解即振动信息，公式为

${\mathrm{\&phi;}}_k^j=\frac{{{\mathrm{WT}}_r}^T{\mathrm{WT}}_p}{{{\mathrm{WT}}_r}^T{\mathrm{WT}}_r}$

式中，WT_r ^T为r点的小波系数变换矩阵的转置；WT_p为第p个响应点的小波系数变换矩阵；WT_r为r点的小波系数变换矩阵。

依次计算各个响应点对参考点的最佳最小二乘解，即可得第k阶固有频率对应的振型。