蔣 強(qiáng)，俞 躍，葉凌偉，郝 元

〈測(cè)量技術(shù)〉

基于修正自相關(guān)算法的THz-TDS聚乙烯測(cè)厚反卷積算法

蔣 強(qiáng)1,2，俞 躍3，葉凌偉1,2，郝 元4

（1. 浙江省特種設(shè)備科學(xué)研究院，浙江 杭州 310018；2. 浙江省特種設(shè)備安全檢測(cè)技術(shù)研究重點(diǎn)實(shí)驗(yàn)室，浙江 杭州 3100183；3. 中國(guó)特種設(shè)備檢測(cè)研究院，北京 100029；4. 中國(guó)計(jì)量大學(xué)，浙江 杭州 310018）

在聚乙烯材料管道生產(chǎn)制造過程中，管道厚度精確測(cè)量是影響管道生產(chǎn)質(zhì)量的一個(gè)重要技術(shù)問題，太赫茲非金屬厚度透射測(cè)量過程中，通過測(cè)量?jī)纱瓮干湫盘?hào)的光程時(shí)間差與材料折射率計(jì)算得到樣品厚度參數(shù)，為了準(zhǔn)確提取兩次透射信號(hào)的光程時(shí)間差值，需要對(duì)原始時(shí)域信號(hào)進(jìn)行信號(hào)表征。通過高斯濾波反卷積得到脈沖響應(yīng)函數(shù)表征時(shí)域信號(hào)已經(jīng)成為了一種成熟且有效的技術(shù)手段。但一方面由于太赫茲波在聚乙烯材料中的傳播存在明顯的頻散現(xiàn)象，同時(shí)其透射率高回波信號(hào)較弱，導(dǎo)致信噪比較低，在脈沖響應(yīng)函數(shù)提取過程中放大了干擾信號(hào)的比重，造成信號(hào)失真、淹沒，無(wú)法準(zhǔn)確提取信號(hào)。本文通過對(duì)帶通濾波后的時(shí)域信號(hào)采用修正型自相關(guān)算法進(jìn)行自相關(guān)性分析，再結(jié)合高斯濾波反卷積運(yùn)算，得到改進(jìn)后的脈沖響應(yīng)函數(shù)。旨在解決太赫茲非金屬測(cè)量領(lǐng)域時(shí)域信號(hào)表征問題。改進(jìn)算法實(shí)驗(yàn)結(jié)果顯示加強(qiáng)了信號(hào)的清晰度，解決了因雜波信號(hào)干擾造成的脈沖響應(yīng)信號(hào)失真或淹沒情況，顯著提高了脈沖響應(yīng)信號(hào)信噪比。

太赫茲；反卷積；信號(hào)表征；自相關(guān)函數(shù)；信噪比

0 引言

近年來由于天然氣的廣泛使用，管道作為天然氣最重要的運(yùn)輸載體，被大量運(yùn)用，聚乙烯管道因制造與安裝成本低、耐腐蝕性強(qiáng)等優(yōu)良特性更是成為首選[1]。天然氣因其易燃易爆等特點(diǎn)極易產(chǎn)生安全隱患，因此對(duì)天然氣管道的質(zhì)量控制就顯得尤為重要。特別是管道的厚度參數(shù)控制，是行業(yè)內(nèi)一個(gè)難點(diǎn)。現(xiàn)有的厚度檢測(cè)手段主要為超聲波厚度檢測(cè)，而超聲波[2]由于聚乙烯材料隨環(huán)境（溫、濕度）變化性質(zhì)改變較大，導(dǎo)致超聲波在聚乙烯材料中傳播速度不同等原因，造成測(cè)量結(jié)果誤差較大，難以保證精確度，更無(wú)法作為質(zhì)量合格的檢驗(yàn)標(biāo)準(zhǔn)。由于太赫茲波對(duì)非金屬材料可穿透性較強(qiáng)，因此被大量應(yīng)用于非金屬材料厚度參數(shù)求取，為聚乙烯管道厚度檢測(cè)提供有力的補(bǔ)充。

Lionel Duvillaret利用透射樣品與樣品中反射的回波折射率差值求取厚度參數(shù)[3]。王秀敏提出一種誤差理論分析方法測(cè)定樣品厚度，確定厚度范圍再選用厚度參考值，計(jì)算所有頻率中的折射率，再重新選用厚度參考值重復(fù)以上工作，直到所有厚度計(jì)算完成，得出對(duì)厚度的誤差曲線，選用誤差最小的點(diǎn)，即為樣品厚度[4]。李麗娟使用太赫茲單點(diǎn)厚度提取模型，根據(jù)太赫茲傳播理論建立模型提出算法，測(cè)量出樣品的厚度，但僅針對(duì)膠層厚度進(jìn)行實(shí)驗(yàn)分析[5]。陸慶華通過對(duì)比，提出利用太赫茲光譜成像原理針對(duì)藥物的厚度進(jìn)行檢測(cè)[6]?？軐捲贚ionel Duvillaret的基礎(chǔ)上，通過分析太赫茲波在樣品不同頻率中的頻譜特性，分析出有效頻段與無(wú)效頻段，并且排除無(wú)效頻段來進(jìn)行厚度測(cè)量，提高了測(cè)量精度[7]。劉子燁提出使用太赫茲光子混頻連續(xù)波透射成像的方法，利用其相位信息實(shí)現(xiàn)對(duì)樣品厚度的測(cè)量[8]。

目前太赫茲非金屬材料厚度求取研究中多以厚度改進(jìn)算法為研究方向，測(cè)試材料也多以反射率較高的樣品為主，對(duì)于反射率較低的材料并不能準(zhǔn)確地通過兩次透射樣品的時(shí)域信號(hào)峰值計(jì)算得到光程時(shí)間差值，無(wú)法為厚度計(jì)算提供準(zhǔn)確光程時(shí)間差值參數(shù)。因此，太赫茲信號(hào)信噪比低已經(jīng)成為工程應(yīng)用中一個(gè)亟需解決的難題。

反卷積技術(shù)作為一種成熟有效的信號(hào)表征方法被國(guó)內(nèi)外學(xué)者廣泛應(yīng)用于信號(hào)處理問題。Zeiler[9-10]首先提出了反卷積（Deconvolution）的概念。蔡川等[11]通過反卷積技術(shù)對(duì)聚乙烯/納米銀復(fù)合材料空間電荷測(cè)量信號(hào)的恢復(fù)發(fā)現(xiàn)，高斯濾波器可以有效抑制恢復(fù)信號(hào)的震蕩，提高測(cè)量空間電荷的分辨率。吳超一等[12]采用反卷積算法對(duì)通過脈沖電聲法測(cè)量固體電介質(zhì)中空間電荷測(cè)量波形進(jìn)行了恢復(fù)，建立了波形衰減模型。Tanaka等[13]使用反卷積技術(shù)首次在界面電荷波形的恢復(fù)方面取得了成功。劉明亮等[14]簡(jiǎn)介了用卷積運(yùn)算實(shí)現(xiàn)反卷積的原理、方法和實(shí)例；其次，提出了該方法的病態(tài)問題，并給出了解決的辦法。邱程[15]等為實(shí)現(xiàn)超寬帶時(shí)域測(cè)量系統(tǒng)的超分辨測(cè)量，將反卷積算法應(yīng)用在測(cè)量系統(tǒng)數(shù)據(jù)處理過程中。根據(jù)測(cè)量系統(tǒng)中輸入脈沖及輸出信號(hào)的特點(diǎn)，選取高斯濾波器對(duì)反卷積處理帶來的病態(tài)問題進(jìn)行了修正，但是系統(tǒng)信噪比較低。Parruck等[16]發(fā)現(xiàn)在反卷積工程中因?yàn)橐恍﹨^(qū)域中的反卷積運(yùn)算誤差會(huì)使得()在計(jì)算過程中被放大，造成結(jié)果存在一系列類似函數(shù)的尖刺。這些尖刺幅值很大，在反Fast Fourier transform（FFT）變換后很可能成為時(shí)域信號(hào)的主要成分，從而掩蓋真實(shí)的輸入信號(hào)，造成頻域反卷積失敗。扈羅全等[17]使用基于Tikhonov正則化方法的反卷積算法得到了真實(shí)的空間電荷分布，利用數(shù)值實(shí)驗(yàn)研究了噪聲對(duì)反卷積算法的影響，結(jié)果表明，在無(wú)噪或者低噪環(huán)境下，反卷積算法能夠非常好地計(jì)算出電介質(zhì)中的空間電荷分布；在處理有噪數(shù)據(jù)時(shí)，反卷積的結(jié)果受到明顯的影響，但仍然有較高的計(jì)算精度。姚若河等[18]提出了時(shí)域卷積反演算法，不僅避免了傳統(tǒng)的傅里葉變換法所存在的缺點(diǎn)，且對(duì)于實(shí)際問題中由于輸入和輸出的起始元測(cè)量值往往幅度很小而信噪比也很小的情況，該算法能給出理想的結(jié)果。

自相關(guān)函數(shù)作為一種簡(jiǎn)單有效的信號(hào)提取方法也被大量應(yīng)用于信號(hào)處理。1953年貝爾實(shí)驗(yàn)室利用磁帶記錄儀最早實(shí)現(xiàn)了的相關(guān)檢測(cè)系統(tǒng)。高晉占等[19]描述了利用自相關(guān)從隨機(jī)噪聲中提取周期信號(hào)的方法。許峰等[20]提出一種基于相關(guān)檢測(cè)的電能質(zhì)量信號(hào)定位方法，能有效定位電能質(zhì)量信號(hào)發(fā)生的起止時(shí)間。胡大紅等[21]利用自相關(guān)函數(shù)識(shí)別時(shí)間序列變量的性質(zhì)和形態(tài)，檢驗(yàn)時(shí)間序列過程是否為零均值過程及其所建模型的適應(yīng)性。劉正平等[22]提出將自相關(guān)檢測(cè)與混沌理論相結(jié)合，檢測(cè)微弱周期脈沖信號(hào)的檢測(cè)方法。仿真結(jié)果表明在相同系統(tǒng)實(shí)驗(yàn)背景下，相結(jié)合的檢測(cè)方案優(yōu)于不加自相關(guān)的檢測(cè)方法。呂文秋等[23]從理論上證明了多重自相關(guān)處理后信號(hào)的信噪比，比相同取樣序列情況下相關(guān)處理的信噪比有所提高，比單次自相關(guān)具有更高的檢測(cè)能力，同時(shí)在獲得相同輸出信噪比時(shí)，重自相關(guān)檢測(cè)可以減小系統(tǒng)的復(fù)雜程度和運(yùn)算時(shí)間。范曉志等[24]提出，噪聲在自相關(guān)運(yùn)算中的變化是隨機(jī)的，因此，可以通過多次測(cè)量取平均值來消除。自相關(guān)算法有很多可以改進(jìn)的地方，這主要視方法所應(yīng)用的環(huán)境，若環(huán)境噪聲較強(qiáng)且復(fù)雜，則運(yùn)用多重自相關(guān)和求平均的方法可以明顯改善算法的效果，對(duì)于一般噪聲環(huán)境可直接利用自相關(guān)法自身的去噪能力。

天然氣管道因考慮安全因素，厚度參數(shù)較大，在太赫茲厚度測(cè)量過程中因太赫茲技術(shù)的局限，使得太赫茲?rùn)z測(cè)系統(tǒng)信噪比較低，回波接收信號(hào)弱，造成主觀選取回波時(shí)間差值的誤差增大，為厚度精確求取增加了難點(diǎn)。因此表征樣本時(shí)域信號(hào)，提高回波信號(hào)信噪比顯得尤為重要。本文使用了自相關(guān)函數(shù)對(duì)高斯濾波反卷積算法進(jìn)行改進(jìn)，以適應(yīng)太赫茲非金屬材料厚度測(cè)量過程中，時(shí)域信號(hào)信噪比較低造成的反卷積結(jié)果不理想等問題。同時(shí)對(duì)自相關(guān)函數(shù)進(jìn)行修正，以適應(yīng)太赫茲時(shí)域信號(hào)的特性。

1 實(shí)驗(yàn)原理

本文使用基于脈沖的太赫茲時(shí)域光譜技術(shù)對(duì)樣品（聚乙烯P100標(biāo)準(zhǔn)管道試快樣品）采用透射式檢測(cè)方法測(cè)量以獲取實(shí)驗(yàn)數(shù)據(jù)。如圖1所示，太赫茲波發(fā)射源產(chǎn)生太赫茲波，首先經(jīng)過首次透射在太赫茲波接收器接收到第一個(gè)脈沖信號(hào)，稱為主波信號(hào)。在太赫茲波首次透射經(jīng)過樣品的同時(shí)在樣品底面發(fā)生反射，太赫茲波在樣品底面與樣品表面經(jīng)過兩次反射后經(jīng)過第二次透射在太赫茲波接收端接收到第二次透射信號(hào)，稱為一次回波信號(hào)。在無(wú)樣品放入的情況下即測(cè)得唯一的主波信號(hào)，稱為空氣參考信號(hào)。通過準(zhǔn)確測(cè)量主脈沖信號(hào)與一次回波信號(hào)時(shí)間差值就可以精確的計(jì)算出樣品厚度。

由單點(diǎn)測(cè)厚模型可知：

式中：d為樣品厚度；n為材料的折射率；c為光在空氣中的傳播速度；DT為樣品主脈沖與一次回波信號(hào)的時(shí)間差值。

2 算法改進(jìn)原理

2.1 測(cè)量信噪比

如公式(2)所示，以一次回波信號(hào)開始與結(jié)束為兩個(gè)節(jié)點(diǎn)，分別向兩個(gè)節(jié)點(diǎn)左右擴(kuò)展100個(gè)脈沖時(shí)間點(diǎn)，共200個(gè)雜波信號(hào)求取均值，以一次回波峰峰值與雜波均值比值的最大值作為信噪比以及濾波效果最優(yōu)化的判別依據(jù)：

式中：SNR表示信噪比；signmm表示信號(hào)峰峰值；signl表示回波信號(hào)從信號(hào)峰峰值所在點(diǎn)向左取值；signr表示回波信號(hào)從信號(hào)峰峰值所在點(diǎn)向右取值；var表示取平均。

2.2 基于盲目反卷積的厚度測(cè)量

太赫茲時(shí)域測(cè)量系統(tǒng)可以簡(jiǎn)化為如圖2所示的信號(hào)模型。

圖2 系統(tǒng)輸入輸出與脈沖響應(yīng)的關(guān)系

如圖2所示，設(shè)()為系統(tǒng)的沖擊響應(yīng)函數(shù)，則有：

()＝()*() (3)

式中：()、()分別為系統(tǒng)時(shí)域輸入輸出信號(hào)。

()＝()() (4)

式中：()、()、()分別為()、()、()頻域變換；表示有效頻率。

()＝IFFT[()/()] (5)

式中：()、()分別為()、()頻域變換，()為系統(tǒng)的沖擊響應(yīng)函數(shù)，表示有效頻率。

由于系統(tǒng)屬于盲目反卷積系統(tǒng)，所以需要對(duì)輸入信號(hào)進(jìn)行假設(shè)。太赫茲時(shí)域系統(tǒng)是以皮秒激光器為輻射源信號(hào)，因此可以忽略信號(hào)在空氣中的損失。本文通過以空氣參考時(shí)域信號(hào)模擬系統(tǒng)輸入信號(hào)得到輸入信號(hào)()。由圖2可知，當(dāng)輸入信號(hào)的()信噪比較低時(shí)，()在有效信號(hào)區(qū)域的信號(hào)強(qiáng)度會(huì)因?yàn)榫垡蚁┎牧蠈?duì)輻射的吸收與反射等造成信號(hào)回波較弱，電場(chǎng)強(qiáng)度很小甚至為零，通過運(yùn)算得到的()就會(huì)很小，使信號(hào)淹沒在雜波信號(hào)中無(wú)法識(shí)別。要利用反卷積算法來提高系統(tǒng)的分辨率，首先必須解決信號(hào)的信噪比較低帶來的信號(hào)損失問題。

2.3 基于帶通濾波的信噪比改進(jìn)方法

由式(5)可知，當(dāng)輸入信號(hào)的()帶寬有限，()在一個(gè)或多個(gè)頻率區(qū)域的幅值很小甚至為零，通過運(yùn)算得到的()就會(huì)突然增大，成為一系列的尖刺。這些尖刺幅值很大，在逆傅里葉變換后很可能成為時(shí)域信號(hào)的主要成分，從而掩蓋真實(shí)的輸入信號(hào)，造成頻域反卷積失效。同樣如果在某些頻點(diǎn)上，輸入和輸出信號(hào)都很小甚至為零，那么得到的結(jié)果就完全是一個(gè)不確定的噪聲信號(hào)。這也就是反卷積運(yùn)算普遍存在的病態(tài)性的問題。

要利用反卷積算法來提高系統(tǒng)的分辨率，其次必須解決反卷積算法帶來的病態(tài)問題。

太赫茲時(shí)域測(cè)量系統(tǒng)的傳遞函數(shù)主要呈現(xiàn)低通特性，為了消除()的尖刺，需要設(shè)計(jì)一種低通濾波器()來對(duì)其進(jìn)行處理，從而得到穩(wěn)定、平滑的估值。由于測(cè)量系統(tǒng)雜波信號(hào)為高斯白噪聲信號(hào)，而且高斯信號(hào)在頻譜上表現(xiàn)和低通濾波器類似，所以在反卷積過程中采用合適的高斯濾波器可以有效地保留沖激響應(yīng)信號(hào)同時(shí)去除高頻的尖刺。

反卷積算法高斯濾波器濾波處理能夠有效剔除顯著雜波干擾信號(hào)，這使得信號(hào)變得平滑。但無(wú)法解決信號(hào)比較低帶來的信號(hào)損失問題，無(wú)法提高信號(hào)信噪比。在太赫茲聚乙烯材料厚度測(cè)量過程中，由于硬件系統(tǒng)信噪比有效區(qū)間的限制，在太赫茲波區(qū)間內(nèi)（0.1THz～10THz）存在大量信噪比較低頻段所帶來的雜波干擾信號(hào)，這使得信號(hào)本身帶有大量的無(wú)效頻段太赫茲信號(hào)，為了濾除無(wú)效的信號(hào)頻率區(qū)間，需要對(duì)信號(hào)進(jìn)行預(yù)處理操作。

運(yùn)用巴特沃斯濾波器帶通特性對(duì)信號(hào)步進(jìn)式帶通濾波，先確定帶通濾波器上下限，采取穩(wěn)定一邊，另一邊0.01THz的步進(jìn)值進(jìn)行信噪比計(jì)算，直到操作邊出現(xiàn)信噪比轉(zhuǎn)折點(diǎn)即確定濾波上截止頻率（下截止頻率），重復(fù)操作另一邊，即可得出相對(duì)最優(yōu)化濾波區(qū)間。濾除無(wú)效頻率區(qū)間，得到較理想的時(shí)域信號(hào)。結(jié)果顯示巴特沃斯濾波器能夠有效濾除無(wú)效頻段太赫茲信號(hào)，顯著提高信號(hào)信噪比：

(＋0)＝(＋0)(＋0) (6)

式中：0表示無(wú)效頻率。

￠(＋0)＝(＋0)()＝()×＋(0)×(9)

式中：代表帶通濾波器；、表示權(quán)值。

2.4 修正自相關(guān)算法

自相關(guān)算法是一種衡量單一信號(hào)時(shí)域序列在不同時(shí)間點(diǎn)的相關(guān)程度的工具。通過與移位之后的自身信號(hào)遍歷相乘即可得出各個(gè)時(shí)域信號(hào)之間的相關(guān)系數(shù)?？梢愿忧逦崛‰[藏或淹沒在雜波信號(hào)中的有效信號(hào)，提高信號(hào)信噪比。根據(jù)相關(guān)系數(shù)定義可以得出自相關(guān)函數(shù)在時(shí)間序列中的表達(dá)式：

式中：()為自相關(guān)系數(shù)；12為時(shí)間序列方差；為序列均值。

假設(shè)信號(hào)的復(fù)信號(hào)形式為：

()＝ej(2p＋)＋()＝1,2,3… (11)

式中：為信號(hào)頻率；、分別為幅度和初相；()為高斯白噪聲序列，其均值為0，方差為22。對(duì)()進(jìn)行自相關(guān)運(yùn)算，得：

式中：為序列長(zhǎng)度；*()為對(duì)輸入信號(hào)共軛；￠()為含有噪聲項(xiàng)的自相關(guān)運(yùn)算結(jié)果，其表達(dá)式為：

當(dāng)較大時(shí)，根據(jù)概率論中的中心極限定理，￠()可以看作近似高斯分布，其均值為0，方差等于(222s22＋s24)/因此，()的信噪比為：

從式(14)可以看出，相比于原信號(hào)，信噪比提高了/(2＋22/2)倍，相關(guān)輸出()的信噪比與序列長(zhǎng)度成正比與2/成反比。2/越小代表雜波信號(hào)越平穩(wěn)，則復(fù)合信號(hào)越平穩(wěn)。

算法1為偽代碼實(shí)現(xiàn)，其中第1行包括各類常量定義、變量定義與賦初值操作，第2～28行遍歷各個(gè)待布局邏輯功能，將各項(xiàng)約束寫入模型，第29行定義目標(biāo)函數(shù)，最后調(diào)用求解器對(duì)模型進(jìn)行求優(yōu).由表1可知，P=2為常數(shù)，同時(shí)對(duì)于目前FPGA芯片規(guī)模來說F，W，H通常也很小[8]，因此建立布局模型的時(shí)間復(fù)雜度為O(N2).由于MILP求優(yōu)的時(shí)間復(fù)雜度為指數(shù)階，因而UPRFloor模型的求優(yōu)時(shí)間復(fù)雜度和文獻(xiàn)[10]均為O(2N).

自相關(guān)函數(shù)的使用前提是平穩(wěn)序列，即在統(tǒng)計(jì)學(xué)自相關(guān)函數(shù)中要求期望值為0，方差為1，且方差和期望不隨時(shí)間的變化而變化，而太赫茲時(shí)域信號(hào)序列不能等同于平穩(wěn)信號(hào)序列，所以應(yīng)用傳統(tǒng)自相關(guān)函數(shù)計(jì)算自相關(guān)系數(shù)時(shí)并不能準(zhǔn)確表征原始信號(hào)。由式(15)可知，在自相關(guān)算方法中，由于補(bǔ)零操作移位之后的時(shí)間序列的均值會(huì)逐漸縮小并趨近于0，也隨著減小并趨近于零。時(shí)間序列在不斷的補(bǔ)零操作下會(huì)成為絕對(duì)意義上的零均值序列，因此只有時(shí)變的序列均值可以反映序列真實(shí)方差的變化。只有時(shí)變的方差才能反映移位后離散序列的離散程度。因此，結(jié)合太赫茲時(shí)域信號(hào)隨時(shí)間變量改變對(duì)應(yīng)期望與方差也隨之改變的特點(diǎn)，進(jìn)而采用修正自相關(guān)公式進(jìn)行處理：

修正自相關(guān)公式如下：

式中：()為修正自相關(guān)系數(shù)；1、2分別為原序列與延時(shí)后時(shí)間序列標(biāo)準(zhǔn)差；1、2分別為原序列與延時(shí)后時(shí)間序列均值。

對(duì)比公式(10)(16)，可以發(fā)現(xiàn)自相關(guān)算法（公式(10)）使用了同一個(gè)均值與方差對(duì)移位前后的時(shí)域信號(hào)序列進(jìn)行處理，修正自相關(guān)的算法（公式(16)）使用各自的均值與方差對(duì)公式進(jìn)行表達(dá)。隨著自相關(guān)序列的延遲時(shí)間增加，脈沖函數(shù)的的主信號(hào)會(huì)被后面的遞補(bǔ)雜波信號(hào)取代為信號(hào)序列，由于雜波信號(hào)的均值為0，時(shí)移序列的均值也會(huì)相應(yīng)減小，因此當(dāng)序列長(zhǎng)度N已經(jīng)選定時(shí)，信號(hào)的信噪比也會(huì)相應(yīng)提高。

綜上，本文改進(jìn)反卷積算法如圖3所示。

1）通過實(shí)驗(yàn)獲取非金屬太赫茲時(shí)域光譜系統(tǒng)測(cè)量時(shí)間序列，以空氣參考信號(hào)為基準(zhǔn)對(duì)空氣參考信號(hào)和樣本信號(hào)進(jìn)行歸一化處理；

2）以空氣參考時(shí)域信號(hào)模擬系統(tǒng)輸入信號(hào)得到輸入信號(hào)()，太赫茲時(shí)域光譜系統(tǒng)測(cè)得樣本時(shí)域輸出信號(hào)()。

3）通過經(jīng)典濾波器（巴特沃斯濾波器）對(duì)樣本信號(hào)進(jìn)行濾波預(yù)處理，得到最優(yōu)化濾波樣本時(shí)域信號(hào)。

4）分別對(duì)空氣參考信號(hào)與最優(yōu)化后的時(shí)域信號(hào)進(jìn)行自相關(guān)性分析；得到各自的自相關(guān)系數(shù)時(shí)間序列。

5）通過對(duì)空氣參考自相關(guān)系數(shù)時(shí)間序列和樣本信號(hào)自相關(guān)系數(shù)時(shí)間序列進(jìn)行高斯濾波反卷積運(yùn)算，提取出信號(hào)的脈沖響應(yīng)函數(shù)。

圖3 反卷積算法改進(jìn)流程圖

3 實(shí)驗(yàn)與結(jié)果分析

3.1 實(shí)驗(yàn)裝置

如圖4所示，實(shí)驗(yàn)采用太赫茲時(shí)域光譜（THz-TDS）系統(tǒng)，飛秒激光設(shè)備參數(shù)為：波長(zhǎng)800nm、脈寬100fs、頻率80MHz。系統(tǒng)原理為鈦寶石飛秒激光器經(jīng)過分束鏡將飛秒激光分為泵浦光、探測(cè)光。泵浦光通過可變延遲線入射到光導(dǎo)天線晶體上激發(fā)出太赫茲脈沖；探測(cè)光經(jīng)過多次反射后與太赫茲脈沖一同通過探測(cè)晶體，用于探測(cè)出太赫茲脈沖的電場(chǎng)強(qiáng)度隨時(shí)間變化情況。實(shí)驗(yàn)環(huán)境溫度為25℃（室溫），環(huán)境濕度為30%。

3.2 算法改進(jìn)實(shí)驗(yàn)

本文實(shí)驗(yàn)數(shù)據(jù)來源于聚乙烯P100樣品測(cè)試結(jié)果，其中時(shí)域信號(hào)取1200組數(shù)據(jù)，自相關(guān)系數(shù)圖取半圖譜，即向右平移1024個(gè)單位，脈沖響應(yīng)函數(shù)與自相關(guān)系數(shù)圖對(duì)應(yīng)取1024組。時(shí)域圖橫坐標(biāo)為歸一化時(shí)間，單位為0.1ps，縱坐標(biāo)為歸一化電場(chǎng)強(qiáng)度，單位為V，自相關(guān)系數(shù)圖橫坐標(biāo)為歸一化時(shí)間，單位為0.1ps，縱坐標(biāo)為自相關(guān)系數(shù)，單位為常數(shù)，脈沖響應(yīng)函數(shù)圖橫坐標(biāo)為歸一化時(shí)間，單位為0.1ps，縱坐標(biāo)歸一化電場(chǎng)強(qiáng)度，單位為V。為了統(tǒng)一化圖示，便于觀測(cè)，橫坐標(biāo)統(tǒng)一標(biāo)示為time/ps－1，縱坐標(biāo)統(tǒng)一標(biāo)示為power。

已知A為8mm厚度聚乙烯樣品，圖5(a)、(b)分別為樣品濾波信噪比值分布圖、濾波器頻域響應(yīng)圖。結(jié)合表1所示通過濾波頻段選取方法得出A樣品的最佳濾波區(qū)間為0.84～1.09THz，其中基準(zhǔn)濾波信噪比值為138.5674。0.84THz頻率處的濾波信噪比值144.5，1.09THz頻率處的濾波信噪比值達(dá)到最大為533.2。

圖6所示為空氣參考信號(hào)與樣本信號(hào)的時(shí)域圖，圖中所示實(shí)線為樣本信號(hào)，虛線為空氣參考信號(hào)，標(biāo)注點(diǎn)為樣本一次回波信號(hào)所在。由表1所示可得樣本原始信號(hào)回波峰峰值信噪比為4.7525，信噪比較低。不易從圖中直接獲得樣本信號(hào)主回波時(shí)間差值。由于采用透射式測(cè)量方法，主波峰值滯后142ps，其中參考信號(hào)峰值在125ps處，樣品時(shí)域信號(hào)主波峰值位于267ps處，樣品回波信號(hào)峰值點(diǎn)滯后主波峰值點(diǎn)80.2ps，位于1070ps處。

圖7所示為巴特沃斯濾波處理后的樣品時(shí)域信號(hào)，圖中標(biāo)注為樣本一次回波峰值點(diǎn)，由表2所示可得樣本巴特沃斯濾波后信號(hào)回波峰峰值信噪比為15.2329，是原始信號(hào)回波信噪比的3.21倍。由此可以說明巴特沃斯濾波器對(duì)回波信號(hào)具有顯著提高信號(hào)信噪比的作用。

圖8、9、10所示分別為參考信號(hào)、樣本信號(hào)自相關(guān)系數(shù)圖、樣本信號(hào)巴特沃斯濾波處理后自相關(guān)系數(shù)圖，如圖9所示自相關(guān)系數(shù)在80.2ps處最大，與圖1中樣本時(shí)域信號(hào)主回波峰值時(shí)間差值相吻合，即在右移802ps個(gè)時(shí)間單位后，出現(xiàn)回波自相關(guān)系數(shù)峰值點(diǎn)。由表2所示可得樣本自相關(guān)系數(shù)信噪比為5.9462，是原始信號(hào)回波信噪比的1.25倍。信噪比提高不明顯，說明自相關(guān)系數(shù)無(wú)法明顯提高信號(hào)信噪比。由圖10與表2所示可得樣本巴特沃斯濾波后自相關(guān)系數(shù)信噪比為28.7159，是原始信號(hào)回波信噪比的6.04倍。信噪比提高明顯，再次說明巴特沃斯濾波器對(duì)信噪比提高效果顯著，同時(shí)也說明自相關(guān)系數(shù)在表征時(shí)域信號(hào)時(shí)存在大量雜波干擾，無(wú)法準(zhǔn)確反映信號(hào)特征。

圖11、12、13所示分別為空氣參考信號(hào)巴特沃斯濾波處理后修正自相關(guān)系數(shù)圖、樣本修正自相關(guān)系數(shù)圖、樣本信號(hào)巴特沃斯濾波處理后修正自相關(guān)系數(shù)圖，由圖10與圖11、圖12與表2中4、5、6所示可得樣本修正自相關(guān)系數(shù)信噪比為26.7989，樣本信號(hào)巴特沃斯濾波后自相關(guān)系數(shù)信噪比為28.7159，樣本巴特沃斯濾波后修正自相關(guān)系數(shù)信噪比為30.8711，考慮到過程中的誤差，三者信噪比基本一致，可以得出自相關(guān)系數(shù)在相對(duì)最優(yōu)化濾波之后的效果等同于修正自相關(guān)系數(shù)，而后者更能體現(xiàn)原始信號(hào)的特征。

圖4 太赫茲時(shí)域光譜檢測(cè)系統(tǒng)

圖5 有效頻段濾波

表1 信噪比比對(duì)結(jié)果

圖6 空氣參考信號(hào)與樣本信號(hào)時(shí)域信號(hào)

圖7 巴特沃斯濾波處理后的樣品時(shí)域信號(hào)

圖8 空氣參考信號(hào)自相關(guān)系數(shù)圖

圖9 樣本信號(hào)自相關(guān)系數(shù)圖

圖10 樣本信號(hào)巴特沃斯濾波處理后自相關(guān)系數(shù)圖

圖11 空氣參考信號(hào)巴特沃斯濾波處理后修正自相關(guān)系數(shù)圖

圖12 樣本信號(hào)修正自相關(guān)系數(shù)圖

圖13 樣本信號(hào)巴特沃斯濾波處理后修正自相關(guān)系數(shù)圖

表2 信噪比比對(duì)結(jié)果

圖14、15、16所示分別為脈沖響應(yīng)函數(shù)（未處理）、脈沖響應(yīng)函數(shù)（濾波后自相關(guān)算法處理）、脈沖響應(yīng)函數(shù)（濾波后修正自相關(guān)算法處理）。表2中7、8、9所示依次對(duì)應(yīng)3種方法的處理結(jié)果。信噪比分別為2.9827、20.9324、26.7791。分別為原始信號(hào)信噪比的0.63倍、4.40倍、5.63倍。由結(jié)果所示，直接通過高斯濾波反卷積算法得到的脈沖響應(yīng)函數(shù)的結(jié)果顯示雜波干擾無(wú)法有效去除，無(wú)法直接適用于太赫茲非金屬反射式測(cè)量方式，經(jīng)有效濾波后自相關(guān)算法處理后的脈沖響應(yīng)函數(shù)與經(jīng)有效濾波后修正自相關(guān)算法處理后的脈沖響應(yīng)函數(shù)均能達(dá)到提高信噪比、表征時(shí)域信號(hào)的目的，但后者效果更加突出，更能體現(xiàn)時(shí)域信號(hào)的特征。

圖14 脈沖響應(yīng)函數(shù)（未處理）

圖15 脈沖響應(yīng)函數(shù)（濾波后自相關(guān)算法處理）

圖16 脈沖響應(yīng)函數(shù)（濾波后修正自相關(guān)算法處理）

4 結(jié)論

本文為解決太赫茲非金屬材料測(cè)量領(lǐng)域時(shí)域信號(hào)表征問題。采用基于修正自相關(guān)算法的THz-TDS聚乙烯測(cè)厚反卷積算法對(duì)太赫茲聚乙烯材料厚度測(cè)量時(shí)域數(shù)據(jù)進(jìn)行信號(hào)處理。本算法通過對(duì)帶通濾波后的時(shí)域信號(hào)采用自相關(guān)算法進(jìn)行自相關(guān)性分析，再結(jié)合高斯濾波反卷積運(yùn)算，得到改進(jìn)后的脈沖響應(yīng)函數(shù)。改進(jìn)算法實(shí)驗(yàn)結(jié)果顯示加強(qiáng)了信號(hào)的清晰度，解決了因雜波信號(hào)干擾造成的脈沖響應(yīng)信號(hào)失真或淹沒情況，顯著提高了脈沖響應(yīng)信號(hào)信噪比。有效解決了高斯反卷積算法直接應(yīng)用于非金屬太赫茲測(cè)量時(shí)域信號(hào)處理無(wú)法有效提取信號(hào)時(shí)域特征信息的問題。另外，提出一種自相關(guān)函數(shù)的修正算法，通過與自相關(guān)算法進(jìn)行結(jié)果比對(duì)，明顯可得自相關(guān)修正算法效果優(yōu)于自相關(guān)算法，更適合時(shí)域信號(hào)的信號(hào)特征。且在實(shí)驗(yàn)過程中發(fā)現(xiàn)修正自相關(guān)函數(shù)出現(xiàn)信號(hào)拖尾性，需要進(jìn)一步研究。

[1] 劉慶兆. 脈沖輻射場(chǎng)診斷技術(shù)[M]. 北京: 科學(xué)出版社, 1994.

LIU Qingzhao.[M]. Beijing: Science Press, 1994.

[2] John G P, Dimitris G M. 數(shù)字信號(hào)處理: 原理、算法與應(yīng)用[M]. 張曉林, 譯. 北京: 電子工業(yè)出版社, 2004 .

John G P, Dimitris G M.[M]. Translated by ZHANG Xiaolin. Beijing: Electronics Industry Press, 2004.

[3] Duvillaret L, Garet F, Coutaz J. Highly precise determination of optical constants and sample thickness in terahertz time-domain spectroscopy[J]..., 1999, 38(2): 409-415.

[4] 王秀敏. 太赫茲光譜分析中用誤差理論確定樣品厚度的研究[J]. 激光與紅外, 2009, 39(8): 884-886.

WANG Xiumin. Study on the determination of sample thickness using error theory in terahertz spectrum analysis[J]., 2009, 39(8): 884-886.

[5] 李麗娟, 周明星, 任姣姣. 基于太赫茲時(shí)域光譜的膠層厚度均勻性檢測(cè)[J]. 激光與紅外, 2014, 44(7): 801-804.

LI Lijuan, ZHOU Mingxing, REN Ye. Detection of uniformity of adhesive layer thickness based on terahertz time domain spectroscopy[J]., 2014, 44(7): 801-804.

[6] 陸慶華, 陳玉潔, 嚴(yán)盈富. 薄膜包衣厚度測(cè)量方法分析[J]. 南昌航空大學(xué)學(xué)報(bào), 2014, 28(4): 76-82.

LU Qinghua, CHEN Yujie, YAN Yingfu. Analysis of measurement methods of film coating thickness[J]., 2014, 28(4): 76-82.

[7] 寇寬, 趙國(guó)忠, 劉英. 利用太赫茲時(shí)域光譜同時(shí)確定樣品厚度和折射率[J]. 中國(guó)激光, 2015, 42(8): 0815001.1-081500.6.

KOU Kuan, ZHAO Guozhong, LIU Ying. Simultaneous determination of sample thickness and refractive index using terahertz time domain spectroscopy[J]., 2015, 42(8): 0815001.1-081500.6.

[8] 劉子燁, 劉建軍, 洪治. 基于光子混頻連續(xù)太赫茲波的厚度檢測(cè)[J]. 激光技術(shù), 2016, 40(4): 496-499.

LIU Zizhen, LIU Jianjun, HONG Zhi. Thickness detection of continuous terahertz waves based on photon mixing[J]., 2016, 40(4): 496-499.

[9] Zeiler MD, Krishnan D, Taylor GW, et al. Deconvolutional networks[C]//(),, 2010: 2528-2535.

[10] Zeiler M D, Taylor G W, Fergus R. Adaptive deconvolutional networks for mid and high level feature learning[C]//(),, 2011: 2018-2025.

[11] 蔡川, 李旭光, 尹毅, 等. 空間電荷測(cè)量信號(hào)恢復(fù)中的頻域反卷積技術(shù)及實(shí)現(xiàn)[J]. 電工技術(shù)學(xué)報(bào), 2009, 24(10): 165-169.

CAI Chuan, LI Xuguang, YIN Yi, et al. Frequency domain deconvolution technology and implementation in space charge measurement signal recovery[J]., 2009, 24(10): 165-169.

[12] 吳超一, 鐘力生, 王霞, 等. 脈沖電聲法空間電荷測(cè)量波形恢復(fù)的數(shù)據(jù)處理[J]. 中國(guó)電機(jī)工程學(xué)報(bào), 2005(16): 137-140.

WU Chaoyi, ZHONG Lisheng, WANG Xia, et al. Data Processing of Waveform Recovery of Pulsed Electroacoustic Space Charge Measurement[J]., 2005(16): 137-140.

[13] Tanaka Y, Hanawa K, Suzuki K, et al. Attenuation recovery technique for acoustic wave propagation in PEA method[C]//,2001, 2001: 407-410.

[14] 劉明亮, 蔡永泉, 饒敏, 等. 用卷積運(yùn)算實(shí)現(xiàn)反卷積[J]. 電子學(xué)報(bào), 2000(5): 111-112.

LIU Mingliang, CAI Yongquan, RAO Min, et al. Deconvolution with Convolution Operation[J]., 2000(5): 111-112.

[15] 邱程, 付紅衛(wèi), 王欣, 等. 反卷積算法在超寬帶時(shí)域測(cè)量中的應(yīng)用[J]. 探測(cè)與控制學(xué)報(bào), 2013, 35(3): 80-84.

QIU Cheng, FU Hongwei, WANG Xin, et al. Application of deconvolution algorithm in ultra-wideband time domain measurement[J]., 2013, 35(3): 80-84.

[16] Parruck B, Riad S M. An optimization criterion for iterative deconvolution[J]., 1983, 32(1): 137-140.

[17] 扈羅全, 鄭飛虎, 張冶文. 用于電介質(zhì)中空間電荷分布測(cè)量的Tikhonov反卷積算法[J]. 計(jì)算物理, 2004(5): 432-438.

LUO Luoquan, ZHENG Feihu, ZHANG Yewen. Tikhonov deconvolution algorithm for measuring space charge distribution in dielectrics[J]., 2004(5): 432-438.

[18] 姚若河. 時(shí)域卷積反演的一個(gè)新算法[J]. 桂林電子工業(yè)學(xué)院學(xué)報(bào), 1988(2): 38-42.

YAO Ruohe. A New Algorithm for Time Domain Convolution Inversion[J]., 1988(2): 38-42.

[19] 高晉占. 微弱信號(hào)檢測(cè)[M]. 北京: 清華大學(xué)出版社, 2004.

GAO Jinzhan.[M]. Beijing: Tsinghua University Press, 2004.

[20] 許峰, 李開成, 王可. 一種基于相關(guān)算法的電能質(zhì)量信號(hào)定位方法[J].電測(cè)與儀表, 2015, 52(3):1-5.

XU Feng, LI Kaicheng, WANG Ke. A Power Quality Signal Positioning Method Based on Correlation Algorithm[J]., 2015, 52(3): 1-5.

[21] 胡大紅, 姚志鵬. 自相關(guān)函數(shù)在時(shí)間序列模型識(shí)別中的應(yīng)用[J]. 襄樊學(xué)院學(xué)報(bào), 2011, 32(5): 15-18.

HU Dahong, YAO Zhipeng. Application of autocorrelation function in time series model identification[J]., 2011, 32(5): 15-18.

[22] 劉正平, 臧觀建, 陳俊杰. 微弱周期脈沖信號(hào)的自相關(guān)-混沌系統(tǒng)聯(lián)合檢測(cè)方法[J]. 噪聲與振動(dòng)控制, 2007(5): 63-65.

LIU Zhengping, ZANG Guanjian, CHEN Junjie. Auto-correlation-chaotic system detection method for weak periodic pulse signal[J]., 2007(5): 63-65.

[23] 呂文秋. 強(qiáng)噪聲中周期信號(hào)的數(shù)字重相關(guān)檢測(cè)[J]. 廈門大學(xué)學(xué)報(bào): 自然科學(xué)版, 1990, 29(1): 34-38．

LYU Wenqiu. Digital recorrelation detection of periodic signals in strong noise[J].: Natural Science Edition, 1990, 29(1): 34-38.

[24] 范曉志, 王長(zhǎng)廣, 黃曉紅, 等. 基于小波分析和多重自相關(guān)法的微弱信號(hào)檢測(cè)技術(shù)[J]. 計(jì)算機(jī)應(yīng)用與軟件, 2007, 24(5):40-41, 44.

FAN Xiaozhi, WANG Changguang, HUANG Xiaohong, et al. Weak signal detection technology based on wavelet analysis and multiple autocorrelation methods[J]., 2007, 24(5): 40-41.

Deconvolution Algorithm of THz-TDS Polyethylene Thickness Measurement Based on Modified Autocorrelation Algorithm

JIANG Qiang1,2，YU Yue3，YE Lingwei1,2，HAO Yuan4

(1. Zhejiang Special Equipment Science Research Institute, Hangzhou 341000, China; 2. Zhejiang Key Laboratory of Special Equipment Safety Testing Technology Research, Hangzhou 341000, China;3. China Special Equipment Inspection Institute, Beijing 100029, China; 4.China Measurement University, Zhejiang 310018, China)

In polyethylene pipe manufacturing processes, the precise measurement of the pipe thickness is an important technical problem, which affects the quality of the pipes produced. In the terahertz non-metal thickness transmission measurement process, the sample thickness parameter is obtained by measuring the optical path time difference of the two transmitted signals and material refractive index. To accurately extract the optical path time difference of the two transmitted signals, signal characterization of the original time domain signal is required, which is a mature and effective technical means to characterize the signals impulse response functions. An impulse response function can be obtained using a deconvolution technique. However, when the echo signal is weak or the system signal-to-noise ratio is low, the interference signal increases, causing distortion and flooding of the signal, and thus, the signal cannot be accurately extracted. In this study, we first pre-process the time domain signal with a Butterworth filter and perform an improved self-precision analysis on the processed signal. Second, we perform Gaussian deconvolution processing for autocorrelation signals. Finally, we obtain an improved impulse response function. The improved method aims to solve the problem of signal characterization in the field of terahertz non-metal measurement. The improved algorithm enhances the clarity of the signal and solves the distortion or flooding of the impulse response signal caused by the interference of the clutter signal. Consequently, the signal-to-noise ratio of the impulse response signal is improved.

terahertz, Gaussian filter deconvolution, signal characterization, improved autocorrelation function, signal to noise ratio

TB324

A

1001-8891(2020)05-0473-10

2019-04-14；

2020-04-02.

蔣強(qiáng)（1992-），男，碩士，主要從事太赫茲無(wú)損檢測(cè)技術(shù)與信號(hào)處理。E-mail: 1125138001@qq.com。

國(guó)家重點(diǎn)研發(fā)計(jì)劃項(xiàng)目（2017YFF0209704）。