一种超声骨刀的谐振频率跟踪方法与流程

1.本发明涉及超声换能器技术领域，具体是一种超声骨刀的谐振频率跟踪方法。


背景技术：

2.超声换能器是一种能把高频电能转化为机械能的装置，利用材料的压电效应将电信号转换为机械振动，是一种能量转换器件，由于其体积小、效率高，广泛应用于医疗、工业加工及军事等领域。超声换能器自身有一个固定的谐振频率，只有当换能器在谐振频率下工作时，才能保证输出功率最大、效率最高，达到最佳的功率输出效果；但在使用过程中随着温度以及环境等变化，谐振频率会发生偏移，超声换能器工作在非谐振频率时会出现能量转化效率低、手柄发热、输出功率不可控等问题，进而造成无功消耗大、换能器寿命缩短，甚至影响安全运行。为了使超声换能器工作在谐振频率，需要对谐振频率进行实时跟踪，使驱动频率自动准确且快速地调节至换能器的谐振频率。在现有技术中常采用的有以下几类频率跟踪方法：
3.1.阻抗圆/导纳圆特征点选取法：
4.如中国专利公开了申请号为202010791986.4名称为一种超声骨刀动态阻抗反馈方法，其包括电压霍尔传感器、电流霍尔传感器、幅相检测及低通滤波电路，该专利通过在工作中不断计算阻抗值，采用最小二乘法拟合出超声换能器r-x模型的阻抗圆图，获取得到实际工况下的谐振点，并将驱动频率通过pid算法快速调节至新谐振点，由此实现频率跟踪。虽然该专利在一定程度上能够找到新的谐振点，但是拟合出阻抗圆图需要足够多的样本点，也就是需要在足够宽范围的驱动频率下获取特征阻抗值，实际上工作过程中驱动频率不适合发生大范围的变动。
5.另一专利名称为超声换能器机械谐振频率的快速追踪装置和方法(申请号：201910147405.0)，该专利采用电流传感器和电压传感器获取换能器两端的电流和电压，接着将两路信号经过滤波放大电路，计算出此时换能器地复导纳大小，采样三次，三个样本值确定一个阻抗圆，从而解算出换能器的谐振频率。虽然该专利能够在一定程度上能快速地找到谐振点，但是具有局限性，要求换能器的导纳特性曲线是单一的，就是只能有一个波峰和一个波谷，否则就无法确定三个样本值同时处于同一导纳圆上，而不在同一导纳圆上的点就没有关联性。
6.2.锁相法：
7.如专利一种数字式超声波发生器及其自动锁频方法(申请号：201610611587.9)，其采用获取换能器负载两端的电压波形和电流波形，将两者波形进行乘法计算，得到的直流波形的大小来反应电压电流的波形相位差，认为相位差越小就越接近最佳谐振点，零相位时即处于谐振状态。虽然该专利在一定程度上能够找到谐振点，但是需要有个前提是调谐匹配后换能器的谐振点是零相位的，由于加工工艺过程的差异或者随着长时间使用发生老化，在某些工况下阻抗存在整体的偏容性或偏感性，这样一来零相位的频率点就不存在，因此并不能找到谐振工作点，此方法可靠性较低。
8.3.回路电流检测法：
9.如发明专利名称为基于直流电流确定谐振频率的超声换能器及其方法(申请号：201910149126.8)，该专利提出将整流电路输出的直流电输入至buck电路，buck电路输出的直流电经逆变器处理并转换成交流电，并将交流电调谐处理后输出至超声换能器中进行驱动；其中控制单元将频率不变占空比可变的信号输出至buck电路；将频率和占空比均可变的信号作为逆变器的pwm控制信号；在此过程中采集整流电路输出的直流电的电流值，找到最大的电流值对应的pwm控制信号的频率，即为超声换能器的谐振频率。虽然该专利在一定程度上能够找到谐振点，谐振状态时虽然阻抗值相对较低，但阻抗越低电流越大的点并不一定是谐振点，例如在阻抗特性曲线的极感性和极容性点也会出现最低阻抗，而工作在这两个点上并不能有很好的工作效率。
10.以上各类方法虽然该专利在一定程度上能够找到谐振点，但都具有其局限性，如何快速且准确地确定超声换能器的最佳工作频率，是当前超声技术领域中需要解决的问题。
11.鉴于以上所述，目前现有的技术和方法对超声骨刀谐振频率的快速选取和跟踪仍存在误差较大的缺陷，针对现有存在的技术缺点，本发明提供一种快速而准确的超声骨刀谐振频率跟踪方法。


技术实现要素：

12.本发明的目的在于解决超声骨刀谐振频率准确选取和快速跟踪的问题，基于卡尔曼滤波算法，提出一种超声换能器谐波分量估计方法来完成对谐振频率的进一步选取和与跟踪，能够显著提高超声骨刀谐振频率的选取精度和跟踪速度。
13.为解决上述技术问题，本发明采用如下技术方案；
14.一种超声骨刀的谐振频率跟踪方法，所述方案包括：
15.步骤s1：建立模态频率附近的线性系统模型。设定线性系统的模型为：
[0016][0017]
式中，r为第r次状态向量，为k时刻的状态向量，为k时刻的观测向量，为k-1时刻的系统传输矩阵，为k时刻的量测矩阵，和分别为过程噪声和量测噪声，且二者为互不相干的高斯白噪声，均值和协方差矩阵分别满足：
[0018][0019]
式中，q为状态噪声均值，r为观测噪声均值；
[0020]
步骤s2，基尔霍夫定律离散并等效替代线性系统模式。状态向量分别取
[0021]
式为k时刻第r次谐波分量的幅值，fs为采样频率，ω为激励频率，为相位角；
[0022]
下式中，利用超声换能器在模态频率附近的机电等效电路模型，建立给予频率为
ω的方波激励信号下超声换能器的端部信号电学状态模型，可得：
[0023][0024]
式中，m为谐波最高次数，利用电学状态模型输出提供第r次谐波分量的幅值，观测向量取其中系统的初始量与状态噪声q和观测噪声r均互不相关，且初始状态服从高斯正态分布，其先验均值和协方差矩阵分别为：
[0025][0026]
步骤s3，引入修正因子的抗扰模型算法。提取模态频率附近的机电模型输出的真有功功率因数值(pf)，真有功功率因数值其中p为有功功率，q为无功功率，d为失真无功功率，β在一定程度上反应了超声换能器的工作状态，能够作为抗扰模型的修正因子。
[0027]
首先，建立pf预测模型：
[0028]ypf(i+1)
＝k(i+1)+b
[0029]
[k,b]
t
＝r(i,h)
[0030]
其中，y
pf(i+1)
为超声换能器的pf值与上一时刻时刻pf值的差值在i+1时刻的预测值；r(
·
,
·
)为最小二乘函数；向量i＝(i,i-1,
…
,i-3)
t
,i≥4表示连续量测四次的次数向量；向量h＝(y
pf(i)
,y
pf(i-1)
,
…
,y
pf(i-3)
)
t
,i≥4表示量测四次，当y
pf(i+1)
≥
±
0.1时，β取1；否则β取0；
[0031]
对y
pf(i)
，有：
[0032]ypf(i)
＝δpf
i-δpf
i-1
[0033]
式中，δpfi为i时刻超声换能器的pf值，δpf
i-1
为i-1时刻的pf值；
[0034]
最后，建立pf抗扰模型：
[0035][0036]
pf抗扰模型算法提高了谐波分量幅值最优估计卡尔曼滤波器的鲁棒性，将其数据代入步骤s2的等效替代线性系统模式。
[0037]
步骤s4，根据不同阶次谐波分量幅值的最优估计值计算k时刻的总谐波畸变率(thd)k，
[0038][0039]
式中，为基波分量幅值的最优估计值。
[0040]
最后，建立损失函数。根据实际，pf值的理想值为1，(thd)k的理想值为0，设在k时
刻频率为f时的pf值为(pf)
k,f
，总谐波畸变率为(thd)
k,f
，损失值为l
k,f
，有：
[0041]
l
k,f
＝w
(pf)
[1-(pf)
k,f
]2+w
(thd)
[0-(thd)
k,f
]2[0042]
式中w
(pf)
为(pf)值损失项的权重，w
(thd)
为(thd)值项的权重。
[0043]
上述方法中，结合基尔霍夫定律和卡尔曼滤波算法，建立了损失函数；在频率选择过程中，当处于非谐振频率时，谐波谐波畸变率(thd)k会偏大，应增大选频步进值，当靠近谐振频率时，(thd)k会减小，应缩小选频步进值；在工作过程中，在模态频率附近小范围δf内，根据pf值与(thd)k的变化趋势，总是朝着使l
k,f
值减小的趋势调节频率，使超声系统动态谐振，实现谐振频率的自适应跟踪。
[0044]
有益效果
[0045]
本发明公开的一种超声骨刀谐振频率的跟踪方法，其相比现有技术而言的有益效果在于：对谐波分量幅值的估计采用了卡尔曼滤波算法，同时兼顾功率因数和谐波畸变率作为约束条件，提高了滤波算法的鲁棒性，实现自适应跟踪，使超声骨刀的频率跟踪更为快速和准确，增大了输出效率和使用效果；同时，对谐波分量幅值的观测能够实时监测超声换能器和工作端部的健康状况，消除工作过程中因器件耐久度下降而引起事故的安全隐患，大大提高了超声骨刀的使用安全性。
附图说明
[0046]
图1是一种超声骨刀的谐振频率跟踪方法的主要电路拓扑框图；
[0047]
图2是本发明提供的超声骨刀谐振频率跟踪方法的流程示意图；
[0048]
图3是本发明提供的自适应频率跟踪方法效果图；
[0049]
图4是模型观测与自适应跟踪的频率偏差对比图；
[0050]
图5是考虑失真无功功率的功率关系图；
[0051]
图6是阻抗匹配后的超声换能器等效电路图。
具体实施方式
[0052]
下面结合附图和实施例对本发明进行更加详细的描述。
[0053]
本发明实例中，一种基于谐振分量的超声换能器谐振频率选择方法，该方法的具体电路包括微控制器、功率驱动电路、阻抗匹配网络、电信号采集模块、信号处理模块(如图1所示)；
[0054]
其中微控制器通过pwm接口与功率驱动电路电性连接；功率驱动电路与阻抗匹配网络为电性连接，将微处理器发出的小信号转换为所需的功率驱动信号；功率驱动电路与阻抗匹配网络电性连接，阻抗匹配网络与超声换能器的正负极电性连接，阻抗匹配网络具体为电感并联或变压器并联；电信号采集模块与阻抗匹配网络具体为经霍尔传感器连接或采样电阻连接；信号处理模块将采集到的电信号进行预处理，预处理包括放大，缩小，逻辑运算或去噪，最终信号传递给微控制器，微控制器将信号进行分析处理，其中包括实现算法运行，最后根据分析结果调整pwm输出，完成反馈。
[0055]
本发明提供的方法流程示意图参见图2，该方法具体包括以下步骤：
[0056]
步骤s1：获取超声换能器中的原始电信号，提取原始电信号中的各次谐波量测。霍尔传感器分别对电路电压信号和流过换能器的电流信号进行采样；超声换能器中电信号为
周期函数且满足狄利克雷(dirichlet)条件，因此可将电压分解为傅里叶级数形式使用离散傅里叶算法提取出电信号中各次谐波分量的幅值；从而可以求解出基波分量以及各次谐波分量。
[0057]
超声换能器的线电压可表示为:
[0058][0059]
其中u0为直流分量，k为谐波次数，uk为第k次电压谐波分量的幅值，m为谐波的最高次数，为基频，θ为相位
[0060]
对采集得到得的有限采样点数为n的电压序列{u(n)},其离散傅里叶变换为
[0061][0062]
设构造的复序列为{x(n)},即x(n)＝u(n)+ji(n),且0≤n≤n-1，复序列{x(n)}的fft为
[0063]
x(k)＝fft[u(n)]+jfft[i(n)]＝u(k)+ji(k
[0064]
x(k)＝fft[x(n)]＝xr(k)+jxi(k)0≤k≤n-1
[0065]
其中xr(k)为x(k)的实部，xi(k)为x(k)的虚部。
[0066]
u(n)}和{i(n)}对应的频谱可表示为:
[0067][0068][0069]
设u
k(t)
为被测信号u(t)的第k次谐波，则u
k(t)
可用向量表示为可用向量表示为同时
[0070]
最后可得出各次谐波分量的幅值为:
[0071][0072][0073]
步骤s2：根据k时刻第r次谐波幅值的量测计算k时刻第r次谐波幅值的最优估
计值对已知的超声换能器模态频率附近的机电模型，选取电压谐波为状态量，其过程误差服从高斯分布；选取经电压采集模块为观测器，所计算得到的谐波幅值为量测量，其测量误差满足统计学上时间序列的白噪声，建立基于卡尔曼滤波的谐波幅值最优估计模型。
[0074]
首先建立模态频率附近的线性系统模型，推导随时间递推的信号预测模型。
[0075]
当换能器等效电路中的l1和c1发生谐振时，其支路上只剩下动态电阻r1，此时等效为与静态电容c0的并联，匹配电感l将抵消电路中的容性部分，阻抗匹配后的超声换能器等效电路图如图6所示：
[0076]
换能器的总阻抗表达式为：
[0077]
z＝z1//z2//z3[0078]
其中，
[0079]
z1＝jωl
[0080][0081][0082]
根据欧姆定律超声换能器线电压中时间递推的预测信号离散形式可表示为：
[0083][0084]
式中uk为预测信号，ts为采样时间，fs为采样频率，m为谐波的最高次数，，ω表示频率，表示相位角；当r＝1时，表示基波分量的幅值，当r≠1时，为第r次谐波的幅值，vk是均值为零、协方差为是均值为零、协方差为的高斯白噪声。
[0085]
选取状态向量xk为
[0086]
xk＝[x
1,k x
2,k
...x
2r-1,k x
2r,k
...x
2m-1,k x
2m,k
]
t
[0087]
其中其中其中
[0088]
将时间递推预测信号转换成如下的状态空间方程形式
[0089][0090]
其中：
[0091][0092]
h＝[10
…
10
…
10]
[0093]
q＝diag([β
1 β2...β
2r-1
β
2r
...β
2m-1
β
2m
])
[0094]
q＝[ε]
[0095]
式中：ar表示为系统矩阵；w
k-1
表示过程噪声，且为协方差为q，均值为零的高斯白噪声；h表示观测矩阵；vk表示观测噪声，且为协方差为r，均值为零的高斯白噪声；式中β1β2...β
2r-1
β
2r
...β
2m-1
β
2m
和ε是待定值，需要进行
[0096]
卡尔曼滤波时在线估计出来。
[0097]
接着引入修正因子β，建立pf抗扰模型基础上的鲁棒卡尔曼滤波算法。常规的卡尔曼滤波算法中的过程噪声和量测噪声都假设满足高斯分布而且是独立的，通过在给定量测的前提下，递归得到状态的后验估计。但在超声骨刀的工作过程中，存在另一种噪声的情况，大部分时间为零，这种噪声称为稀疏噪声；例如传感器故障，观测异常，甚至一些故意干扰，这些都会造成稀疏噪声的产生，稀疏噪声不常见，但普遍存在，如果利用常规的卡尔曼滤波器去做最优估计，严重影响跟踪精度，甚至使滤波严重发散；所以严格地说超声骨刀属于包含稀疏噪声的超声系统，为了更好地进行谐振频率动态跟踪，需要设计一个鲁棒性强的滤波方法来解决这个问题。由于非线性负载，传输线路上会产生谐波失真，当谐波失真大到一定的程度，就不可忽略了，而超声换能器又是一种对谐波极为敏感的器件，因此真有功功率因数值(pf)在一定程度上反应了超声换能器的工作状态，能够作为抗扰模型的修正参考；提取模态频率附近的机电模型输出的真有功功率因数值(pf)，真有功功率因数值因其中p为有功功率，q为无功功率，d为失真无功功率(如图5所示)。
[0098]
建立pf预测模型：
[0099]ypf(i+1)
＝k(i+1)+b
[0100]
[k,b]
t
＝r(i,h)
[0101]
其中，y
pf(i+1)
为超声换能器的pf值与上一时刻时刻pf值的差值在i+1时刻的预测值；r(
·
,
·
)为最小二乘函数；向量i＝(i,i-1,
…
,i-3)
t
,i≥4表示连续量测四次的次数向量；向量h＝(y
pf(i)
,y
pf(i-1)
,
…
,y
pf(i-3)
)
t
,i≥4表示量测四次，当y
pf(i+1)
≥
±
0.1时，β取1；否则β取0；
[0102]
对y
pf(i)
，有：
[0103]ypf(i)
＝δpf
i-δpf
i-1
[0104]
式中，δpfi为i时刻超声换能器的pf值，δpf
i-1
为i-1时刻的pf值；
[0105]
引出pf抗扰模型：
[0106][0107]
式中为k-1时刻谐波幅值的鲁棒卡尔曼预测值；为k-1时刻的最优估计值，为k-1时刻谐波幅值的模型预测值，β为修正因子。
[0108]
建立卡尔曼滤波算法，需要先给定初值和初始协方差矩阵即
[0109][0110][0111]
式中，e表示期望值，var表示方差；然后建立先验谐波幅值预测值(先验估计)与后验谐波估计值(后验估计)的联系，为的进一步预测，即先验估计；为的最优估计，即后验估计。
[0112]
k时刻的一步预测与k-1时刻的最优估计之间的关系为：
[0113][0114]
此时系统预测结果已经更新，接着更新的协方差矩阵
[0115][0116]
式中，是对应的协方差矩阵，是对应的协方差矩阵，表示的转置矩阵，为系统噪声的协方差矩阵；而后验估计的最终输出表达式为：
[0117][0118]
上式实际上是以卡尔曼增益作为权值矩阵的“预测—校正”公式，既包含基于时间更新基础上的先验估计，由具有量测更新基础上的后验估计，所以采用卡尔曼滤波后最优估计值明显优于预测值是实际数据和预测之间的误差。
[0119]
卡尔曼增益矩阵的递推公式如下：
[0120][0121]
式中，是预测值的协方差矩阵，为量测噪声的协方差矩阵；
[0122]
为了使卡尔曼滤波算法递推下去，还需要更新的协方差矩阵
[0123]
[0124]
式中，i为单位矩阵。
[0125]
联合上述所得到的结果，形成环路反馈，可以得到基于卡尔曼滤波的谐波幅值最优估计算法如下式所示：
[0126][0127]
式中，k为时刻序列；r为谐波阶次；为给定初值；为给定初始协方差矩阵；为k时刻谐波幅值的一步预测值；为k-1时刻的系统传输矩阵；为k-1时刻的系统噪声；为的协方差矩阵；为的协方差矩阵；为的协方差矩阵；为k时刻的卡尔曼增益矩阵；为k时刻的量测矩阵；为量测噪声的协方差矩阵；为k时刻谐波幅值的鲁棒卡尔曼修正值；β为修正因子；为k时刻谐波幅值的最优估计值；i为单位矩阵；
[0128]
步骤s3：根据不同阶次的谐波幅值的最优估计值计算k时刻的谐波畸变率(thd)k。由于阶次越高，谐波分量越小，因而本发明实施例根据实际需求，取到第7次谐波；
[0129][0130]
对于电压信号，电压的谐波畸变率表示为(thdu)k：
[0131][0132]
对于电流信号，电压的谐波畸变率表示为(thdi)k：
[0133][0134]
步骤s4：根据损失函数值，判断超声换能器的工作状态，自适应跟踪频率。当超声换能器处于谐振状态工作时，其谐波畸变率最小，且功率因数(pf)相对较高；当超声换能器处于非谐振状态工作时，传输线路中谐波分量幅值增大，功率因数降低。pf值的理想值为1，(thd)k的理想值为0，设在k时刻频率为f时的pf值为(pf)
k,f
，总谐波畸变率为(thd)
k,f
，损失
值为l
k,f
，有：
[0135]
l
k,f
＝w
(pf)
[1-(pf)
k,f
]2+w
(thd)
[0-(thd)
k,f
]2[0136]
式中，w
(pf)
为(pf)值损失项的权重，w
(thd)
为(thd)值项的权重。
[0137]
在频率选择过程中，当处于非谐振频率时，谐波谐波畸变率(thd)k会偏大，应增大选频步进值，当靠近谐振频率时，(thd)k会减小，应缩小选频步进值；在工作过程中，在模态频率附近小范围δf内，根据pf值与(thd)k的变化趋势，总是朝着使l
k,f
值减小的趋势动态调节频率，使超声系统动态谐振，实现谐振频率的自适应跟踪。