魏 娟,張志杰,趙晨陽,李巖峰

(中北大學儀器科學與動態(tài)測試教育部重點實驗室,太原 030051)

一定當量的炸藥引爆后,能在極短的時間內(nèi)釋放出大量的能量,產(chǎn)生高溫,同時釋放大量氣體,向周圍介質施加高壓,接著高壓氣浪繼續(xù)向外膨脹并壓縮周圍空氣介質,使其壓力、密度、溫度等物理量發(fā)生突躍式變化,從而形成沖擊波[1-2]。沖擊波的超壓峰值、正壓作用時間和比沖量是估計炸藥毀傷效果的主要參數(shù)[3]。

試驗中通常利用沖擊波存儲測試系統(tǒng)對沖擊波信號直接進行采集,但系統(tǒng)中的壓力傳感器動態(tài)響應特性不能滿足測試系統(tǒng)的需求,有效帶寬不能覆蓋信號的頻譜。在沖擊波信號的陡峭前沿的激勵下,在諧振頻率點附近信號的幅值大幅度增加,測試數(shù)據(jù)出現(xiàn)激烈震蕩,進而產(chǎn)生較大的動態(tài)誤差,影響測試精度[4-6]。因此需要對壓力傳感器及測試系統(tǒng)的動態(tài)特性進行補償,提高傳感器的響應速度、降低超調量、擴寬工作頻帶以達到動態(tài)測試需求[7]。

1 動態(tài)補償濾波器設計

1.1 動態(tài)補償原理

在實際測量中,被測信號總是不斷變化的,通過測試系統(tǒng)傳感器的輸出結果來獲取、估計被測信號,這就要求傳感器能夠實時地、無失真地跟蹤被測信號的變化過程,因此需要研究傳感器的動態(tài)響應特性[8]。傳感器的動態(tài)特性是指傳感器對隨時間變化的輸入量的響應特性。動態(tài)特性好的傳感器,其輸入量隨時間變化的曲線與被測量隨同一時間變化的曲線一致或近似[9]。實際動態(tài)校準系統(tǒng)中,根據(jù)“標準”輸入信號,觀測傳感器的響應信號來估計傳感器的響應特性。

傳感器動態(tài)特性補償?shù)哪康氖鞘箓鞲衅鞯挠行Чぷ黝l帶能夠包含被測信號的所有頻譜分量,進而滿足動態(tài)測試需求。本文采用逆向濾波的方式(如圖1),在傳感器系統(tǒng)后接補償濾波器,達到擴寬頻帶的目的。根據(jù)QR分解估計逆向濾波器模型的階數(shù),改進粒子群(PSO)算法估計濾波器系數(shù),從而得到濾波器傳遞函數(shù)。將所得到的濾波器進行硬件實現(xiàn),用于實際沖擊波測試中,達到實時動態(tài)補償目的。

圖1 動態(tài)補償原理框圖

1.2 算法原理

動態(tài)補償濾波器的模型估計首先需辨識模型階數(shù)。壓力傳感器系統(tǒng)一般等效為二階系統(tǒng),但在試驗或實驗中,由于輸入信號的激勵作用,它的輸出信號在諧振點處容易發(fā)生震蕩[10],引入了動態(tài)誤差。因此以二階系統(tǒng)作為補償模型的最終估計偏差較大。

本文以激波管產(chǎn)生的階躍信號作為“標準”輸入信號,對被校準傳感器進行激勵,由采集電路獲取響應數(shù)據(jù)并進行處理,得到建模所用樣本數(shù)據(jù)。對樣本數(shù)據(jù)進行QR分解,求取殘差平方和,判別最優(yōu)階數(shù)。構建實驗數(shù)據(jù)矩陣D={(u(k),y(k))|k=1,2,…,N}

式中:v為待檢驗的階數(shù),N為數(shù)據(jù)長度,對矩陣進行QR分解,可以得到

R′是一個m×m維的上三角方陣。計算R′陣對角線偶元素平方和數(shù)值即為各階差分方程模型對應于最小二乘法估計的殘差平方和[11]。通過比較各階模型殘差平方和的大小,確定補償模型適合的階數(shù)。

濾波器的系數(shù)由改進粒子群(PSO)算法得到,粒子群算法是一種基于群體智能的全局優(yōu)化計算技術,具有高精度的穩(wěn)定性、并行性和全局搜索能力。在此算法中,每個粒子通過不斷更正空間方位在D維空間中尋找最佳方案。每個粒子搜尋過程如圖2所示。

圖2 粒子搜尋過程示意圖

每個粒子都有獨立的位置Xi=(xi1,xi2,…,xiD)和速度Vi=(vi1,vi2,…,viD),并且這兩個因素隨時間和空間不斷變化,各個粒子按照式(1)(2)更新其速度和位置[12]。

vid(n+1)=wvid(n)+c1r1d(n)(pid-xid(n))
+c2r2d(n)(pgd-xid(n))

(1)

xid(n+1)=xid(n)+vid(n+1)

(2)

式中:w為慣性權,c1,c2為加速系數(shù),r1d,r2d為在[0,1]內(nèi)均勻分布的隨機數(shù),n為當前迭代次數(shù)。選擇合適的慣性權值w使全局和局部搜索之間達到動態(tài)平衡后,能夠具有較快的收斂速度[13],故本文中選取線性遞減權值(LDIW)策略,慣性權值按照式(3)進行更新:

(3)

式中:M為最大允許迭代次數(shù),i為當前的迭代次數(shù)。

計算各個粒子的適應值與當前最優(yōu)位置的適應值進行比較,根據(jù)適應度函數(shù)的運算結果,在D維空間中不斷更新自己的最優(yōu)位置Pbest=(pi1,pi2,…,piD)和全局最優(yōu)位置Pg=(pg1,pg2,…,pgD)。本文中適應度函數(shù)(式4)采用均方誤差值進行比較,具體的優(yōu)化流程如圖3所示。

(4)

式中:yi為預測值,ki是理想輸出,N為樣本數(shù)目。

圖3 改進粒子群算法優(yōu)化過程

1.3 模型搭建

本文建模所用的實驗樣本數(shù)據(jù)采樣率為2M,選取6 000個點和“理想”階躍信號構建逆濾波模型,辨識逆向濾波器的階數(shù)和系數(shù)。圖4為殘差平方和隨階數(shù)變化示意圖,從圖中可以明顯看出殘差平方和在2階之后明顯下降,結合硬件的實現(xiàn)難易程度考慮選取6階。

圖4 殘差平方和隨階數(shù)變化示意圖

圖5 補償前后信號對比圖

在補償濾波器階數(shù)確定之后,將粒子群訓練次數(shù)設置為4 000次,初始粒子群個數(shù)為40個,數(shù)據(jù)采樣頻率為2M,長度為6 000個點,wstart=0.95,wend=0.4,進行訓練,得到補償濾波器的傳遞函數(shù)H(z)和補償前后信號對比圖,如圖5所示。從圖5可明顯看出補償后信號的超調量降到5%左右,上升時間3 μs,平臺壓力穩(wěn)定,補償效果良好。

H(z)=(0.214 8+0.392 2z-1+0.094 1z-2+0.375 0z-3+

0.205 1z-4+0.147 4z-5)/(2-0.320 3z-1-0.157 1z-2+

0.860 6z-3-0.357 4z-4-0.598 6z-5)

2 IIR數(shù)字濾波器量化精度設計

補償濾波器為IIR遞歸模型,這樣的結構在硬件實現(xiàn)過程中,需要考慮有限字長效應,保證量化后的系統(tǒng)極點都位于單位圓內(nèi),滿足設計的穩(wěn)定性[14-15]。

數(shù)字濾波器量化精度的設計很大程度上決定了濾波器的性能指標。理論上,數(shù)字濾波器的系數(shù)應當采用無限精度來表示,但實際中用硬件來實現(xiàn),將系數(shù)量化成有限位數(shù),將引起濾波器頻率特性的改變,甚至導致濾波器的不穩(wěn)定[16]。從量化角度考慮,系數(shù)的位數(shù)可以選擇足夠大來保證系統(tǒng)的穩(wěn)定性,但由于硬件條件的約束,系數(shù)位數(shù)不可能太大,需要通過仿真確定濾波器系數(shù)字長。此處將得到的系數(shù)統(tǒng)一量化以滿足后續(xù)計算要求,且保持系統(tǒng)的穩(wěn)定性。

本文將不同的系數(shù)量化位數(shù)和輸出量化位數(shù)進行比較,從圖6、圖7可以看出,系數(shù)及運算字長均會對濾波器的性能產(chǎn)生影響,綜合考慮資源占用率和數(shù)據(jù)精度,本文選取12 bit系數(shù),16 bit輸出量化來滿足數(shù)字濾波器的濾波效果以及頻響穩(wěn)定性。

圖6 IIR濾波器系數(shù)有限字長效應

圖7 IIR濾波器輸出數(shù)據(jù)有限字長效應

3 IIR數(shù)字濾波器軟核設計

3.1 IIR濾波器設計原理

根據(jù)已知,動態(tài)補償濾波器的離散傳遞函數(shù)為

(5)

b0y(m)+b1y(m-1)+…+bny(m-n)=a0x(m)+
a1x(m-1)+…+anx(m-n)

(6)

因IIR濾波器具有反饋結構,在實現(xiàn)零點系數(shù)及極點系數(shù)的運算時,應嚴格滿足時序要求,這一結構特點限制了系統(tǒng)的實現(xiàn)速度[17]。為了提高系統(tǒng)的運行速度,零極點系數(shù)的運算采用全并行結構,采用移位相加法實現(xiàn)常系數(shù)乘法運算。傳遞函數(shù)H(z)可看成由兩個FIR型濾波器構成,即式(7)、式(8)中的零點部分(zero parallel)和極點部分(pole parallel)兩個部分。

Zero(n)=a0x(m)+a1x(m-1)+…+anx(m-n)

(7)

Pole(n)=b1y(m-1)+…+bny(m-n)

(8)

b0y(m)=Zero(n)-Pole(n)=Ysum

(9)

式(7)、式(8)中:延遲數(shù)據(jù)y(m-n)和z(m-n)用寄存器保存,與常系數(shù)an和bn的乘法運算用左移替代,即將常系數(shù)分解成多個2的N次冪數(shù)相加的形式,將乘法轉換為移位及加法操作。同時用移位運算替代式(9)中除b0操作(此處量化后的b0=1 024,故通過右移10位來實現(xiàn)),得到當前時刻的輸出值z(m)。

整個IIR濾波器的閉環(huán)求取過程只在求取Zsum的減法器,以及移位運算來實現(xiàn)除法運算的過程中完成。濾波器在求取Zero(n)和Pole(n)、Zsum信號的過程中通過增加寄存器字長來實現(xiàn)全精度運算,出現(xiàn)運算誤差的環(huán)節(jié)只存在于除法運算(即z(m)的求取),此處通過仿真確定輸出數(shù)據(jù)位數(shù)精度,減小對輸出結果的影響。

3.2 仿真結果分析

將建模所用的實測數(shù)據(jù)保存成12位的二進制文本作為測試文件的輸入信號,用文件IO函數(shù)讀取,并保存輸出數(shù)據(jù)(數(shù)據(jù)長度16位的txt文件)。用MATLAB分析理想輸出數(shù)據(jù)和仿真數(shù)據(jù)曲線,可知擬合效果良好。

仿真圖8中,當rst復位信號無效,clk時鐘信號為上升沿開始采樣,din[11:0](AD采集到的傳感器數(shù)字信號)給xin[12:0]賦值。xin_reg[5:0]是零點部分Zero(n)輸入信號寄存器,將輸入信號xin[12:0]依次延遲保存下來。mult_reg[5:0]用來保存aix(m-i)值,xout[30:0]是零點部分計算結果Zero(n)。極點部分pole(n)同理,yin[15:0]保存最終計算結果y(m)。

圖8 IIR濾波器仿真圖

4 實測結果分析

沖擊波測試系統(tǒng)的電路設計框圖如圖9所示,在激勵信號的作用下,傳感器將測試信號轉換成電壓信號,通過調理信號電路進行濾波和放大,將輸入信號的電平控制在A/D的測試范圍內(nèi),在FPGA的控制下,A/D量化后的數(shù)據(jù)流經(jīng)過補償濾波器的處理,其結果保存到存儲器SDRAM中,通過USB接口和上位機(Labview平臺),讀取處理后的有效信號,進行下一步的效應分析。

圖9 設計框圖

用激波管對所設計的濾波器模型進行重復性動態(tài)測試,并對輸出的二進制文件讀取和分析得到對比圖,如圖10所示,得到補償前后的動態(tài)參數(shù)值,如表1所示。可明顯看出,該方法補償效果明顯,動態(tài)參數(shù)改善明顯,能為實際沖擊波測試提供準確數(shù)據(jù)。

圖10 補償前后信號對比圖

表1 動態(tài)補償前后參數(shù)對比

從實驗結果可以看出,該補償系統(tǒng)提高了系統(tǒng)的響應時間,降低了超調量,同時衰減度也近似約等于1,滿足系統(tǒng)穩(wěn)定性的要求,同時提高了系統(tǒng)的有效帶寬。從圖11中可看出測試系統(tǒng)由補償前87.1 kHz擴展到172.4 kHz,滿足沖擊波有效帶寬100 kHz(±3 dB)的要求,高頻處的噪聲影響沒有放大,也減小了諧振點對輸出信號的影響,動態(tài)特性補償效果明顯,滿足實際測試沖擊波的要求。

圖11 補償前后幅頻特性圖

5 結論

基于傳感器動態(tài)補償原理,利用激波管實驗數(shù)據(jù),結合QR分解和改進粒子群算法,設計了動態(tài)補償濾波器。根據(jù)仿真結果可以看出,該補償模型很好地修正了由于傳感器有效帶寬不足引起的動態(tài)誤差。在考慮數(shù)字濾波器的有限字長效應的基礎上,選取合適的系數(shù)和輸出位數(shù),同時補償濾波器硬件實現(xiàn)采用全并行單反饋結構,避免了設計復雜的反饋電路,保證運算精度的同時提高了運算速度,實現(xiàn)了對傳感器動態(tài)誤差實時修正的目的。實驗結果表明,該濾波器補償效果明顯且切實可行。

參考文獻:

[1] 何志文. 沖擊波超壓測試系統(tǒng)動態(tài)特性研究[D]. 太原:中北大學,2014.

[2] 劉一江,孟立凡,張志杰,等. 沖擊波測試系統(tǒng)中傳感器動態(tài)補償裝置[J]. 傳感技術學報,2012,(11):1516-1521.

[3] 童曉. 爆炸場沖擊波壓力測量及數(shù)據(jù)處理方法研究[D]. 南京:南京理工大學,2015.

[4] 賴富文,張志杰,張建宇,等. 基于動態(tài)特性補償?shù)臎_擊波測試數(shù)據(jù)處理方法[J]. 爆炸與沖擊,2015(6):871-875.

[5] 佘天莉. 測振傳感器的動態(tài)特性補償研究[D]. 中國地震局工程力學研究所,2006.

[6] 李悅,鐘新躍. 一種傳感器動態(tài)誤差實時修正方法及其FPGA實現(xiàn)[J]. 核電子學與探測技術,2012(9):1112-1116.

[7] 吳健,張志杰,王文廉. 傳感器動態(tài)誤差高速并行修正方法及其FPGA實現(xiàn)[J]. 傳感技術學報,2012,25(1):67-71.

[8] Rivera-Mejía J,Villafuerte-Arroyo J E,Vega-Pineda J,et al. Comparison of Compensation Algorithms for Smart Sensors with Approach to Real-Time or Dynamic Applications[J]. IEEE Sensors Journal,2015,15(12):7071-7080.

[9] 張海龍,劉一江,馬鐵華,等. 基于DSP和IIR的傳感器動態(tài)特性改善單元[J]. 傳感技術學報,2013,26(9):1254-1257.

[10] 軒春青,軒志偉,陳保立. 基于最小二乘與粒子群算法的壓力傳感器動態(tài)補償方法[J]. 傳感技術學報,2014,27(10):1363-1367.

[11] 黃俊欽. 測試系統(tǒng)動力學[M]. 北京:國防工業(yè)出版社,2013:12.

[12] 陳達榮. 基于粒子群優(yōu)化算法的自適應濾波器電路設計[D]. 西安:.西安電子科技大學,2013.

[13] 劉楊,田學鋒,詹志輝. 粒子群優(yōu)化算法慣量權重控制方法的研究[J]. 南京大學學報(自然科學版),2011,47(4):364-371.

[14] 羅海. 基于FPGA的高速IIR數(shù)字濾波器設計與實現(xiàn)[D]. 電子科技大學,2007.

[15] 曾菊容. 基于FPGA的IIR數(shù)字濾波器的設計與實現(xiàn)[D]. 西南交通大學,2008.

[16] Wisniewski M,Wcislik M. Digital Equalizer for Data Acquisition Path,Constructed Using IIR Filters[J]. IFAC-Papers on Line,2016,49(25):342-345.

[17] 杜勇.數(shù)字濾波器的MATLAB與FPGA實現(xiàn)[M]. 北京:電子工業(yè)出版社,2014:8.