CN113822354B - 基于贝叶斯反演算建模的微纳米探头动态特性补偿方法 - Google Patents

Abstract

本发明公开了一种基于贝叶斯反演算建模的微纳米探头动态特性补偿方法，其特征是通过微纳米探头系统的阶跃响应实测数据，采用自适应模型辨识法得到探头系统的传递函数模型；根据系统传递函数模型参数推导确定贝叶斯反演算补偿模型的观测回归矩阵和状态转移矩阵，获得动态补偿模型的状态方程和观测方程；采用参考先验分析法确定状态向量的初始先验分布，利用贝叶斯递推算法估计状态向量实时分布参数，得到补偿后的实时输入待估计量。本发明用于实现微纳米探头系统的动态特性补偿，能够快速有效地减小探头阶跃输出的超调量，缩短探头响应时间，以低成本改善微纳米探头系统的动态特性。

Description

基于贝叶斯反演算建模的微纳米探头动态特性补偿方法

技术领域

本发明涉及微纳米探头动态特性补偿技术领域，更具体地说是一种基于贝叶斯反演算建模的微纳米探头动态特性补偿方法。

背景技术

动态特性不理想是微纳米探头系统动态测量误差的重要来源，严重制约探头系统测量速度和精度的提升。微纳米探头系统结构复杂、样式多变，难以通过理论建模或结构改进实现动态特性提升。为了在提升系统测量速度同时保证测量精度，对其进行动态特性补偿是较为经济有效的办法，是在其测量输出端串联动态补偿器，要求动态补偿器模型与系统的动态特性模型互逆，从而实现动态特性补偿。

动态补偿器的设计依据主要有基于输出数据和基于系统模型两类。前者虽然无需辨识系统模型，但输出数据的误差会对动态补偿器的设计精度和稳健性产生较大的不良影响。因此，基于系统模型设计传感器的动态补偿器仍然是目前最常用的办法。为了获得动态补偿器，常用的补偿方法有：零极点配置法、最小二乘法和神经网络法等。其中，零极点配置法基于传感器的动态模型，能够通过输入输出数据准确地获知传感器的参数，但是对于动态特性复杂的传感器，很难获得高精度的参数模型。最小二乘法易于实现、应用广泛，但在辨识过程中存在不确定度干扰、数据饱和及参数偏差等问题。神经网络法收敛速度慢、容易陷入局部极值，而且需要训练网络，程序运行时间较长。现有技术中微纳米探头系统动态补偿方法普遍存在着实时性差和建模精度低的问题，不适用于对微纳米探头系统进行在线动态特性补偿。

发明内容

本发明是为避免上述现有技术所存在的不足，提供一种基于贝叶斯反演算建模的微纳米探头动态特性补偿方法，以实现快速、准确地建立探头系统的动态补偿模型，经济、有效地改善系统动态测试性能的发明目的。

为实现发明目的，本发明采用如下技术方案：

本发明基于贝叶斯反演算建模的微纳米探头动态特性补偿方法的特点是按如下步骤进行：

步骤1、由动态阶跃响应实验获得微纳米接触式探头系统的阶跃响应实测数据，利用所述阶跃响应实测数据，根据自适应模型参数辨识法建立探头系统的输入输出传递函数模型，获得由式(1)表征的n阶线性常系数微分方程：

式(1)中：

y_t为t时刻的探头系统输出，以

表征y_t的i阶导数，i＝1,2,…,n；

u_t为t时刻的探头系统输入，为待测量；{a_n-1,a_n-2,...,a₀,b}为模型参数序列；

步骤2、令x_1,t＝y_t,

将式(1)转化为由式(2)表征的状态空间方程：

式(2)中：

y_t和x_1,t均表示t时刻的探头系统输出，即：y_t＝x_1，t；

y_t+1和x_1,t+1均表示t+1时刻的探头系统输出，即：y_t+1＝x_1，t+1；

以x_j,t表示y_t的j-1阶导数，j＝2,…,n，因此有：

和x_2,t均为y_t的一阶导数；

和x_n,t均为y_t的n-1阶导数；

为y_t的n阶导数；

以

表示探头系统在动态响应中的采样时间间隔，设定

由前后相邻的k+1时刻和k时刻的状态值近似估计获得t时刻的一阶导数，将式(2)转化为由式(3)表征的方程组：

式(3)中：

u_k为k时刻的探头系统输入，为待估计量；

x_1,k表示k时刻的探头系统输出，y_k+1和x_1,k+1均表示k+1时刻的探头系统输出，即：y_k+1＝x_1,k+1；

以x_j,k表示k时刻状态参量x_j-1,k的单位时间增量，将x_j,k作为y_t的j-1阶导数的估计值；

以x_j,k+1表示k+1时刻状态参量x_j-1,k+1的单位时间增量，j＝2,…,n；

由式(3)分别获得式(4)和式(5)：

y_k+1＝x_1,k+1 (4)

式(4)表征测试序列的观测方程；式(5)表征状态序列方程组；

步骤3、采用线性逼近建模方法通过拟合获得未知的待估计量参数规律，是将u_k表示成一阶线性贝叶斯动态模型，获得由式(6)表征的观测方程和由式(7)表征的状态方程组：

观测方程：u_k＝μ_k+v_uk (6)

状态方程：

式(6)和式(7)中：

μ_k、α_k是u_k在k时刻的两个状态参数，分别是水平值及水平增量；

μ_k-1、α_k-1是u_k在k-1时刻的两个状态参数，分别是水平值及水平增量；

v_uk是k时刻u_k的观测噪声，ω_μk是k时刻μ_k的状态噪声，ω_αk是k时刻α_k的状态噪声；

v_uk、ω_μk和ω_αk相互独立且各自服从正态分布；

步骤4、令θ_k+1＝[x_1，k+1,…,x_n-1，k+1,x_n，k+1,μ_k,α_k]^T，θ_k+1是由各状态参数组成的k+1时刻的状态向量矩阵，联立式(4)、式(5)、式(6)和式(7)，并引入k+1时刻的观测噪声ν_y,k+1和状态向量噪声ω_θ,k+1获得由式(8)表征的贝叶斯反演算补偿模型的观测方程和由式(9)表征的状态方程：

观测方程：y_k+1＝[1 0 0 … 0]_1×(n+2)θ_k+1+ν_y,k+1 (8)

状态方程：

式(8)中：

以F^T表示观测回归矩阵，F^T＝[1 0 0 ... 0]_1×(n+2)；

ν_y,k+1服从均值为0、方差为V_k+1的正态分布，即：ν_y,k+1～N[0，V_k+1]；

式(9)中：以G表示状态转移矩阵，即：

ω_θ,k+1服从均值为0、方差为W_k+1的正态分布，即：ω_θ,k+1～N[0，W_k+1]；

对式(8)和式(9)进行简化获得由式(10)表征的贝叶斯反演算补偿模型：

步骤5、采用参考先验分析法和贝叶斯递推算法确定θ_k+1的后验分布，得到状态参数μ_k的最佳估计值，由此获得补偿后待估输入量u_k的最佳估计值，实现微纳米探头动态特性补偿。

本发明基于贝叶斯反演算建模的微纳米探头动态特性补偿方法的特点也在于：

以N表征数据量，将k时刻的取值范围0到N-1划分为第一阶段和第二阶段，第一阶段k取0,1,…,n+2，第二阶段k取n+3,n+4,…,N-1；将所述步骤5中为获得状态参数μ_k的最佳估计值按第一阶段和第二阶段分别进行：

第一阶段：k＝0,1,…,n+2，根据探头系统最初测得的n+3个阶跃响应测试序列，按如下方式采用参考先验分析法确定状态向量θ_k+1的先验分布及参数；

测量数据误差服从正态分布，状态噪声方差W_k+1＝0，并满足由式(11)表征的参考先验联合概率密度函数：

p(θ₁,V|D₀)∝V^-1 (11)

式(11)中：

以p(θ₁,V|D₀)表示状态参数θ₁和观测方差V的参考先验联合概率密度函数；

以V为观测方差，以V^-1表示观测方差V的倒数；

以D₀表示所有有效先验信息的集合；以∝表示正比于；

由式(12)表征观测数据的似然函数：

p(y_k+1|θ_k+1,V,D_k)∝V^-(1/2)exp{-0.5V^-1(y_k+1-F^Tθ_k+1)²} (12)

式(12)中：

以p(y_k+1|θ_k+1,V,D_k)表示y_k+1的似然函数；

以V^-(1/2)表示观测方差V的平方根的倒数；

以exp表示以无理数e为底的指数函数；

以D_k表示所有有效信息的集合，D_k＝{y_k,D_k-1}＝{y_k,···,y₁,D₀}；

结合式(12)并根据贝叶斯公式和探头系统模型，分别获得状态向量θ_k+1和观测方差V在k+1时刻的由式(13)表征的先验联合概率密度函数和由式(14)表征的后验联合概率密度函数：

式(13)和(14)中：

以p(θ_k+1,V|D_k)表示状态向量θ_k+1和观测方差V的先验联合概率密度函数，

以p(θ_k+1,V|D_k+1)表示状态向量θ_k+1和观测方差V的后验联合概率密度函数；

以H_k+1，h_k+1，λ_k+1，γ_k+1，L_k+1，l_k+1，L_k，l_k，δ_k+1和δ_k表示运算过程中产生的计算量，分别为：

H_k+1＝(G^-1)^TL_kG^-1，h_k+1＝(G^-1)^Tl_k，L_k+1＝H_k+1+FF^T，l_k+1＝h_k+1+Fy_k+1；

γ_k+1＝γ_k+1，λ_k+1＝δ_k，δ_k+1＝λ_k+1+y_k+1 ²；

各初始值为零，即：H₁＝0，h₁＝0，λ₁＝0，γ₀＝0；

由此获得由式(15)表征的θ_k+1的后验边缘分布和由式(16)表征的

的后验边缘分布：

式(15)中：

(θ_k+1|D_k+1)表示θ_k+1的后验边缘分布；～表示服从分布；

表示自由度为n_k的T分布，m_k+1为均值矩阵，C_k+1为方差矩阵，且：

S_k+1＝d_k+1/n_k+1；

式(16)中：

表示

的后验边缘分布；

Γ[n_k+1/2,d_k+1/2]表示自由度为n_k+1的Γ分布，d_k+1为尺度参数，且：

n_k+1＝γ_k+1-n-2，d_k+1＝δ_k+1-(l_k+1)^Tm_k+1；

第二阶段：k＝n+3,n+4,…,N-1，采用折扣因子法确定由式(17)表征的递推算法中k+1时刻W_k+1的值：

W_k+1＝GC_kG^T(ρ^-1-1) (17)

式(17)中：

C_k为θ_k+1后验边缘分布的方差矩阵；

ρ为折扣因子，折扣因子ρ是以观测数据预测方差最小为准则通过最优化搜索以确定；

根据探头系统输出y_k+1通过贝叶斯动态模型递推算法获得由式(18)表征的θ_k+1和V的后验联合概率密度函数：

p(θ_k+1,V|D_k+1)∝p(θ_k+1,V|D_k)p(y_k+1|θ_k+1,V,D_k) (18)

由式(18)获得由式(19)表征的θ_k+1的后验边缘分布和由式(20)表征的

的后验边缘分布：

θ_k+1的后验边缘分布均值矩阵m_k+1中相应的元素即为θ_k+1＝[x_1，k+1,…,x_n-1，k+1,x_n，k+1,μ_k,α_k]^T中各状态参数的最佳估计值，所述最佳估计值中包含探头系统输入量的最佳估计值，所述探头系统输入量的最佳估计值即为通过反演补偿获得的探头系统输入量的最佳估计值。

与已有技术相比，本发明有益效果体现在：

1、本发明在传递函数或微分方程已知的基础上直接建立贝叶斯反演算动态补偿模型，实现了快速、准确地探头系统动态特性补偿。

2、本发明基于贝叶斯统计原理，融合历史先验信息和当前样本信息进行参数递推估计，不仅所需样本少，而且还保障了动态补偿输出结果的可靠性与稳定性。

3、本发明的动态补偿方法可由软件编程实现，无需改变探头系统的原有结构，只要将该方法串联在系统的测量数据处理环节中，就能经济、有效地提高探头系统的动态测试性能。

附图说明

图1是补偿前的微纳米接触式探头系统阶跃响应实测数据图。

图2是基于贝叶斯反演算补偿探头系统阶跃响应实测数据的效果图。

具体实施方式

本实施例关于基于贝叶斯反演算建模的微纳米探头动态特性补偿方法，其微纳米探头是指微纳米接触式探头系统，其动态特性补偿是指针对微纳米接触式探头系统的实测输出数据进行动态补偿，以满足快速测量的需求；其动态特性补偿方法是按如下步骤通过动态阶跃响应实验、模型参数辨识以及动态补偿器构造，实现微纳米探头动态特性补偿：

步骤1、由动态阶跃响应实验获得微纳米接触式探头系统的阶跃响应实测数据，利用阶跃响应实测数据，根据自适应模型参数辨识法建立探头系统的输入输出传递函数模型，获得由式(1)表征的n阶线性常系数微分方程：

式(1)中：

y_t为t时刻的探头系统输出，以

表征y_t的i阶导数，i＝1,2,…,n；

u_t为t时刻的探头系统输入，为待测量；{a_n-1,a_n-2,…,a₀,b}为模型参数序列；

步骤2、为了得到贝叶斯反演算补偿模型，需要先把线性常系数微分方程转换成状态空间方程，然后再进行近似地离散化处理，获得测试序列的观测方程和状态序列方程组。因此，首先，假设x_1,t,x_2,t,…,x_n,t均连续可导，令

将式(1)转化为由式(2)表征的状态空间方程：

式(2)中：

y_t和x_1,t均表示t时刻的探头系统输出，即：y_t＝x_1，t；

y_t+1和x_1,t+1均表示t+1时刻的探头系统输出，即：y_t+1＝x_1，t+1；

以x_j,t表示y_t的j-1阶导数，j＝2,…,n，因此有：

和x_2,t均为y_t的一阶导数；

和x_n,t均为y_t的n-1阶导数；

为y_t的n阶导数；

考虑到探头系统在动态响应或扫描测量过程中，采样频率较高，采样时间间隔较短，以Δt表示探头系统在动态响应中的采样时间间隔，设定Δt→0，由前后相邻的k+1时刻和k时刻的状态值近似估计获得t时刻的一阶导数，将式(2)转化为由式(3)表征的方程组：

式(3)中：

u_k为k时刻的探头系统输入，为待估计量；

x_1,k表示k时刻的探头系统输出，y_k+1和x_1,k+1均表示k+1时刻的探头系统输出，即：y_k+1＝x_1,k+1；

以x_j,k表示k时刻状态参量x_j-1,k的单位时间增量，将x_j,k作为y_t的j-1阶导数的估计值；

以x_j,k+1表示k+1时刻状态参量x_j-1,k+1的单位时间增量，j＝2,…,n；

由式(3)分别获得式(4)和式(5)：

y_k+1＝x_1,k+1 (4)

式(4)表征测试序列的观测方程；式(5)表征状态序列方程组；

步骤3、贝叶斯反演算补偿的主要目的是获得输入量u_t的最佳估计值，因此，实际测量时，公式(5)中的输入量u_k也是待估计的未知量。为了反映待估计量的一般特性，采用线性逼近建模方法通过拟合获得未知的待估计量参数规律，是将u_k表示成一阶线性贝叶斯动态模型，获得由式(6)表征的观测方程和由式(7)表征的状态方程组：

观测方程：u_k＝μ_k+v_uk (6)

状态方程：

式(6)和式(7)中：

μ_k、α_k是u_k在k时刻的两个状态参数，分别是水平值及水平增量；

μ_k-1、α_k-1是u_k在k-1时刻的两个状态参数，分别是水平值及水平增量；

v_uk是k时刻u_k的观测噪声，ω_μk是k时刻μ_k的状态噪声，ω_αk是k时刻α_k的状态噪声；

v_uk、ω_μk和ω_αk相互独立且各自服从正态分布。如果对被测量有一些其他的已知先验信息，如被测量是一个平面，则可以采用贝叶斯常均值模型来拟合；如已知被测量是一个曲面，则可以考虑与曲面特性相一致的多项式模型来表征u_k。

步骤4、根据式(4)表征的测试序列的观测方程、式(5)表征的状态序列方程组以及输入量u_k的一阶线性贝叶斯动态模型，建立贝叶斯反演算补偿模型的观测方程和状态方程。

令θ_k+1＝[x_1，k+1,…,x_n-1，k+1,x_n，k+1,μ_k,α_k]^T，θ_k+1是由各状态参数组成的k+1时刻的状态向量矩阵，考虑到测量过程中不可避免会受到随机干扰的影响，联立式(4)、式(5)、式(6)和式(7)，并引入k+1时刻的观测噪声ν_y,k+1和状态向量噪声ω_θ,k+1获得由式(8)表征的贝叶斯反演算补偿模型的观测方程和由式(9)表征的状态方程：

观测方程：y_k+1＝[1 0 0 … 0]_1×(n+2)θ_k+1+ν_y,k+1(8)

状态方程：

式(8)中：

以F^T表示观测回归矩阵，F^T＝[1 0 0 … 0]_1×(n+2)；

ν_y,k+1服从均值为0、方差为V_k+1的正态分布，即：ν_y,k+1～N[0，V_k+1]；

式(9)中：

以G表示状态转移矩阵，即：

ω_θ,k+1服从均值为0、方差为W_k+1的正态分布，即：ω_θ,k+1～N[0，W_k+1]；

对式(8)和式(9)进行简化获得由式(10)表征的贝叶斯反演算补偿模型：

步骤5、采用参考先验分析法和贝叶斯递推算法确定θ_k+1的后验分布，得到状态参数μ_k的最佳估计值，由此获得补偿后待估输入量u_k的最佳估计值，实现微纳米探头动态特性补偿。

具体实施中，以N表征数据量，将k时刻的取值范围0到N-1划分为第一阶段和第二阶段，第一阶段k取0,1,…,n+2，第二阶段k取n+3,n+4,…,N-1；将步骤5中为获得状态参数μ_k的最佳估计值按第一阶段和第二阶段分别进行：

第一阶段：k＝0,1,…,n+2，根据探头系统最初测得的n+3个阶跃响应测试序列，按如下方式采用参考先验分析法确定状态向量θ_k+1的先验分布及参数；

在参考先验分析中，探头系统共有n+3个待估参数，即θ_k+1＝[x_1，k+1,…,x_n-1，k+1,x_n，k+1,μ_k,α_k]^T和V，需要采用探头前n+3个动态测量数据得到合适的初始分布。由于每一个参数实际上只有一个观测值能够提供相关信息，所以不可能估计或探测出参数的任何变化，因此合理假定状态噪声方差W_k+1＝0。

测量数据误差服从正态分布，状态噪声方差W_k+1＝0，并满足由式(11)表征的参考先验联合概率密度函数：

p(θ₁,VD₀)∝V^-1 (11)

式(11)中：

以p(θ₁,VD₀)表示状态参数θ₁和观测方差V的参考先验联合概率密度函数；

以V为观测方差，以V^-1表示观测方差V的倒数；

以D₀表示所有有效先验信息的集合；以∝表示正比于；

由式(12)表征观测数据的似然函数：

p(y_k+1|θ_k+1,V,D_k)∝V-^(1/2)exp{-0.5V^-1(y_k+1-F^Tθ_k+1)²} (12)

式(12)中：

以p(y_k+1|θ_k+1,V,D_k)表示y_k+1的似然函数；

以V^-(1/2)表示观测方差V的平方根的倒数；

以exp表示以无理数e为底的指数函数；

以D_k表示所有有效信息的集合，D_k＝{y_k,D_k-1}＝{y_k,···,y₁,D₀}；

结合式(12)并根据贝叶斯公式和探头系统模型，分别获得状态向量θ_k+1和观测方差V在k+1时刻的由式(13)表征的先验联合概率密度函数和由式(14)表征的后验联合概率密度函数：

式(13)和(14)中：

以p(θ_k+1,V|D_k)表示状态向量θ_k+1和观测方差V的先验联合概率密度函数，

以p(θ_k+1,V|D_k+1)表示状态向量θ_k+1和观测方差V的后验联合概率密度函数；

以H_k+1，h_k+1，λ_k+1，γ_k+1，L_k+1，l_k+1，L_k，l_k，δ_k+1和δ_k表示运算过程中产生的计算量，分别为：

H_k+1＝(G^-1)^TL_kG^-1，h_k+1＝(G^-1)^Tl_k，L_k+1＝H_k+1+FF^T，l_k+1＝h_k+1+Fy_k+1；

γ_k+1＝γ_k+1，λ_k+1＝δ_k，δ_k+1＝λ_k+1+y_k+1 ²；

各初始值为零，即：H₁＝0，h₁＝0，λ₁＝0，γ₀＝0；

由此获得由式(15)表征的θ_k+1的后验边缘分布和由式(16)表征的

的后验边缘分布：

式(15)中：

(θ_k+1|D_k+1)表示θ_k+1的后验边缘分布；～表示服从分布；

T_nk[m_k+1,C_k+1]表示自由度为n_k的T分布，m_k+1为均值矩阵，C_k+1为方差矩阵，且：

S_k+1＝d_k+1/n_k+1；

式(16)中：

表示

的后验边缘分布；

Γ[n_k+1/2,d_k+1/2]表示自由度为n_k+1的Γ分布，d_k+1为尺度参数，且：

n_k+1＝γ_k+1-n-2，d_k+1＝δ_k+1-(l_k+1)^Tm_k+1；

第二阶段：k＝n+3,n+4,…,N-1，采用折扣因子法确定由式(17)表征的递推算法中k+1时刻W_k+1的值：

W_k+1＝GC_kG^T(ρ^-1-1) (17)

式(17)中：

C_k为θ_k+1后验边缘分布的方差矩阵；

ρ为折扣因子，折扣因子ρ是以观测数据预测方差最小为准则通过最优化搜索以确定；

按照参考先验分析法，根据探头系统输出y_k+1通过贝叶斯动态模型递推算法获得由式(18)表征的状态向量θ_k+1和观测方差V的后验联合概率密度函数：

p(θ_k+1,V|D_k+1)∝p(θ_k+1,V|D_k)p(y_k+1|θ_k+1,V,D_k) (18)

由式(18)获得由式(19)表征的θ_k+1的后验边缘分布和由式(20)表征的

的后验边缘分布：

θ_k+1的后验边缘分布均值矩阵m_k+1中相应的元素即为θ_k+1＝[x_1，k+1,…,x_n-1，k+1,x_n，k+1,μ_k,α_k]^T中各状态参数的最佳估计值，最佳估计值中包含探头系统输入量的最佳估计值，探头系统输入量的最佳估计值即为通过反演补偿获得的探头系统输入量的最佳估计值。

利用本发明方法按如下过程针对公开号为CN104457613A的发明专利说明书中提出的一种三维微纳米接触触发式探头进行动态特性补偿，其结构简单、装调方便，分辨力达到1nm，重复性为8.2nm，允许触碰范围为±10μm。

首先对探头系统进行阶跃响应实验测试，采样间隔为2.5×10^-4s，采样时长为4s，对探头系统输出数据进行归一化处理得到用于建模的测试数据，探头系统的阶跃响应实测数据如图1所示。

根据阶跃响应实测数据，采用自适应模型辨识法确定探头系统为典型的二阶系统，即n取值为2，探头系统输入输出传递函数模型G(s)由式(21)所表征：

式(21)中，s为复频率。

由式(21)计算得到探头系统的阻尼比ζ为0.0099、固有频率ω_n为158Hz、以及由式(22)表征的二阶线性常系数微分方程：

令x_1,t＝y_t,

将式(22)转化成由式(23)表征的方程组：

将式(6)和式(7)代入到式(23)中，可得由式(24)表征的贝叶斯反演算模型：

第一阶段：设定k＝0,1,2,3,4，假定测量数据误差服从正态分布，状态噪声方差W_k+1＝0，并满足由式(11)表征的参考先验概率密度函数，由式(12)表征的观测数据的似然函数：

在参考先验分析中，二阶探头系统共有5个待估参数，即θ_k+1＝[x_1，k+1,x_2，k+1,μ_k,α_k]^T和V，需要采用探头前5个动态测量数据得到合适的初始分布。由于每一个参数实际上只有一个观测值能够提供相关信息，所以不可能估计或探测出参数的任何变化，因此合理假定状态噪声方差W_k+1＝0。

按照参考先验分析法，根据贝叶斯公式、探头系统模型以及探头前5个动态测量数据，分别获得状态向量θ_k+1和观测方差V在k+1时刻的由式(13)表征的先验联合概率密度函数和由式(14)表征的后验联合概率密度函数，再分别获得由式(15)表征的θ_k+1的后验边缘分布和由式(16)表征的

的后验边缘分布，并作为贝叶斯动态模型递推算法所需的先验信息。

第二阶段：设定k＝5,6,…,15999，采用折扣因子法确定由式(17)表征的递推算法中k+1时刻W_k+1的值。

根据探头系统输出y_k+1通过贝叶斯动态模型递推算法获得由式(18)表征的θ_k+1的后验联合概率密度函数；由式(18)获得由式(19)表征的θ_k+1的后验边缘分布和由式(20)表征的

的后验边缘分布。θ_k+1的后验边缘分布均值矩阵m_k+1中相应的元素即为θ_k+1＝[x_1，k+1,x_2，k+1,μ_k,α_k]^T中各状态参数的最佳估计值，最佳估计值中包含探头系统输入量的最佳估计值，探头系统输入量的最佳估计值即为通过反演补偿获得的探头系统输入量的最佳估计值。

微纳米探头系统动态特性补偿效果展示：

微纳米接触式探头系统的阶跃响应实测数据补偿后效果如图2所示，补偿前后探头系统的动特性参数如表1所示。可见：补偿后探头系统阶跃响应的超调量从99.15％降低到40.05％，响应时间由2.4825s缩短到0.1720s。从应用的角度出发，补偿后探头阶跃输出的响应时间很短，在实际测量过程中，一般会在输出信号稳定之后开始采集有效数据，所以在响应时间之前的输出信号的超调量不会对最终结果产生影响。因此，该贝叶斯反演算补偿结果满足动态测量的需求，探头系统的动态特性得到了较大改善。

表1贝叶斯反演算补偿前后效果比较

动特性参数	动态补偿前	贝叶斯反演算
			上升时间t<sub>r</sub>/s	0.0016	0.0043
峰值时间t<sub>p</sub>/s	0.0033	0.0067
			响应时间t<sub>s</sub>/s	2.4825	0.1720
超调量σ<sub>p</sub>/％	99.15％	40.05％

尽管已经示出和描述了本发明的实施例，对于本领域的技术人员，可以在不脱离本发明的原理和精神的情况下对实施例进行多种变化、修改、替换和变型，本发明的范围由所附权利要求及其等同物限定。

Claims (1)

1.一种基于贝叶斯反演算建模的微纳米探头动态特性补偿方法，其特征是按如下步骤进行：

步骤1、由动态阶跃响应实验获得微纳米接触式探头系统的阶跃响应实测数据，利用所述阶跃响应实测数据，根据自适应模型参数辨识法建立探头系统的输入输出传递函数模型，获得由式(1)表征的n阶线性常系数微分方程：

式(1)中：

y_t为t时刻的探头系统输出，以

表征y_t的i阶导数，i＝1,2,…,n；

u_t为t时刻的探头系统输入，为待测量；{a_n-1,a_n-2,...,a₀,b}为模型参数序列；

步骤2、令

将式(1)转化为由式(2)表征的状态空间方程：

式(2)中：

y_t和x_1,t均表示t时刻的探头系统输出，即：y_t＝x_1，t；

y_t+1和x_1,t+1均表示t+1时刻的探头系统输出，即：y_t+1＝x_1，t+1；

以x_j,t表示y_t的j-1阶导数，j＝2,…,n，因此有：

和x_2,t均为y_t的一阶导数；

和x_n,t均为y_t的n-1阶导数；

为y_t的n阶导数；

以Δt表示探头系统在动态响应中的采样时间间隔，设定Δt→0，由前后相邻的k+1时刻和k时刻的状态值近似估计获得t时刻的一阶导数，将式(2)转化为由式(3)表征的方程组：

式(3)中：

u_k为k时刻的探头系统输入，为待估计量；

x_1,k表示k时刻的探头系统输出，y_k+1和x_1,k+1均表示k+1时刻的探头系统输出，即：y_k+1＝x_1,k+1；

以x_j,k表示k时刻状态参量x_j-1,k的单位时间增量，将x_j,k作为y_t的j-1阶导数的估计值；

以x_j,k+1表示k+1时刻状态参量x_j-1,k+1的单位时间增量，j＝2,…,n；

由式(3)分别获得式(4)和式(5)：

y_k+1＝x_1,k+1 (4)

式(4)表征测试序列的观测方程；式(5)表征状态序列方程组；

步骤3、采用线性逼近建模方法通过拟合获得未知的待估计量参数规律，是将u_k表示成一阶线性贝叶斯动态模型，获得由式(6)表征的观测方程和由式(7)表征的状态方程组：

观测方程：u_k＝μ_k+v_uk (6)

状态方程：

式(6)和式(7)中：

μ_k、α_k是u_k在k时刻的两个状态参数，分别是水平值及水平增量；

μ_k-1、α_k-1是u_k在k-1时刻的两个状态参数，分别是水平值及水平增量；

v_uk是k时刻u_k的观测噪声，ω_μk是k时刻μ_k的状态噪声，ω_αk是k时刻α_k的状态噪声；

v_uk、ω_μk和ω_αk相互独立且各自服从正态分布；

步骤4、令

θ_k+1是由各状态参数组成的k+1时刻的状态向量矩阵，联立式(4)、式(5)、式(6)和式(7)，并引入k+1时刻的观测噪声ν_y,k+1和状态向量噪声ω_θ,k+1获得由式(8)表征的贝叶斯反演算补偿模型的观测方程和由式(9)表征的状态方程：观测方程：y_k+1＝[1 0 0 … 0]_1×(n+2)θ_k+1+ν_y,k+1 (8)

状态方程：

式(8)中：

以F^T表示观测回归矩阵，F^T＝[1 0 0 … 0]_1×(n+2)；

ν_y,k+1服从均值为0、方差为V_k+1的正态分布，即：ν_y,k+1～N[0，V_k+1]；

式(9)中：以G表示状态转移矩阵，即：

ω_θ,k+1服从均值为0、方差为W_k+1的正态分布，即：ω_θ,k+1～N[0，W_k+1]；

对式(8)和式(9)进行简化获得由式(10)表征的贝叶斯反演算补偿模型：

步骤5、采用参考先验分析法和贝叶斯递推算法确定θ_k+1的后验分布，得到状态参数μ_k的最佳估计值，由此获得补偿后待估输入量u_k的最佳估计值，实现微纳米探头动态特性补偿：

以N表征数据量，将k时刻的取值范围0到N-1划分为第一阶段和第二阶段，第一阶段k取0,1,…,n+2，第二阶段k取n+3,n+4,…,N-1；将所述步骤5中为获得状态参数μ_k的最佳估计值按第一阶段和第二阶段分别进行：

第一阶段：k＝0,1,…,n+2，根据探头系统最初测得的n+3个阶跃响应测试序列，按如下方式采用参考先验分析法确定状态向量θ_k+1的先验分布及参数；

测量数据误差服从正态分布，状态噪声方差W_k+1＝0，并满足由式(11)表征的参考先验联合概率密度函数：

p(θ₁,V|D₀)∝V^-1 (11)

式(11)中：

以p(θ₁,V|D₀)表示状态参数θ₁和观测方差V的参考先验联合概率密度函数；

以V为观测方差，以V^-1表示观测方差V的倒数；

以D₀表示所有有效先验信息的集合；以∝表示正比于；

由式(12)表征观测数据的似然函数：

p(y_k+1|θ_k+1,V,D_k)∝V^-(1/2)exp{-0.5V^-1(y_k+1-F^Tθ_k+1)²} (12)

式(12)中：

以p(y_k+1|θ_k+1,V,D_k)表示y_k+1的似然函数；

以V^-(1/2)表示观测方差V的平方根的倒数；

以exp表示以无理数e为底的指数函数；

以D_k表示所有有效信息的集合，D_k＝{y_k,D_k-1}＝{y_k,…,y₁,D₀}；

结合式(12)并根据贝叶斯公式和探头系统模型，分别获得状态向量θ_k+1和观测方差V在k+1时刻的由式(13)表征的先验联合概率密度函数和由式(14)表征的后验联合概率密度函数：

式(13)和(14)中：

以p(θ_k+1,V|D_k)表示状态向量θ_k+1和观测方差V的先验联合概率密度函数，

以p(θ_k+1,V|D_k+1)表示状态向量θ_k+1和观测方差V的后验联合概率密度函数；

以H_k+1，h_k+1，λ_k+1，γ_k+1，L_k+1，l_k+1，L_k，l_k，δ_k+1和δ_k表示运算过程中产生的计算量，分别为：

H_k+1＝(G^-1)^TL_kG^-1，h_k+1＝(G^-1)^Tl_k，L_k+1＝H_k+1+FF^T，l_k+1＝h_k+1+Fy_k+1；

γ_k+1＝γ_k+1，λ_k+1＝δ_k，δ_k+1＝λ_k+1+y_k+1 ²；

各初始值为零，即：H₁＝0，h₁＝0，λ₁＝0，γ₀＝0；

由此获得由式(15)表征的θ_k+1的后验边缘分布和由式(16)表征的

的后验边缘分布：

式(15)中：

(θ_k+1|D_k+1)表示θ_k+1的后验边缘分布；～表示服从分布；

T_nk[m_k+1,C_k+1]表示自由度为n_k的T分布，m_k+1为均值矩阵，C_k+1为方差矩阵，且：

S_k+1＝d_k+1/n_k+1；

式(16)中：

表示

的后验边缘分布；

Γ[n_k+1/2,d_k+1/2]表示自由度为n_k+1的Γ分布，d_k+1为尺度参数，且：

n_k+1＝γ_k+1-n-2，d_k+1＝δ_k+1-(l_k+1)^Tm_k+1；

第二阶段：k＝n+3,n+4,…,N-1，采用折扣因子法确定由式(17)表征的递推算法中k+1时刻W_k+1的值：

W_k+1＝GC_kG^T(ρ^-1-1) (17)

式(17)中：

C_k为θ_k+1后验边缘分布的方差矩阵；

ρ为折扣因子，折扣因子ρ是以观测数据预测方差最小为准则通过最优化搜索以确定；

根据探头系统输出y_k+1通过贝叶斯动态模型递推算法获得由式(18)表征的θ_k+1和V的后验联合概率密度函数：

p(θ_k+1,V|D_k+1)∝p(θ_k+1,V|D_k)p(y_k+1|θ_k+1,V,D_k) (18)

由式(18)获得由式(19)表征的θ_k+1的后验边缘分布和由式(20)表征的

的后验边缘分布：

θ_k+1的后验边缘分布均值矩阵m_k+1中相应的元素即为θ_k+1＝[x_1，k+1,…,x_n-1，k+1,x_n，k+1,μ_k,α_k]^T中各状态参数的最佳估计值，所述最佳估计值中包含探头系统输入量的最佳估计值，所述探头系统输入量的最佳估计值即为通过反演补偿获得的探头系统输入量的最佳估计值。

Images