基于CZT的變時間窗時變參數在線辨識方法及應用

鐘鴻豪 白文艷 黃萬偉

北京航天自動控制研究所，北京100854

由于建模在實際情況下不可避免地存在誤差，特別是某些參數隨著實際系統的工作條件、構型、故障等情況而發生變化，這些變化通常無法預知或測量，需要通過在線參數辨識來獲得[1]。

在線辨識受到算法計算時間和復雜程度的限制，目前只有少數幾種辨識方法可用于實時在線辨識。典型的有遞推最小二乘法、卡爾曼濾波法等，其中遞推最小二乘方法又分為時域和頻域2種形式[2]。時域的辨識方法，在測量噪聲較小的情況下，能夠得到較好的辨識結果，但是如果辨識方程中的微分項存在較大的測量噪聲，由于數值微分的影響，辨識效果會很差。頻域的遞推辨識方法的噪聲抑制能力強、計算量小、內存需求固定，在處理信號的微分時也非常簡便。但是，由于頻域方法辨識參數需要一定數據量的累積，以往的頻域辨識方法，都假定參數在辨識過程中固定不變，對于參數時變的情形適應性相對較差[3]。

本文提出的基于線性調頻Z變換(CZT)的變時間窗時變參數在線辨識方法，在保證有效減小測量噪聲影響的基礎上，利用滑動時間窗設計，提高了頻域方法對參數變化的敏感程度，使系統能夠較快地敏感參數變化，增強了頻域辨識方法對于參數時變情形的適應性，有效提高了傳統辨識方法對時變參數的在線辨識能力。最后，針對飛機系統執行機構的故障情況進行了仿真驗證，驗證了方法的有效性。

1 問題描述

考慮如下的系統模型：

(1)

其中x(t)=[x1(t),x2(t),…,xn(t)]T為狀態變量，可以通過計算得到，w(t)為過程噪聲。y(t)為量測量，v(t)為量測噪聲，θ(t)=[θ1(t),θ2(t),…,θn(t)]T表示待辨識的時變參數向量。

實際應用中，是在各個離散采樣時刻tk(k=1,2,…,N)，進行采樣計算。則對于每個采樣時刻tk，式(1)可以改寫為

Δxn(k)=θ1(k)x1(k)+θ2(k)x2(k)+
…+θn(k)xn(k)+w(k)
y(k)=xn(k)+v(k)

(2)

在時域辨識方法中，一般需要通過數值微分得到等式中的微分項，即

Δxn(k)≈(y(k)-y(k-1))/2Δt=
(xn(k)+v(k)-xn-1(k-1)-v(k-1))/2Δt=
(xn(k)-xn-1(k-1))/2Δt+
(v(k)-v(k-1))/Δt

(3)

其中Δt為采樣時間步長。因此可知，數值微分對于測量噪聲會有一個近似1/Δt倍的放大作用。

而在頻域中，對于在有限時間區間[t1,tk]內的信號x(k)，它的有限傅里葉變換可以定義為

(4)

可以近似為

而它的微分為

(5)

可以避免進行數值微分，進而避免數值微分對測量噪聲的放大作用。

另外，使用線性調頻Z變換，可以細化待分析頻段的頻譜，獲得很高的頻譜分辨率。因此，可以選擇一個離噪聲較遠的關心頻帶，進一步減弱噪聲的影響，從而獲得更好的辨識效果[4]。

然而傳統頻域方法，由于頻域方法本身需要數據量的累積，限制了頻域方法的辨識速度，導致辨識速度較慢。

本文的工作就是針對時變參數的在線辨識中，傳統時域和頻域辨識方法的問題，提出基于CZT的變時間窗時變參數在線辨識方法，實時準確地辨識待辨識參數。

2 基于CZT的變時間窗時變參數在線辨識方法

為了保證算法的去噪能力和計算的穩定，基于CZT的變時間窗時變參數在線辨識方法，需要預設一個較大的時間窗，作為缺省時間窗，比如4s的時間窗。這樣在參數變化不大的情況下，使用缺省的時間窗，可以利用較大的數據量，提高辨識算法的去噪能力和穩定性。

另一方面，為了提高頻域方法的辨識速度，在參數發生變化時，快速更新辨識結果，本方法設計了基于方差閾值的重啟算法。

首先，設定一個關心頻段[f0,f1)，并在頻段內選擇等間隔分布的M個頻率點，即fi=f0+iΔf,i=0,1,…,M-1。在基于時間窗的遞推CZT的基礎上，對于每個時刻tk=kΔt，首先計算出每個量的CZT：

(6)

(7)

(8)

(9)

(10)

否則認為該時刻參數發生了改變。下一個時刻重啟參數辨識。

此處，為了保證數據的穩定，設定一個較小的休眠時間tso，比如ts0=0.2s，使得參數值在[iΔt,iΔt+ts0]時間內，不發生變化，沿用上一刻的參數估計值，并且系統不再判斷方差變化，保證系統不會在休眠時間內再次重啟。

具體的基于CZT的變時間窗時變參數在線辨識流程如下：

step1:初始化參數。

step2: 判斷是否為初次啟動或重啟動時刻(startf的值是否為1)。

如果startf的值為1，進入step3，否則進入step4。

step3:啟動時刻的初始化。

首先記錄啟動時刻，令p=i。然后對每個頻率點(fi=f0+iΔf,i=0,1,…,M-1)的CZT變換賦初值。將startf置為0，計算休眠時間ts=Δt，累加序列長度length=1。

step4:判斷累加序列長度是否小于L。

如果累加序列長度小于L，進入step5，否則進入step6。

step5：利用CZT的遞推算法，累加計算，并根據休眠時間進行判斷。

累加求解每個頻率點(fi=f0+iΔf,i=0,1,…,M-1)的CZT變換值。計算休眠時間ts=ts+Δt，累加序列長度length=length+1。

然后判斷當前休眠時間和設定休眠總時間ts0的大小關系。

如果當前休眠時間大于設定休眠總時間ts0，則進入step7。

step6： 更新時間窗。

由于時間窗的長度為L，則CZT變換只采用最近L次采樣數據，在上一時刻的基礎上，舍棄記憶中最老的信息，加入當前最新的數據，進而更新每個頻率點的CZT變換值。然后，進入step7。

step7：計算方差，并進行判斷。

3 仿真試驗及分析

為了驗證所提出方法的有效性，針對F-16[6]飛機執行機構故障的條件下，對氣動參數進行了辨識。

下面簡要列出用于辨識的飛機縱向短周期運動的狀態空間模型，具體推導過程可參見文獻[7]：

(11)

其中，C1(t)、C2(t)與C3(t)為飛機氣動力系數，B1(t)、B2(t)與B3(t)為飛機氣動力矩系數；α(t)為飛機的攻角，ωz1(t)為飛機繞z1軸旋轉的角速度，δφ(t)為飛機的舵偏，這幾個狀態量和控制輸入均可測量或計算得到；e(t)為測量噪聲。

試驗選取子系統，進行辨識:

(12)

4s前氣動參數的真值為：

C1=-0.5709,C2=0.9638,C3=-0.0692
B2=-4.8765,B1=-0.9552,B3=-9.0689

4s時，執行機構突然出現故障，控制能力下降50%，即此時B3=-4.5344，其他條件不變。實驗中，時域帶遺忘因子的最小二乘法，遺忘因子選取為0.98，基于CZT的時間窗傳統頻域辨識方法，時間窗選取為L=400，即為2s。

3.1 無噪聲情況下辨識結果

在無噪聲情況，將基于CZT的變時間窗時變參數在線辨識方法與基于CZT的時間窗傳統頻域辨識方法對比，結果如圖1～3。

圖1 無噪聲情況下氣動參數B1(t)辨識結果

圖2 無噪聲情況下氣動參數B2(t)辨識結果

圖3 無噪聲情況下氣動參數B3(t)辨識結果

3.2 噪聲情況下辨識結果

為了驗證方法的去噪聲能力，在輸出信號中加入頻率為20Hz、信噪比為17dB的正弦噪聲，同樣可以得到2種方法的辨識結果如圖4～6。

圖4 噪聲情況下氣動參數B1(t)辨識結果

圖5 噪聲情況下氣動參數B2(t)辨識結果

圖6 噪聲情況下氣動參數B3(t)辨識結果

可以看出，在參數突變(系統發生故障)的情況下，仍能快速平滑地得到參數辨識結果(1s時誤差均在10%以內)，而且穩定時的誤差均在1%以內，驗證了方法的快速性和去噪性。

4 總結

提出了一種基于CZT的變時間窗時變參數在線辨識方法，不但能夠有效減小測量噪聲的影響，而且能夠較快地敏感參數變化，增強了頻域辨識方法對于參數時變情形的適應性，可以實現時變參數的在線辨識。并且針對飛機系統的執行機構故障情況下的氣動參數，進行了辨識，驗證了方法的適應性和準確性。