基于相關分析法的脈沖響應辨識實驗設計

劉 恒, 孫 晉, 邰凡彬, 張易晨

(南京信息工程大學 電子與信息工程學院, 江蘇 南京 210044)

系統辨識是控制領域的一個十分重要分支。在多數情況下,被控對象的模型是不知道的,或者在正常運行期間模型的參數可能發生變化。在控制領域經常會碰到這樣的問題,并且模型未知的被控對象很難得到合理控制。為了能更好地分析和控制被控對象,首先需要建立它的數學模型,即系統辨識[1]。單位脈沖響應函數中包含了系統所有的動力特性參數,提取系統的單位脈沖函數或頻率響應函數是獲取系統動態特性參數、進行系統辨識以及對系統進行有效監測與精確控制的前提[2]。相關分析法是目前普遍受到重視的一種辯識脈沖響應的方法,當系統存在噪聲時,利用相關分析法辯識脈沖響應可以收到比較滿意的效果[3]。相關分析法脈沖響應辨識廣泛應用于大氣溫室效應建模[4]、地球物理電磁方法勘探的研究[5]、電機勵磁繞線變形[6]等。目前,相關分析法辨識脈沖響應主要基于Matlab軟件,通過軟件產生m序列信號,激勵待辨識系統,通過m序列信號和待辨識系統的輸出計算,得到相關分析法辨識的脈沖響應[7-9]。仿真有一定的空洞性,不利于學生對辨識方法的各個節點響應的理解。鮮有文獻報道硬件辨識實驗[10]。本實驗結合電子信息類學生的專業基礎,利用仿真軟件和實際硬件完成了實驗介紹,將相關分析法辨識系統脈沖響應仿真通過硬件系統來實現,提高學生對系統辨識中m序列信號、卷積運算等相關知識的理解。

1 相關分析法原理

相關分析法辨識系統脈沖響應是經典辨識方法中的一種。如果有2個時間函數,其中一個總是在任何時刻的值總是以某種確定關系依賴于另一個函數的值,則這2個函數是時間相關的。m 序列的特性類似于白噪聲,常在系統辨識過程中作為輸入信號,通過對輸入輸出信號的互相關運算,就可以得到待測系統的脈沖響應。

m序列是循環周期為Np·Δt、自相關函數近似于δ函數的一種隨機序列,其統計特性近似于白噪聲。m序列的自相關函數Rm(τ)為:

(1)

式中,a為m序列的幅值,τ為時間常數,Δt為移位寄存器移位周期,Np為m序列的長度。Np足夠大(選擇Np使m序列的循環周期大于系統的過渡過程時間)。

設數據的采樣時間等于m序列移位脈沖周期Δt,則Wiener-Hopf方程寫成離散形式[11]為

(2)

當m序列的循環周期Np·Δt大于系統的過渡過程時間tr時,脈沖響應在時間大于Np·Δt后基本上衰減為零。

(3)

實際應用中,式(3)中,j小于Np,式(3)同樣可表示為

(4)

式(4)中,M(j-τ)為離散后第j-τ個時刻的m序列信號,y(j)為第j時刻的系統激勵后的測量輸出。

系統輸入和輸出的時刻τ的相關函數可以通過對輸入輸出測量獲得。結合式(1)和(2),將式(3)展開為

(5)

假設:

對于一個穩定的待辨識系統來說,C是一個有界常數,而且很小(當Np充分大時),根據式(4)和(5),則有[11]:

(6)

待辨識的二階系統的傳遞函數為

時間常數T1=2.0 ms,T2=1.0 ms,靜態增益K=5.0。二階系統可以轉換為

(7)

式中，?02=1/(T1·T2)=500 000，

很顯然,系統阻尼ξ大于1,系統過阻尼,不存在超調量。對應的單位階躍響應h(t)為

(8)

式(8)中,M1和M2分別為:

隨著t的增大,h(t)的后兩項趨近于0,所以對應的穩態值為1,假定從t=0開始,達到穩定值的98%為過渡時間tr,根據式(8)可以得到:

(9)

求解式(9),得到過渡時間tr=9.2 ms。利用Multisim軟件搭建待辨識二階系統,二階系統由運算放大器和電阻、電容構成,搭建的電路原理見圖1。

圖1 二階待辨識系統電路原理圖

系統在?=0 rad/s時,s=j?=0,對應二階系統增益為5,此時幅頻曲線幅度為20log(5)=13.979 4 dB,見圖2。二階系統有低通濾波特性,系統的帶寬約為72 Hz,見圖3。

圖2 系統對應的幅度-頻率響應曲線(低頻)

圖3 系統帶寬仿真估算

用m序列作輸入信號辨識脈沖響應的步驟:

(1) 預估被辨識系統的過渡時間tr和系統的最高工作頻率fmax,tr和fmax是選擇m序列的重要依據,實驗中tr=9.2 ms,fmax=72Hz。

(2) 在系統辨識中,激勵輸出數據必須盡可能多的包含系統動態特性的信息,這和輸入信號的選擇有很大關系。當系統的頻率特性接近低通濾波特性時,m序列的參數Np和Δt應滿足下述條件:

(10)

(Np-1)·Δt≥tr

(11)

根據式(10)有Δt≤1/216=4.63 ms。對于m序列的階數N與參數Np之間存在如下關系：Np=2N-1。 利用m序列來激勵二階系統,通過系統的輸出及激勵信號的運算來達到辨識。根據相關分析法的理論推導,m序列階數越高,對應的辨識精度越高,但辨識算法的卷積計算的時間越長,同時采樣的周期也更長,使得辨識響應變慢。實驗綜合考慮后,結合式(10)和(11),選擇了6階m序列信號來激勵二階系統,Np=63,Δt=0.2 ms。通過對激勵信號和響應信號進行采樣,利用辨識算法得到系統的脈沖響應曲線。

2 辨識電路設計

辨識電路包括m序列產生電路和脈沖響應辨識電路。前者包括時鐘電路、移位寄存器及控制電路、序列變換電路；后者包括二階系統、電平抬升電路、主控電路、按鍵、數字顯示屏,見圖4。

圖4 辨識電路硬件電路模塊

555定時器產生時鐘信號,移位寄存器及控制電路由雙D觸發器及邏輯門構成,序列變換由運算放大器和電阻組成。待辨識二階系統由運算放大器和電容、電阻構成,電平抬升電路均由運算放大器和電阻構成,主控制器為STM32F103單片機最小系統,數字顯示屏為TFT屏SSD1289,8位獨立按鍵。

2.1 M序列產生電路

時鐘電路見圖5，由555定時器LM555及直插電阻、電容、二極管構成。輸出的方波信號高電平為+5 V,低電平為0 V。芯片采用5 V電源供電,調節電阻R1、R2、R3及電容C1就可以改變方波的頻率和占空比。其中輸出高電平時間tw1≈0.7(R1+R2+R3)C1,低電平時間tw2≈0.7(R2+R3)C1。設計的輸出方波頻率為5 kHz,對應周期為0.2 ms,仿真輸出波形見圖6,對應的硬件測試輸出波形見圖7,硬件測試與仿真一致,但電源提高到6 V,便于后續驅動觸發器。

圖5 時鐘電路

圖6 仿真時鐘波形

m序列信號由1和0構成,通過移位寄存器和控制電路來實現。邏輯“1”的個數比邏輯“0”多1。一種6階M序列信號,由6個移位寄存器組成。對應的m序列信號長度為63,序列信號的邏輯為{1,1,1,1,1,1,0,0,0,0,0,1,0,0,0,0,1,1,0,0,0,1,0,1,0,0,1,1,1,1,0,1,0,0,0,1,1,1,0,0,1,0,0,1,0,1,1,0,1,1,1,0,1,1,0,0,1,1,0,1,0,1,0}。對于6階m序列的反饋系數(八進制)為103,147,155,有3種實現方式,對應二進制數為1000011、1100111、1101101。7位分別對應D0、Q1、Q2、Q3、Q4、Q5、Q6,0表示沒有參與反饋,1表示參與了反饋。3種控制輸出D0= Q5⊕Q6,D0= Q1⊕Q4⊕Q5⊕Q6,D0= Q1⊕Q3 ⊕Q4⊕Q6。在本項目中,沒有采用異或門,采用與非門來實現異或門的邏輯功能。設計中,采用了2個2輸入1輸出與非門單元,1個8輸入1輸出與非門,1個3輸入1輸出與非門。沒有全部采用二輸入一輸出的與非門,便于更好地可靠連線,保障信號發生器的功能。圖8為由3片雙D觸發器和邏輯門電路構成的移位寄存器及控制電路，圖9為其仿真輸出波形,周期約為12.6 ms,剛好是周期0.2 ms的63倍,圖10為實測的輸出波形,與設計的期望一致。

圖7 實測時鐘輸出波形

圖8 移位寄存器及控制電路

圖9 仿真m序列信號波形

圖10 實測m序列信號

輸出的m序列信號高電平為+5 V,低電平為0 V,波形不具有幅度的對稱性。為此增加了m序列變換電路,見圖11,為線性電路。實驗設定幅值為+1 V和-1 V,仿真波形見圖12,實測波形見圖13,與仿真設計的一致。

圖11 m序列變換電路

圖12 仿真m序列變換后的波形

圖13 實測m序列變換后的波形

2.2 信號激勵及抬升電路

系統脈沖響應辨識需對激勵信號及響應信號進行卷積運算,需利用單片機等數字處理器來完成[12]。由于激勵信號及響應信號電平不在A/D采樣要求范圍內,需對2個信號進行抬升處理。激勵信號的幅值為-1 V～+1 V,響應信號的幅值為-0.5 V～+0.88 V,利用分壓電路和同相加法器將2個信號分別抬升1.45 V,滿足電平0～3.3 V范圍。圖14為對應的電平抬升電路,電路由電阻串聯分壓電路、電壓跟隨器電路組成,抬升電路參考電壓VR=5×4.1/14.1=1.454 V。抬升后m序列信號仿真電壓范圍+0.453 V～+2.453 V,硬件實測電壓范圍為+0.44 V～+2.40 V,電壓范圍略有變小,誤差為2%。

圖14 電壓抬升電路原理圖

變換后的m序列激勵二階系統的響應需能完全反映二階系統的特性。圖11輸出的SSS信號端連接圖1的SSS端。圖15是仿真對應圖1的Vos端輸出波形及電平抬升后的波形,系統激勵輸出波形具有周期性,周期為12.6 ms,激勵后二階系統的輸出電壓幅值范圍為-0.55 V～+0.77 V,利用圖14電路波形抬升了1.45 V后,電壓幅值為+0.9 V～+2.22 V。硬件實測激勵后的波形幅值為-0.5 V～+0.8 V,抬升后的實測波形,電壓幅值范圍為+0.9 V～+2.28 V,滿足設計預期。

圖15 仿真激勵輸出抬升前、后的波形

2.3 辨識算法及實現

電平抬升后,利用STM32F103主機自帶A/D模塊對抬升后的激勵信號和抬升后的系統輸出電壓進行兩通道采樣。在辨識算法嵌入式程序實現前,用Matlab軟件進行了仿真設計,程序由5個自主編寫和軟件自帶的其他函數組成,函數correlationAnalysis.m,函數transferFcn.m,函數estValue.m,函數ImpEstErr.m,函數theoreticalValue.m。

打開Matlab軟件,在主窗口輸入以下語句:

clear；% 清空工作空間的變量數據

clc； % 清除命令行的命令及多余的字符信息

然后打開correlationAnalysis.m函數,并在Matlab軟件選單里點擊“運行”控件,函數就開始運行,得到辨識均方誤差和辨識結果比較曲線。給定的二階系統,辨識結果與理論計算結果均方根最大誤差為0.019 3；對應的辨識結果和理論計算值在一個周期內比較見圖16；曲線非常接近,也驗證了軟件算法的有效性。

圖16 理論計算和辨識算法結果對比

圖17 系統軟件流程圖

圖18為實測及結果顯示,第一個波形為變換后的m序列,第二個波形為系統響應,黑色線為理論值,紅色為辨識算法得到的響應曲線,辨識結果與理論只有一定的偏差,這跟實物電路元件參數與仿真的不一致有關,上位機數據處理得到最大偏差在2%以內。

圖18 實測及結果顯示

3 結論

一種二階系統的脈沖響應辨識綜合教學實驗,包括m序列產生和脈沖響應測試2個子實驗。m序列采用555定時器、D觸發器、邏輯門、變換電路產生,涵蓋了模擬電子技術和數字電子技術課程知識,沒有采用可編程器件,不需復雜的編程和程序下載,增加了實驗的適應性。采用6階m序列信號,一方面考慮辨識的準確性,另一面兼顧了電路的復雜性。增加了m序列變換電路,提高了序列電平的對稱性,擴大了序列的應用范圍,同時將Multisim軟件應用到教學實驗設計中。采用了單片機自帶A/D模塊,不需增加硬件開銷,穩定性好。采用了基于相關分析法來對脈沖響應曲線進行辨識,相比其他算法,電路可實現性強,并通過Matlab軟件建模對算法進行驗證。可擴展藍牙模塊,在上位機上顯示輸入、輸出波形曲線及辨識曲線,方便直觀,學生可以掌握上位機編程。

致謝:感謝中國高等教育學會及浙江天煌科技實業有限公司對實驗設計的支持。

參考文獻(References)

[1] 鄧春龍,許偉明,張昕.相關函數脈沖響應法系統實時辨識[J].微計算機信息,2012,28(3):29-30.

[2] 于亮亮,宋漢文.環境激勵下脈沖響應函數與相關函數的關系[J].噪聲與振動控制,2017,33(3):14-18.

[3] 王燕平,陳昕志.利用相關分析法辨識系統脈沖響應[J].鄭州工業高等專科學校學報,2003,19(2):1-2.

[4] 孫興帥,秦琳琳,吳剛. 基于稀疏有限脈沖響應的溫室溫度系統建模[J].中國科學技術大學學報,2015,45(1):9-16.

1.4.1 心腦血管事件定義 心腦血管病事件定義為發生冠心病事件和（或）腦卒中事件。其中冠心病事件包括急性心肌梗死、冠心病猝死和其他冠心病死亡：腦卒中事件包括出血性腦卒中、缺血性腦卒中，但不包括一過性腦缺血和其他原因引起的腦血管病。在分析時，按冠心病事件和腦卒中事件進行統計分析，腦卒中事件再進一步劃分為缺血性腦卒中事件和出血性腦卒中事件。同一類事件如發生2次或2次以上，以首次發生的事件為終點事件。

[5] 武欣,薛國強,底青云,等.偽隨機編碼源電磁響應的精細辨識[J].地球物理學報,2015,58(8):2792-2802.

[6] 高佳平,沈煜,陳鵬,等.用偽隨機M 序列激勵的繞組變形試驗方法[J].中國電機工程學報,2016, 36(20):5678-5687.

[7] 方俊初,聶啟燕.多次采樣M序列法辨識LTI脈沖相應計算方法[J].河北工程大學學報,2016,33(3):109-112.

[8] 楊憲強.線性參數變化系統辨識方法研究[D]. 哈爾濱:哈爾濱工業大學,2014.

[9] 陳磊.連續系統的傳遞函數模型辨識[D].無錫:江南大學, 2012.

[10] 何華鋒,李元凱,董海迪,等.偽隨機信號(m序列)相關辨識的動態測試方法[J].中國測試,2013,39(4):88-92.

[11] 劉金琨,沈曉蓉,趙龍.系統辨識理論及Matlab仿真[M].北京:電子工業出版社,2013.

[12] 劉恒,吳朝陽,劉建成,等.一種典型閉環PID控制教學實驗設計[J].實驗技術與管理,2017,34(9):42-46.