CN103822696A - 一种变压器绕组状态的诊断方法 - Google Patents

Abstract

本发明公开了一种变压器绕组状态诊断方法，其包括下列步骤：（1）采集变压器绕组各个测点的振动信号；（2）对振动信号进行傅里叶变换，得到变压器绕组的振动频响曲线；（3）振动频响曲线做傅里叶反变换，得到变压器绕组的自由振动信号；（4）使用Morlet小波对自由振动信号进行小波变换；（5）做小波变换时频图；（6）提取小波变换时频图中的小波脊线，各条小波脊线的纵坐标为变压器绕组的各阶固有频率；（7）根据固有频率的变化对变压器绕组状态进行判别。

Description

一种变压器绕组状态的诊断方法

技术领域

本发明涉及一种信号监测方法，尤其涉及一种变压器绕组状态的诊断方法。

背景技术

变压器是电力系统中最重要的设备之一，其运行的稳定性对电力系统安全影响重大。随着我国电网容量的日益增大，短路容量亦随之不断增大，变压器出口短路形成的冲击电流产生的巨大电磁作用力对变压器绕组的机械强度和动稳定性构成了严重的威胁。目前变电站设备及线路的运行环境始终不容乐观，因外部短路造成变压器绕组受冲击而引发的变形，是变压器运行过程中较为常见的故障，其对系统的安全运行造成了很大的威胁。

变压器遭受突发短路后，其绕组可能首先发生松动或轻微变形，通过大量的实验研究分析变压器绕组变形具有累积效应，如果对于松动或变形不能及时发现和修复，那么在变压器的松动或变形累积到一定程度后会使变压器的抗短路能力大幅下降而在遭受较小的冲击电流下也会引发大的事故发生。

绕组的变形一方面会导致机械抗短路电流冲击能力的下降，另一方面也会导致线圈内部局部绝缘距离发生变化，使局部出现绝缘薄弱点，当遇到过电压作用时，绕组有可能发生饼间或匝间短路导致变压器绝缘击穿事故，或者由于局部场强增大而引起局部放电，绝缘损伤部位会逐渐扩大，最终导致变压器发生绝缘击穿事故而引发进一步的事态扩大。

因此，在运行过程中当变压器经历了外部短路事故后或运行一段时间后的常规检修中，如何有效地检测出变压器绕组是否存在松动和变形，从而判断变压器是否需要检修处理显得十分重要，是保障变压器安全运行的一个重要手段，因此变压器绕组变形的检测是目前变压器常规试验项目之一。

目前实际应用的对变压器绕组状态的检测方法主要有以下三种：

1、短路阻抗法

变压器短路阻抗是当负载阻抗为零时变压器内部的等效阻抗，短路阻抗是变压器绕组的漏抗和电阻的矢量和，由于变压器直流电阻相对于漏抗数值很小，因此变压器的短路阻抗反映的主要是变压器绕组的漏抗。由变压器的理论分析可知，变压器漏抗值是由绕组的几何尺寸所决定的，或者说是由绕组的结构决定的，一旦变压器绕组发生变形，从理论上来说变压器的漏抗相应也会发生变化，因此通过对变压器短路阻抗的检测可以间接地反映变压器绕组内部是否发生了变形。

一般情况下，运行中的变压器受到了短路电流的冲击后，或在定期常规检查时要将测得的短路阻抗值与原有的记录进行比较来判断绕组是否发生了变形，如果短路阻抗值变化较大，例如国标中设定为变化超过3%，则可确认绕组有显著变形。

按照有关标准规定，变压器在短路阻抗测试试验中，要求测量每一相的短路阻抗，并把试验后所测量的短路阻抗值与以往试验的数据加以比较，根据其变化的程度，作为判断被试变压器绕组是否合格的重要依据之一。

从实际应用情况来看，短路阻抗法在长期的生产实践中已建立了标准，判据较为明确，在国际电工标准IEC60076-5和GB1095-85中均明确给出了线圈变形程度的判据。但很多情况下这种方法的灵敏度很低，故障的检出率较低，只有在线圈整体变形情况较为严重时才能够得到较明确的反映。

2、频响分析法

频响分析法的基本原理是将变压器绕组视为一个分布参数网络，它由对地电容C、纵向电容K、电感L等分布参数构成一无源线性双端口网络，该网络的特性在频域上可以用传递函数H(jω)来描述。

绕组发生局部机械变形后，其内部的分布电感L、纵向电容K和对地电容C等分布参数会发生相应的变化，从而在网络的传递函数H(jω)上得到反映。因此分析变压器绕组的网络传递函数曲线的变化情况就可以分析内部的网络电参数是否发生变化，从而推断相应的机械结构是否发生了变形，这是频响分析法测试变压器绕组变形的依据和基础。

频响法测试首先将一稳定的正弦扫频电压信号V_i施加到被试变压器绕组的一端，然后同时记录该端口V_i和其它输出端口上的电压V_o,从而得到该被试绕组的一组频响特性曲线，其表达式为

H(jω)＝V_o/V_i

频响法的测试灵敏度较短路阻抗法高，但由于其频响波形的复杂性，对绕组状况的判别需要较多的经验，较难形成明确的定量判据，因此至今没有形成判别标准。

上述两种方法是目前判别变压器绕组状况最常用的，两种方法都是采用电测方法，出发点都是基于变压器绕组发生明显变形的状况下模型中对应的元件电参数发生变化来进行测量判别，这对变压器绕组发生较明显的变形情况较为适宜，但对绕组发生轻微变形，尤其是对变压器绕组存在的相对松动和扭曲变形的状态不能给出较明确的判断，因为这些情况下反映在等效电路模型中的电参数几乎没有变化，其传递函数的变化也就非常小。然而变压器绕组松动或扭曲变形对其抗短路能力有很大的影响，因此研究绕组的状况需有灵敏度更高的方法来进行判别。

3、振动分析法

振动分析法的基本原理是把变压器绕组看作一个机械结构体，则当绕组结构或受力发生任何变化时，都可以从它的机械振动特性变化上得到反映。因此，可通过分析箱壁上的振动信号来对绕组的工作状态进行检测。与前述电气测量法相比较，振动分析法的最大优点是可通过吸附在变压器箱壁上的振动传感器来获得变压器的振动信号，通过分析其振动特性的变化来判断绕组状态的变化情况，只要绕组的机械特性（如结构变形、预紧力松动等）发生变化，都可以从它的机械振动特性变化上得到反映，从而大大提高了检测的灵敏度。此外，将振动传感器置于箱壁上的振动检测与整个强电系统没有直接的连接，对于整个电气系统的正常运行没有任何影响，因此，可发展成为一种较准确、便捷、安全的在线监测方法。

发明内容

本发明的目的是提供一种变压器绕组状态的诊断方法，该方法利用利用振动模态特性对变压器绕组的工作状态进行判别。

为了实现上述发明目的，本发明提供了一种变压器绕组状态诊断方法，其包括下列步骤：

（1）在变压器绕组上放置激振器，向激振器输入信号V_i对变压器绕组进行激振，采用放置在变压器绕组表面的N个振动加速度传感器采集和记录各个测点的振动信号V_oi(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)（傅里叶反变换是本领域内常用的数学方法，因此发明人在此不再进行详细的描述）。

(4)使用Morlet小波对自由振动信号H(t)进行小波变换，得到自由振动信号H(t)的小波变换系数矩阵WT(m,k)，然后对小波变换系数矩阵WT(m,k)取模，计算小波变换系数矩阵WT(m,k)的模值，其中：

$\mathrm{WT}(m,k)\≈\frac{\sqrt{\mathrm{m\π}}}{2}{\mathrm{Ae}}^{-\mathrm{\ϵ}{\mathrm{\ω}}_{n}k}{e}^{-{({\mathrm{m\ω}}_d-{\mathrm{\ω}}_0)}^2/2}{e}^{j({\mathrm{\ω}}_{d}k+\mathrm{\θ})}$

式中，WT(m,k)表示小波变换系数矩阵中第m行和第k列对应的元素，其中小波变换系数矩阵的总行数M对应频率，小波变换系数矩阵的总列数K对应时间；A为自由振动信号H(t)的最大幅值；ε为阻尼；ω_n和ω_d分别为振动信号的无阻尼固有频率和有阻尼固有频率；θ为系统初始相位，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(m,k)的模值做小波变换时频图：即以时间为横坐标轴，频率为纵坐标轴，以小波变换系数矩阵WT(m,k)的模值为显示结果作图，就可以得到小波变换时频图；

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

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

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

6b.在初始时刻t＝t₀，确定运动点X_t对应的位置为(i,j)，即有X_k(t+1)＝(i,j)，其中i和j分别为运动点在小波变换时频图中的位置坐标，且1≤i≤M和1≤j≤K；

6c.在下一时刻t＝t₀+1，若满足1＜i＜m和1＜j＜k，则首先在水平方向以概率P＝50%向左或向右移动运动点X_t(i,j)，即有

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

在垂直方向上，若满足WT(i′，j′)＞WT(i,j′)，则以概率P＝50%的概率向上或向下移动运动点X_t(i,j)，即有

i'＝i+1或i'＝i-1

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

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

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

按照移动概率p进行垂直方向的移动，即X_k(t+1)＝(i′,j′)，按照概率1-p垂直方向不动，即X_k(t+1)＝(i,j')，X_k(t+1)为第k个运动点在t+1时刻的位置，且有1≤k≤n；

若满足i＝1和1＜j＜k、i＝m和1＜j＜k，则以概率P＝50%向左或向右移动运动点X_t(i,j)；

若满足1＜i＜m和j＝1、1＜i＜m和j＝k，则以概率P＝50%向上或向下移动初始运动点X_t(i,j)。

6d.记录新的运动点位置为X_t+1(i',j′)，并根据下式更新度量密度矩阵，即有

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

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

6f.更新系统时间步长t＝t+1和系统当前循环参量T_t＝T_t-1/t²，对所有n个点重复步骤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)根据固有频率的变化对变压器绕组状态进行判别：当变压器绕组的各阶固有频率向低频方向偏移且数值减小为原来的5%及以上时，判定变压器绕组发生松动或变形，此时需要及时进行处理，避免形成重大故障。

进一步地，在本发明所述的变压器绕组状态的诊断方法中，向激振器输入的信号Vi为白噪声信号。

本发明所述的变压器绕组状态的诊断方法通过测试变压器绕组的振动频响曲线和识别变压器绕组的固有频率特性而对变压器绕组的工作状态进行判别，从而准确、高效地判断变压器绕组的工作状态，以便于能够及时发现问题，对变压器进行及时检修。

附图说明

图1显示了本技术方案的实施例中变压器绕组状态良好时的振动频响曲线。

图2显示了本技术方案的实施例中变压器绕组状态恶化时的振动频响曲线。

图3显示了本技术方案的实施例中变压器绕组状态良好时的小波脊线图。

图4显示了本技术方案的实施例中变压器绕组状态恶化时的小波脊线图。

具体实施方式

下面将结合说明书附图和具体实施例对本技术方案所述的变压器绕组状态诊断方法做进一步的详细说明。

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

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

(2)分别对输入的白噪声信号V_i和各个测点的振动信号V_oi进行傅里叶变换，得到变压器绕组的振动频响曲线H(ω)：（图1显示了该变压器绕组状态良好时的振动频响曲线，图2显示了该变压器绕组状态恶化时的振动频响曲线）

$H\left(\mathrm{\&omega;}\right)=\left(\underset{i=11}{\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)。

(4)使用Morlet小波对自由振动信号H(t)进行小波变换，得到自由振动信号H(t)的小波变换系数矩阵WT(m,k)，然后对小波变换系数矩阵WT(m,k)取模，计算小波变换系数矩阵WT(m,k)的模值，其中：

$\mathrm{WT}(m,k)\&ap;\frac{\sqrt{\mathrm{m\&pi;}}}{2}{\mathrm{Ae}}^{-\mathrm{\&epsiv;}{\mathrm{\&omega;}}_{n}k}{e}^{-{({\mathrm{m\&omega;}}_d-{\mathrm{\&omega;}}_0)}^2/2}{e}^{j({\mathrm{\&omega;}}_{d}k+\mathrm{\&theta;})}$

式中，WT(m,k)表示小波变换系数矩阵中第m行和第k列对应的元素，其中小波变换系数矩阵的总行数M对应频率，为512；小波变换系数矩阵的总列数K对应时间，为100；A为自由振动信号H(t)的最大幅值；ε为阻尼；ω_n和ω_d分别为振动信号的无阻尼固有频率和有阻尼固有频率；θ为系统初始相位，e为自然对数的底数，是一常数，取值为2.71828；j为虚数单位，为

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

（5）以时间为横坐标轴，频率为纵坐标轴，以小波变换系数矩阵WT(m,k)的模值为显示结果作图，得到小波变换时频图；

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

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

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

6c.在下一时刻t＝t₀+1，若满足1＜i＜m和1＜j＜k，则首先在水平方向以概率P＝50%向左或向右移动运动点X_t(i,j)，即有

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

在垂直方向上，若满足WT(i′，j′)＞WT(i,j′)，则以概率P＝50%的概率向上或向下移动运动点X_t(i,j)，即有

i'＝i+1或i'＝i-1

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

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

按照移动概率p进行垂直方向的移动，即X_k(t+1)＝(i′,j′)，按照概率1-p垂直方向不动，即X_k(t+1)＝(i,j')，X_k(t+1)为第k个运动点在t+1时刻的位置，且有1≤k≤n；

若满足i＝1和1＜j＜k、i＝m和1＜j＜k，则以概率P＝50%向左或向右移动运动点X_t(i,j)；

若满足1＜i＜m和j＝1、1＜i＜m和j＝k，则以概率P＝50%向上或向下移动初始运动点X_t(i,j)。

6d.记录新的运动点位置为X_t+1(i',j′)，并根据下式更新度量密度矩阵，即有

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

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

6f.更新系统时间步长t＝t+1和系统当前循环参量T_t＝T_t-1/t²，对所有n个点重复步骤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)各个元素的模值为显示结果作图，得到小波脊线。图3显示了该变压器绕组状态良好时的小波脊线图。图4显示了该变压器绕组状态恶化时的小波脊线图。

(7)根据固有频率的变化对变压器绕组状态进行判别：当变压器绕组的各阶固有频率向低频方向偏移且数值减小为原来的5%及以上时，判定变压器绕组发生松动或变形，此时需要及时进行处理，避免形成重大故障。

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

Claims (3)

1.一种变压器绕组状态的诊断方法，其特征在于，包括下列步骤：

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

（2）分别对输入的信号Vi和各个测点的振动信号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(ω)为输入的信号Vi的傅里叶变换；V_oi(ω)为各个测点振动信号V_oi的傅里叶变换；N为振动加速度传感器的个数；

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

（4）使用Morlet小波对自由振动信号H(t)进行小波变换，得到自由振动信号H(t)的小波变换系数矩阵WT(m,k)，然后对小波变换系数矩阵WT(m,k)取模，计算小波变换系数矩阵WT(m,k)的模值，其中：

$\mathrm{WT}(m,k)\&ap;\frac{\sqrt{\mathrm{m\&pi;}}}{2}{\mathrm{Ae}}^{-\mathrm{\&epsiv;}{\mathrm{\&omega;}}_{n}k}{e}^{-{({\mathrm{m\&omega;}}_d-{\mathrm{\&omega;}}_0)}^2/2}{e}^{j({\mathrm{\&omega;}}_{d}k+\mathrm{\&theta;})}$

式中，WT(m,k)表示小波变换系数矩阵中第m行和第k列对应的元素，其中小波变换系数矩阵的总行数M对应频率，小波变换系数矩阵的总列数K对应时间；A为自由振动信号H(t)的最大幅值；ε为阻尼；ω_n和ω_d分别为振动信号的无阻尼固有频率和有阻尼固有频率；θ为系统初始相位，e为自然对数的底数，是一常数，取值为2.71828；j为虚数单位，为

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

（5）根据小波变换系数矩阵WT(m,k)的模值做小波变换时频图；

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

（7）根据固有频率的变化对变压器绕组状态进行判别：当变压器绕组的各阶固有频率向低频方向偏移且数值减小为原来的5%及以上时，判定变压器绕组发生松动或变形。

2.如权利要求1所述的变压器绕组状态的诊断方法，其特征在于，所述向激振器输入的信号Vi为白噪声信号。

3.如权利要求1所述的变压器绕组状态的诊断方法，其特征在于，按照下述步骤采用疯狂爬山法提取小波变换时频图中的小波脊线：

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

6b.在初始时刻t＝t₀，确定运动点X_t对应的位置为(i,j)，即有X_k(t+1)＝(i,j)，其中i和j分别为运动点在小波变换时频图中的位置坐标，且1≤i≤M和1≤j≤K；

6c.在下一时刻t＝t₀+1，若满足1＜i＜m和1＜j＜k，则在水平方向以概率P＝50%向左或向右移动运动点X_t(i,j)，即有

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

在垂直方向上，若满足WT(i',j')＞WT(i,j')，则以概率P＝50%的概率向上或向下移动运动点X_t(i,j)，即有

i'＝i+1或i'＝i-1

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

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

将运动点按照移动概率p进行垂直方向的移动，即X_k(t+1)＝(i',j′)，按照概率1-p垂直方向不动，即X_k(t+1)＝(i,j')；

若满足i＝1和1＜j＜k、i＝m和1＜j＜k，则以概率P＝50%向左或向右移动运动点X_t(i，j)；

若满足1＜i＜m和j＝1、1＜i＜m和j＝k，则以概率P＝50%向上或向下移动初始运动点X_t(i,j)；

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

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

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

6f.更新系统时间步长t＝t+1和系统当前循环参量T_t＝T_t-1/t²，对所有n个点重复步骤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.根据脊点做出小波脊线。

Images