基于MATLAB的PDV系統(tǒng)數(shù)據(jù)處理算法與實(shí)現(xiàn)

房文涵

(航空工業(yè)北京長(zhǎng)城計(jì)量測(cè)試技術(shù)研究所，北京 100095)

0 引言

光子多普勒測(cè)振儀是利用激光多普勒效應(yīng)對(duì)物體振動(dòng)進(jìn)行測(cè)量的一種測(cè)量?jī)x器。它與傳統(tǒng)激光測(cè)振儀相比，具有可以遠(yuǎn)距離測(cè)量、分辨力高、測(cè)量時(shí)間短、抗干擾能力強(qiáng)、響應(yīng)頻帶寬等優(yōu)點(diǎn)，在特殊結(jié)構(gòu)振動(dòng)物體和復(fù)雜環(huán)境的測(cè)量方面受到重視[1]。

理想的多普勒信號(hào)與被測(cè)物體的運(yùn)動(dòng)速度呈線性關(guān)系。光子多普勒振動(dòng)測(cè)量得到的信號(hào)，其幅度和頻率都隨時(shí)間變化而變化，是非理想的多普勒信號(hào)。這個(gè)原始信號(hào)受到很多因素的影響，例如光纖器件、被測(cè)物體等。所以要想從這些非理想信號(hào)中提取出有用的振動(dòng)信息，必須同時(shí)滿足實(shí)時(shí)性、精度、測(cè)速范圍的要求，這無(wú)疑對(duì)PDV系統(tǒng)的信號(hào)處理技術(shù)提出了很高的要求。

1 原理及算法

1.1 條紋法

當(dāng)運(yùn)動(dòng)物體和光源發(fā)生相對(duì)運(yùn)動(dòng)的時(shí)候，從運(yùn)動(dòng)物體表面反射回來(lái)的光會(huì)產(chǎn)生多普勒頻移。光源發(fā)出頻率為f1的光波，被測(cè)物體的運(yùn)動(dòng)速度為u，當(dāng)入射光垂直打在被測(cè)物體表面時(shí)，被測(cè)物體接受到的頻率為F1；接著被測(cè)物體當(dāng)作發(fā)射源，在它的垂直向有一個(gè)接收體，接收到的頻率為F2。此時(shí)被測(cè)物體的運(yùn)動(dòng)速度u表示為

(1)

(2)

只要能找出差頻頻率Δf，就可以計(jì)算出被測(cè)物體的運(yùn)動(dòng)速度u。差頻信號(hào)的每一個(gè)周期對(duì)應(yīng)一個(gè)條紋。找到條紋的相鄰極大值(或極小值)對(duì)應(yīng)的時(shí)間后，那么其時(shí)間差就是差頻周期Tbeat(t)，根據(jù)Δf=1/Tbeat(t)，那么差頻頻率就得到了。代入式(2)可得

(3)

這種方法就是條紋法，即通過(guò)計(jì)算條紋的周期來(lái)計(jì)算被測(cè)物體的運(yùn)動(dòng)速度，進(jìn)而求出被測(cè)物體的振動(dòng)信息。由以上公式可以看出，條紋法求被測(cè)物體運(yùn)動(dòng)速度的關(guān)鍵是要精確地求出每個(gè)信號(hào)極值點(diǎn)所對(duì)應(yīng)的時(shí)間。

1.2 PDV系統(tǒng)測(cè)振原理

光子多普勒速度測(cè)量系統(tǒng)原理如圖1所示。半導(dǎo)體激光器發(fā)出激光，進(jìn)入環(huán)形器的1端口，2端口與自聚焦棒相連作為檢測(cè)系統(tǒng)的探頭，3端口與探測(cè)器相連，并通過(guò)傳輸線與示波器相連，探測(cè)器是系統(tǒng)的光電轉(zhuǎn)化部分，示波器作為檢測(cè)結(jié)果的觀察顯示部分。激光在自聚焦棒的端面上，不可能完全透射，必將有一部分反射回來(lái)，將這路自聚焦棒的端面反射光稱為一次反射光，將透過(guò)自聚焦棒后打在靶體上反射回來(lái)的光，稱為二次反射光[6]。兩次反射的光將會(huì)發(fā)生干涉，并在探測(cè)器上實(shí)現(xiàn)光電轉(zhuǎn)化，轉(zhuǎn)化后的信號(hào)通過(guò)示波器進(jìn)行觀察。

圖1 PDV系統(tǒng)測(cè)振原理圖

在PDV系統(tǒng)中，利用探測(cè)器的平方檢測(cè)特性來(lái)提取差頻信號(hào)，探測(cè)器中的光強(qiáng)I為

(4)

式中：I1和I2分別為參考光和傳感光的光強(qiáng)，若采用交流探測(cè)器探測(cè)，從探測(cè)器輸出端到達(dá)示波器的電壓信號(hào)應(yīng)為

(5)

式中：R為探測(cè)器的電壓靈敏度。在式(5)中，對(duì)由運(yùn)動(dòng)引起的相位變化除以2π并取其整數(shù)部分，就得到示波器上顯示的周期數(shù)(條紋數(shù))N，即

(6)

式中：Int表示取整。對(duì)式(6)兩邊求導(dǎo)，可得物體的運(yùn)動(dòng)速度與條紋數(shù)隨時(shí)間變化率的關(guān)系為

(7)

式中：τ(t)為示波器顯示的干涉周期。實(shí)驗(yàn)分析得到的關(guān)系式(7)與理論分析中的式(3)對(duì)應(yīng)，因此只要確定τ(t)，就可以計(jì)算出物體的運(yùn)動(dòng)速度，進(jìn)而可計(jì)算出運(yùn)動(dòng)位移、加速度及振動(dòng)的頻率、振幅等物理量。

1.3 數(shù)據(jù)處理原理及算法實(shí)現(xiàn)

光子多普勒振動(dòng)測(cè)量得到的多普勒信號(hào)是一個(gè)相對(duì)比較復(fù)雜的信號(hào)，它并不是連續(xù)的信號(hào)，其具體波形取決于流速場(chǎng)中散射粒子的大小、濃度以及通過(guò)測(cè)量體的位置；多普勒信號(hào)的頻率一般都較高，通常在千赫茲的量級(jí)上；多普勒信號(hào)的信噪比隨測(cè)量對(duì)象而變化，速度越高，信噪比越低[2]。

光電探測(cè)器接收的是粒子散射光，其強(qiáng)度本身就比較微弱，而且信號(hào)受到諸如光路系統(tǒng)、雜散光、光檢測(cè)器件噪聲影響；當(dāng)用一般的光子多普勒測(cè)振儀檢測(cè)運(yùn)動(dòng)物體時(shí)，得到的原始多普勒信號(hào)都存在一定程度的噪聲及其他誤差干擾信號(hào)。誤差引起的干擾信號(hào)可根據(jù)理論公式等進(jìn)行部分補(bǔ)償修正[3]。而對(duì)于高頻噪聲和高頻信號(hào)相互混疊的含噪信號(hào)或者非平穩(wěn)信號(hào)的消噪，采用小波變換去噪會(huì)有很好的效果。這是因?yàn)?，噪聲是一種在時(shí)域和頻域上能量分布都比較均勻的信號(hào)，經(jīng)小波變換后其頻譜圖中不會(huì)出現(xiàn)特別大的突出峰；而多普勒信號(hào)的能量則非常集中，故在信號(hào)到達(dá)時(shí)刻和所在頻率段，信號(hào)的能量將有一個(gè)突變，在頻譜上有尖峰突起。

小波去噪算法可以多尺度的分解初始信號(hào)，準(zhǔn)確找到原始信號(hào)的高頻成分和低頻成分，減少對(duì)原始信號(hào)的影響，最大可能地去除噪聲[4]。對(duì)含噪聲的信號(hào)多尺度分解后，再對(duì)各個(gè)尺度上的小波系統(tǒng)進(jìn)行去噪，獲得小波逆變換重構(gòu)函數(shù)，最終得到去噪后的有效信號(hào)。

小波閾值去噪法是針對(duì)噪聲的小波系數(shù)幅值較小的特點(diǎn)，通過(guò)選擇合適的閾值來(lái)抑制噪聲，且保留反映原始的特征尖峰點(diǎn)。該方法計(jì)算量小、速度快、適用性廣。小波去噪一般有三種方法，分別是強(qiáng)制去噪聲法、默認(rèn)閾值去噪聲法、給定閾值去噪聲法[5]。其中，強(qiáng)制去噪聲法將小波分解結(jié)構(gòu)中的高頻系數(shù)全部置為0，濾掉所有高頻信號(hào)部分，然后對(duì)信號(hào)進(jìn)行小波重構(gòu)。這種方法較簡(jiǎn)單，且去噪后的信號(hào)曲線相對(duì)平滑，但容易丟失有用部分的信號(hào)；默認(rèn)閾值去噪聲法，是利用小波函數(shù)生成信號(hào)的默認(rèn)閥值，但這個(gè)閾值是不確定的；給定閾值去噪聲法，其閥值可在實(shí)際去噪過(guò)程中通過(guò)實(shí)驗(yàn)公式獲得，比默認(rèn)閥值的可信度高。

在此基礎(chǔ)上，本文提出了一種速度快、可靠性好的數(shù)據(jù)處理算法。下面對(duì)該算法進(jìn)行詳細(xì)描述：

首先對(duì)獲得的多普勒信號(hào)進(jìn)行移動(dòng)均值平滑處理，通過(guò)設(shè)置合適的平滑度參數(shù)，在保留有效的振動(dòng)信息的同時(shí)，濾掉噪聲信號(hào)。同時(shí)，標(biāo)出各個(gè)時(shí)刻平滑后信號(hào)的最大值，擬合成曲線，用不同顏色標(biāo)出。移動(dòng)均值濾波的原理是，設(shè)置平滑度參數(shù)k1，k1個(gè)數(shù)值一組組取平均值，多普勒信號(hào)中正反向曲線帶有的噪聲就可被消除。示波器采集到的原始信號(hào)見(jiàn)圖2。當(dāng)平滑度參數(shù)過(guò)大時(shí)，有效多普勒信號(hào)會(huì)被部分濾掉，即如圖3所示的密集曲線處顯示大片空白；當(dāng)平滑度參數(shù)過(guò)小時(shí)，多普勒信號(hào)中的大部分噪聲仍存在，即如圖4所示的多普勒曲線上仍疊加有尖刺。只有平滑度參數(shù)大小合適時(shí)，信號(hào)中的大部分噪聲被濾掉，且有效的信息仍被保留，各時(shí)刻信號(hào)最大值可連成平滑的半正弦弧線，如圖5中的綠色曲線所示。

然后，根據(jù)標(biāo)出的各個(gè)相鄰最大值所對(duì)應(yīng)的時(shí)間，求出差頻頻率Δf，代入式(2)得到初步的速度時(shí)間曲線，但在速度曲線的最低點(diǎn)處有少量的尖刺雜峰，見(jiàn)圖6。

圖2 示波器采集的原始信號(hào)圖

圖3 過(guò)大平滑度參數(shù)的濾波曲線圖

圖4 過(guò)小平滑度參數(shù)的濾波曲線圖

圖5 合適平滑度參數(shù)的濾波曲線圖

圖6 初步速度時(shí)間曲線

這時(shí)，需要調(diào)整極小限定值與極大限定值這兩個(gè)參數(shù)，找到雜峰的起始點(diǎn)和結(jié)束點(diǎn)，如圖7所示，將起始點(diǎn)標(biāo)為綠色，結(jié)束點(diǎn)標(biāo)為黑色。消去雜峰起始點(diǎn)(綠點(diǎn))到雜峰結(jié)束點(diǎn)(黑點(diǎn))之間的數(shù)據(jù)，即可去掉雜峰。

圖7 找出雜峰起始和結(jié)束點(diǎn)

去掉雜峰后，得到一系列不連續(xù)的半正弦曲線，將這些曲線進(jìn)行翻轉(zhuǎn)及最小二乘擬合處理，獲得連續(xù)的正弦速度時(shí)間曲線，見(jiàn)圖8，該曲線直觀顯示被測(cè)物體的振動(dòng)信息。

圖8 連續(xù)的正弦速度時(shí)間曲線

2 編程實(shí)現(xiàn)

MATLAB具有非常強(qiáng)的數(shù)據(jù)分析能力，可對(duì)數(shù)據(jù)進(jìn)行小波變換處理。該P(yáng)DV數(shù)據(jù)處理軟件操作便捷，其數(shù)據(jù)處理過(guò)程如圖9所示。

圖9 數(shù)據(jù)處理程序流程圖

1)數(shù)據(jù)讀取

數(shù)據(jù)讀取部分，利用示波器采樣，示波器的輸出文件為csv格式，則選用數(shù)據(jù)文件讀取函數(shù)csvread來(lái)讀取全部數(shù)據(jù)。

2)數(shù)據(jù)初步處理：濾波、平滑等

選用MATLAB中的smooth(I,k1,′loess′)給定閾值去噪聲法來(lái)進(jìn)行數(shù)據(jù)的平滑處理。smooth為移動(dòng)平均濾波，其默認(rèn)對(duì)數(shù)據(jù)周圍的5個(gè)數(shù)據(jù)點(diǎn)做移動(dòng)平均處理；k1是平滑度參數(shù)，可對(duì)指定k1個(gè)數(shù)據(jù)進(jìn)行移動(dòng)平均濾波；smooth有5個(gè)特定的平滑方法，loess方式是加權(quán)的線性最小二乘濾波方式。smooth平滑算法數(shù)據(jù)處理速度快，實(shí)時(shí)性好。通過(guò)設(shè)置合適的平滑度參數(shù)，在保留有效的振動(dòng)信息的同時(shí)，濾掉噪聲信號(hào)。

3)五點(diǎn)法處理得到原始速度曲線

經(jīng)過(guò)初步的平滑濾波后，得到原始速度曲線，只有速度最低點(diǎn)處仍存在雜峰。根據(jù)在每五個(gè)點(diǎn)中找到其中最大最小值，調(diào)整極小限定值k2與極大限定值k3這兩個(gè)參數(shù)，使極小限定值k2的值小于所有的最小值，極大限定值k3的值大于所有的最大值。此時(shí)，極小限定值k2和極大限定值k3的值即為雜峰的起始點(diǎn)和結(jié)束點(diǎn)[7]；將雜峰的起始點(diǎn)和結(jié)束點(diǎn)對(duì)應(yīng)的時(shí)間段清零，即可消去起始點(diǎn)和結(jié)束點(diǎn)之間的數(shù)據(jù)，去除雜峰，得到一系列不連續(xù)的半正弦曲線[8]。最后將這些曲線進(jìn)行翻轉(zhuǎn)及最小二乘擬合處理，獲得最終的連續(xù)正弦曲線。

4)速度加速度曲線等參量的顯示和存儲(chǔ)

利用函數(shù)uiputfile({′*.jpg′,′figure type(*.jpg)′}，可將得到的速度時(shí)間圖與加速度時(shí)間圖，存為jpg格式圖片。利用函數(shù)uiputfile({′*.xls′,′excel(*.xls)′}，可將得到的速度時(shí)間數(shù)據(jù)與加速度時(shí)間數(shù)據(jù)，存為xls格式文件。

3 實(shí)驗(yàn)驗(yàn)證

為了測(cè)試數(shù)據(jù)處理軟件的穩(wěn)定、可靠以及普適性，選取激振臺(tái)作為運(yùn)動(dòng)物體，進(jìn)行測(cè)試實(shí)驗(yàn)。振動(dòng)平臺(tái)由信號(hào)發(fā)生器、功率放大器還有激振臺(tái)組成，激振臺(tái)隨著信號(hào)的輸入發(fā)生振動(dòng)，振動(dòng)的頻率和信號(hào)的頻率相同，振動(dòng)的強(qiáng)度和信號(hào)的強(qiáng)度有關(guān)，即振動(dòng)強(qiáng)度和功率放大器所選擇的放大倍數(shù)有關(guān)。

信號(hào)發(fā)生器產(chǎn)生一個(gè)固定頻率的正弦信號(hào)，經(jīng)功率放大器給激振臺(tái)一個(gè)穩(wěn)定的正弦電壓信號(hào)驅(qū)動(dòng)，使鋼質(zhì)的振子做恒定上下往復(fù)運(yùn)動(dòng)。分別給激振臺(tái)20,80,160,1000 Hz的正弦信號(hào)，讓測(cè)速系統(tǒng)對(duì)其進(jìn)行檢測(cè)，通過(guò)示波器收集數(shù)據(jù)，用處理軟件處理數(shù)據(jù)；并且與經(jīng)校準(zhǔn)的激光測(cè)振儀的檢測(cè)結(jié)果進(jìn)行對(duì)比，見(jiàn)圖10。

圖10 標(biāo)準(zhǔn)激光干涉測(cè)速系統(tǒng)對(duì)比實(shí)驗(yàn)

具體操作為：PDV系統(tǒng)與經(jīng)校準(zhǔn)的激光測(cè)振儀分別對(duì)同一點(diǎn)進(jìn)行振動(dòng)測(cè)量，求經(jīng)該軟件處理得到的振動(dòng)頻率和標(biāo)準(zhǔn)激光干涉儀解調(diào)出的振動(dòng)頻率的相對(duì)偏差，是否在實(shí)驗(yàn)室規(guī)定偏差范圍內(nèi)，來(lái)驗(yàn)證數(shù)據(jù)處理軟件的穩(wěn)定、可靠以及普適性。

不同頻率輸入信號(hào)下的參數(shù)及結(jié)果如下表1所示。將經(jīng)校準(zhǔn)的激光測(cè)振儀解調(diào)出的振動(dòng)頻率設(shè)為標(biāo)準(zhǔn)頻率。光子多普勒測(cè)振系統(tǒng)的一般測(cè)量精度約為5%，經(jīng)校準(zhǔn)的激光測(cè)振儀精度為0.05%，則兩者對(duì)實(shí)驗(yàn)允差可看作5%。

表1 不同頻率下解調(diào)結(jié)果偏差

不同頻率輸入信號(hào)解調(diào)后得到的速度時(shí)間曲線如圖11所示。

圖11 不同頻率信號(hào)的速度曲線圖

根據(jù)表1中的實(shí)驗(yàn)數(shù)據(jù)及得到的速度時(shí)間曲線，可以發(fā)現(xiàn)：經(jīng)該數(shù)據(jù)處理軟件處理得到的振動(dòng)頻率與經(jīng)校準(zhǔn)的激光測(cè)振儀解調(diào)得到的振動(dòng)頻率的偏差很小，在對(duì)比實(shí)驗(yàn)的允差范圍內(nèi)。數(shù)據(jù)處理的頻率范圍是20 Hz～1 kHz。故本數(shù)據(jù)處理軟件算法穩(wěn)定、可靠、便捷且具有普適性。

4 總結(jié)

在條紋法基礎(chǔ)上，利用MATLAB小波函數(shù)smooth給定閾值去噪聲法實(shí)現(xiàn)了對(duì)PDV系統(tǒng)多普勒信號(hào)的去噪平滑處理，得到了理想的正弦信號(hào)曲線。經(jīng)過(guò)激振臺(tái)的測(cè)試試驗(yàn)，PDV數(shù)據(jù)處理軟件能正確地顯示并計(jì)算出物體速度變化，具有良好的可靠性。