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

权利要求书

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中各状态参数的最佳估计值，所述最佳估计值中包含探头系统输入量的最佳估计值，所述探头系统输入量的最佳估计值即为通过反演补偿获得的探头系统输入量的最佳估计值。