段麗麗， 原 達(dá)， 能昌信

(1. 山東師范大學(xué)信息科學(xué)與工程學(xué)院，山東 濟(jì)南 250014；2. 山東省高校智能信息處理重點(diǎn)實(shí)驗(yàn)室，山東工商學(xué)院，山東 煙臺(tái) 264005)

基于DTW距離的探地雷達(dá)數(shù)據(jù)可視化

段麗麗1， 原 達(dá)2， 能昌信2

(1. 山東師范大學(xué)信息科學(xué)與工程學(xué)院，山東 濟(jì)南 250014；
2. 山東省高校智能信息處理重點(diǎn)實(shí)驗(yàn)室，山東工商學(xué)院，山東 煙臺(tái) 264005)

提出一種基于動(dòng)態(tài)時(shí)間彎曲算法距離度量的探地雷達(dá)數(shù)據(jù)可視化方法，利用動(dòng)態(tài)時(shí)間彎曲算法在時(shí)間軸方向上伸縮的優(yōu)越性，結(jié)合可指定類數(shù)的聚類算法對(duì)探地雷達(dá)數(shù)據(jù)進(jìn)行聚類和可視化分析?？捎糜趯?shí)測(cè)的探地雷達(dá)數(shù)據(jù)集，實(shí)驗(yàn)結(jié)果表明，相對(duì)于傳統(tǒng)的聚類算法，本文算法能得到更好的聚類結(jié)果。

k-means；動(dòng)態(tài)時(shí)間彎曲算法；探地雷達(dá)數(shù)據(jù)；可視化

探地雷達(dá)(ground penetrating radar, GPR)作為一個(gè)非常重要的遙感工具，是近幾年適應(yīng)快速、準(zhǔn)確地?zé)o損探測(cè)地下障礙物以及對(duì)地下工程評(píng)價(jià)而迅速發(fā)展起來的新技術(shù)。在土木工程、地雷探測(cè)和環(huán)境應(yīng)用等領(lǐng)域中，運(yùn)用GPR探測(cè)地下埋藏物或非均勻介質(zhì)一直是被關(guān)注的話題。

與一般的成像方法不同，GPR成像是電磁波經(jīng)過反射、相干和疊加后得到的地下介質(zhì)間接表示。電磁波在地下的傳播過程十分復(fù)雜，各種噪聲和雜波的干擾非常嚴(yán)重，使得探測(cè)回波數(shù)據(jù)伴有多種雜波和相干噪聲。地質(zhì)雷達(dá)噪聲信號(hào)的影響是多方面的，有收發(fā)天線件直接耦合波、空氣直達(dá)波、地表反射波、介質(zhì)環(huán)境反射波、噪聲干擾波和目標(biāo)反射波等，這些噪聲和干擾波有來自環(huán)境的也有來自系統(tǒng)的[1]。在成像過程中，為了反映更多的反射波屬性，通常采用寬頻帶進(jìn)行記錄，于是在圖像中不僅含有各種有效屬性，也存有各種復(fù)雜干擾波，增加了專業(yè)人員的辨識(shí)難度。

目前常用的GPR數(shù)據(jù)處理方法有濾波、反褶積、多次疊加、道均衡處理、速度分析、振幅處理技術(shù)等，其中濾波是最常用的有效手段，目的是利用頻譜特性的不同來壓制干擾波。雷達(dá)反射波，是由所有介質(zhì)的反射響應(yīng)組合而成，可以使用多種成像算法反轉(zhuǎn)，例如時(shí)域標(biāo)準(zhǔn)法相投影(standard back projection,SBP)[2]、傅里葉域雷達(dá)成像技術(shù)[3]等。標(biāo)準(zhǔn)的時(shí)域和頻域成像算法對(duì)數(shù)據(jù)采集過程的脈沖響應(yīng)進(jìn)行匹配濾波，以形成圖像。所有這些算法需要精細(xì)的空間采樣和高速率的時(shí)間采樣[4]。因此，GPR數(shù)據(jù)采集是地下成像過程的關(guān)鍵所在。現(xiàn)在大部分GPR處理軟件的時(shí)頻分析功能強(qiáng)大，但是數(shù)據(jù)解釋能力有限。為了更好地反映GPR數(shù)據(jù)包含的特征和規(guī)律，需要對(duì)處理后的數(shù)據(jù)進(jìn)行更進(jìn)一步的成像處理。

在GPR成像中，使用頻率分集通?？梢詫?shí)現(xiàn)良好的深度分辨率，但橫向分辨率難于實(shí)現(xiàn)。Capineri等[5]提出了一種獲得良好分辨力的 GPR成像方法，它是通過對(duì)B掃描數(shù)據(jù)應(yīng)用Hough變換技術(shù)來實(shí)現(xiàn)的。Morrow和Van Genderen[6]通過運(yùn)用反向傳播和共軛梯度反演技術(shù)形成二維和三維的鉆孔雷達(dá)圖像，但是這些技術(shù)有很大的計(jì)算消耗。因此，Ozdemir等[7]提出了一種新的基于合成孔徑雷達(dá)(synthetic aperture radar，SAR)的GPR成像算法，它能夠在深度和橫向方面同時(shí)達(dá)到良好分辨力。該方法使用快速傅里葉變換形成三維圖像，隨后從距離角度域簡(jiǎn)單變換到圖像域。但是該算法在轉(zhuǎn)換域時(shí)存在圖像失真問題，并且假設(shè)土壤介質(zhì)均勻，土壤屬性為已知，具有一定的局限性。Qin和 Wang[8]提出運(yùn)用基爾霍夫遷移方法來處理GPR數(shù)據(jù)，實(shí)現(xiàn)了對(duì)地下物體的定位和再造。雖然基爾霍夫偏移方法具有較強(qiáng)的能力，能輸出單個(gè)或多個(gè)物體的雷達(dá)輸出圖像，但圖像結(jié)果比較簡(jiǎn)單，并且圖像解釋對(duì)專業(yè)經(jīng)驗(yàn)的要求非常高。

數(shù)據(jù)可視化的思想是運(yùn)用計(jì)算機(jī)圖形學(xué)和圖像處理技術(shù)來表達(dá)數(shù)據(jù)，將離散的、雜亂的數(shù)據(jù)經(jīng)過加工處理，利用計(jì)算機(jī)可視化提供的圖像表示出來，直接觀察數(shù)據(jù)中蘊(yùn)涵的現(xiàn)象或規(guī)律，提高解釋的效率、可靠性和精度[9]，以滿足各行各業(yè)非專業(yè)人員的需求。探地?cái)?shù)據(jù)可視化是探地?cái)?shù)據(jù)解釋的基礎(chǔ)，圖形能夠給用戶一種更直觀的感受，方便用戶對(duì)感興趣的數(shù)據(jù)進(jìn)行操作，幫助解釋人員更好地了解地質(zhì)情況[10]。常用的地質(zhì)剖面顯示方法有3種：波形顯示、變面積顯示和密度顯示[11]。

本文提出一種新的基于動(dòng)態(tài)時(shí)間彎曲算法(dynamic time warping, DTW)距離度量的GPR數(shù)據(jù)可視化方法，對(duì) GPR數(shù)據(jù)進(jìn)行可視化分析，在成像的基礎(chǔ)上產(chǎn)生更直觀有效的視覺效果，為 GPR數(shù)據(jù)的解釋提供可靠的理論基礎(chǔ)?；舅悸肥?，首先對(duì)數(shù)據(jù)進(jìn)行濾波和特征提取，然后運(yùn)用DTW距離與k-means聚類相結(jié)合的方法對(duì)特征數(shù)據(jù)集進(jìn)行聚類分析，最后對(duì)類別信息進(jìn)行顏色映射以得到GPR數(shù)據(jù)的可視化圖像(如圖1所示)。

圖1 可視化流程概念圖

1 探地雷達(dá)波形特點(diǎn)

根據(jù)電磁波傳播和反射理論，接收信號(hào)常常由4個(gè)部分組成：天線串?dāng)_、地面反射、目標(biāo)反射和白噪聲?？晒交癁椋?。由圖2可以看出，GPR回波信號(hào)具有非平穩(wěn)性，脈沖信號(hào)非線性衰減等特點(diǎn)，圖中給出的示例波形從左到右存在明顯的非線性衰減的現(xiàn)象。

圖2 探地雷達(dá)波形圖示例

距離度量用于確定時(shí)間序列之間以及時(shí)間序列分類的相似性。歐氏距離是可以使用的一種有效的距離度量。兩個(gè)相等長(zhǎng)度的時(shí)間序列之間的歐氏距離是簡(jiǎn)單地從一個(gè)時(shí)間序列中的 n個(gè)點(diǎn)到另一個(gè)時(shí)間序列中的 n個(gè)點(diǎn)的距離的總和[12]。它的定義為：

歐氏距離適合某些領(lǐng)域[13]，但其結(jié)果往往不直觀。如圖3所示，如果兩個(gè)時(shí)間序列是相同的，但是一個(gè)沿時(shí)間軸稍稍錯(cuò)開，那么歐式距離會(huì)認(rèn)為是不同的，如圖2中標(biāo)注的①、②兩段波形。DTW 被引入來克服這個(gè)限制，忽視時(shí)間序列的時(shí)間維度的移動(dòng)，并給出直觀的距離測(cè)量。它能夠較好地解決歐氏距離難以處理的對(duì)應(yīng)數(shù)據(jù)間沿時(shí)間軸方向的伸縮、彎曲和線性漂移等問題[14]。與歐式距離相比，DTW 可以度量等長(zhǎng)或不等長(zhǎng)的時(shí)間序列[15]，能夠?qū)崿F(xiàn)異步相似性比較，并且DTW不滿足三角不等式。由于GPR波在傳播過程中會(huì)發(fā)生沿時(shí)間方向的吸收衰減，進(jìn)而對(duì)應(yīng)反射數(shù)據(jù)在時(shí)間軸方向上會(huì)產(chǎn)生壓縮，因此，用DTW距離度量能更好地計(jì)算GPR時(shí)間序列之間的相似性。

圖3 歐氏距離與DTW距離比較

2 動(dòng)態(tài)時(shí)間彎曲

DTW是一種通過彎曲時(shí)間軸來更好地對(duì)時(shí)間序列形態(tài)進(jìn)行匹配映射的相似性度量方法。DTW距離的計(jì)算通過彎曲時(shí)間軸獲取兩個(gè)時(shí)間序列間的最小距離，并確定兩個(gè)時(shí)間序列中各個(gè)點(diǎn)的最佳對(duì)應(yīng)關(guān)系。該算法是要尋找一條通過各個(gè)交叉點(diǎn)的從起始點(diǎn)到終止點(diǎn)的最佳路徑，使得該路徑上所有交叉點(diǎn)的幀失真度總和達(dá)到最小[16]。

圖4 DTW算法求最小彎曲距離

最優(yōu)彎曲路徑是一個(gè)具有最小距離的彎曲路徑，彎曲路徑W的距離是：

式中，Dist(wkm,wkn)是彎曲路徑的第k個(gè)元素中的兩個(gè)數(shù)據(jù)點(diǎn)(一個(gè)來自序列R，一個(gè)來自序列T)索引的距離。

采用基于累積距離矩陣的動(dòng)態(tài)規(guī)劃方法來計(jì)算兩序列之間的最小距離的彎曲路徑。一個(gè)M×N的二維距離代價(jià)矩陣 D被創(chuàng)建，D(m,n)的值是時(shí)間序列R′=＜R1,…,Rm＞和序列T′=＜T1,… ,Tn＞之間的最小彎曲距離。D(M,N)包含了時(shí)間序列R和T之間的最小彎曲距離。圖4展示了從 D(1,1)到D(M,N)的一個(gè)最小距離彎曲路徑實(shí)例。為了找到最小距離的彎曲路徑，代價(jià)矩陣的每個(gè)元素都必須被計(jì)算。構(gòu)造的累積距離矩陣為：

其中，D(0,0)=0，D(m,0)=D(0,n)=+∞。

到 D(m,n)的變形路徑必須通過相鄰的三個(gè)元素之一，由于它們之間的彎曲距離是已知的，所有這一切需要的是加當(dāng)前兩點(diǎn)之間的距離Dist(m,n)到這三個(gè)鄰近元素的最小值上。代價(jià)矩陣從下到上、從左向右填充。因此，整個(gè)矩陣填充后，能找到一條從D(1,1)到D(M,N)的彎曲路徑。

3 基于DTW距離的聚類

本文提出了一種改進(jìn)的基于 DTW 距離的GPR數(shù)據(jù)可視化模型，算法步驟如下：

(1) 對(duì)給定數(shù)據(jù)集Data，進(jìn)行滑動(dòng)窗口均值濾波處理，得到數(shù)據(jù)集D，表示為：

其中，Window表示滑動(dòng)窗口，N表示窗口中像素的總個(gè)數(shù)。

(2) 對(duì)數(shù)據(jù)集 D中的每個(gè)元素提取特征向量Pi，1≤i≤row(D)column(D)，則：

(3) 在特征集P中，任意指定k個(gè)初始聚類中心cj，1≤j≤k，設(shè)定迭代次數(shù)t=0。

(5) 更新每個(gè)簇的聚類中心：

其中，nj是簇 Cj中向量的總個(gè)數(shù)。

(6) 如果：

收斂或達(dá)到指定迭代次數(shù)，輸出類別矩陣；否則，轉(zhuǎn)到步驟(4)。

(7) 對(duì)類別矩陣進(jìn)行顏色映射，最終得到數(shù)據(jù)集的可視化成像。

4 實(shí)驗(yàn)結(jié)果及分析

實(shí)驗(yàn)使用真實(shí)數(shù)據(jù)集進(jìn)行分析。所用數(shù)據(jù)由1 238道數(shù)據(jù)組成，每道數(shù)據(jù)由240字節(jié)道頭和道數(shù)據(jù)組成，采樣點(diǎn)個(gè)數(shù)為2 048(如圖5所示)。

首先，對(duì)原始數(shù)據(jù)進(jìn)行濾波預(yù)處理。通過比較各種尺寸窗口的濾波效果，選擇效果較好的參數(shù)來對(duì)GPR數(shù)據(jù)進(jìn)行處理，本實(shí)驗(yàn)經(jīng)驗(yàn)選取較優(yōu)的大小為 7×7的滑動(dòng)窗口濾波，用幅值彩色成像來對(duì)比濾波前后效果，如圖6所示。

由圖5和圖6的初始成像效果可以觀察到，圖像右下方的數(shù)值范圍比較密集，在灰度圖中區(qū)分度很小，而幅值直接映射圖中，圖像右下方有了較為明顯的類別區(qū)分。

由于GPR波在地下的傳播過程中，會(huì)隨著時(shí)間的推移產(chǎn)生振幅衰減，但頻率等特征是不會(huì)發(fā)生改變的。因此，對(duì)濾波后的數(shù)據(jù)提取多維特征p(如圖7)，得到的樣本總數(shù)為2042×1232，且每個(gè)樣本的特征向量的維數(shù)為5，并用k-means聚類對(duì)數(shù)據(jù)進(jìn)行無監(jiān)督聚類處理。經(jīng)過多次實(shí)驗(yàn)選定 k值為5和7，并同時(shí)對(duì)兩種方法進(jìn)行效果比較，如圖8所示。

從圖8可以看出，在聚類k值相同時(shí)，本文方法比傳統(tǒng)k-means方法具有更好的區(qū)分度，聚類效果更符合預(yù)期的結(jié)果。在k值為5時(shí)，兩種方法對(duì)矩形框中數(shù)據(jù)的聚類結(jié)果均為 2類，但是，可以看出，傳統(tǒng)k-means方法得到的圖像有隨著時(shí)間推移衰減的趨勢(shì)，甚至消失，而本文的方法很好地移除了吸收衰減對(duì)成像效果的影響，成像更接近預(yù)期效果。在 k值為7時(shí)，對(duì)比圖 8(c)和圖 8(d)兩圖，發(fā)現(xiàn)本文提出的方法聚類更細(xì)致，能夠發(fā)現(xiàn)數(shù)據(jù)中更多的潛在特征規(guī)律。而同時(shí)研究圖8(b)和圖8(c)，可以推出，對(duì)數(shù)據(jù)集矩形框區(qū)域的數(shù)據(jù)，傳統(tǒng)k-means聚7類與k-means-DTW聚5類的效果是幾乎一致的。因此，無論是在精細(xì)度還是抗衰減方面，本文提出的k-means-DTW聚類均優(yōu)于傳統(tǒng)的基于歐氏距離的k-means聚類，并且實(shí)驗(yàn)驗(yàn)證了所提方法的有效性。

圖5 原始數(shù)據(jù)成像及單道波形圖

圖6 對(duì)原始數(shù)據(jù)依據(jù)幅值進(jìn)行顏色映射

圖7 特征向量p={橫向梯度，縱向梯度，瞬時(shí)幅值，瞬時(shí)頻率，瞬時(shí)相位}

圖8 兩種方法實(shí)驗(yàn)結(jié)果對(duì)比

5 結(jié) 論

本文提出一種新的基于 DTW距離的動(dòng)態(tài)相似性度量方法，并結(jié)合指定類數(shù)的聚類算法對(duì)GPR數(shù)據(jù)進(jìn)行聚類并可視化分析。處理了數(shù)據(jù)在時(shí)間軸方向上的伸縮和平移問題，一定程度上克服了信號(hào)隨時(shí)間衰減對(duì)可視化結(jié)果的影響。將本文算法應(yīng)用于真實(shí)的GPR數(shù)據(jù)集，并與傳統(tǒng)方法進(jìn)行對(duì)比，結(jié)果表明本文方法可以挖掘到更多的GPR數(shù)據(jù)規(guī)律。在今后的工作中，將結(jié)合多種算法進(jìn)一步研究GPR數(shù)據(jù)的多方面的屬性展示。

[1] 張永山. 探地雷達(dá)異物目標(biāo)探測(cè)與數(shù)據(jù)處理研究[D].太原： 中北大學(xué), 2012.

[2] Feng Xuan,SatoM. Pre-stackMigration applied to GPR for landmine detection [J]. Inverse Problems, 2004, 20(6) 99-115.

[3]Stolt R H.Migration by Fourier transform [J]. Geophysics, 1978,43(1)： 23-48.

[4] Gurbuz A C,McClellan J H,Scott W R. CompressiveSensing for GPR imaging [C]//The41stAsilomar Conference onSignals,Systems and Computers. Pacific Grove, USA, 2007： 2223-2227.

[5] Capineri L, Grande P, Temple J A G. Advanced image‐ processing technique for real‐ time interpretation of ground ‐penetrating radar images [J]. International Journal of ImagingSystems and Technology, 1998, 9(1)： 51-59.

[6]Morrow I L, Van Genderen P. A 2-D polarimetric backpropagation algorithm for ground ‐penetrating radar applications [J].Microwave and Optical Technology Letters, 2001, 28(1)： 1-4.

[7] Ozdemir C, LimS, Hao Ling. ASynthetic‐aperture algorithm for ground‐penetrating radar imaging [J].Microwave and Optical Technology Letters, 2004,42(5)：412-414.

[8] Qin Yao, Wang Qifu. Kirchoff migration algorithm for ground penetrating radar data [C]//ComputerScience and Electronics Engineering (ICCSEE), 2012 International Conference on IEEE. Hangzhou, 2012：396-398.

[9] 張二華, 高 林, 馬仁安, 等. 三維地震數(shù)據(jù)可視化原理及方法[J]. CT理論與應(yīng)用研究, 2007, 16(3)：20-28.

[10] 蔡維緯. 二維地震數(shù)據(jù)可視化的研究與實(shí)現(xiàn)[D]. 昆明：云南大學(xué), 2013.

[11] 肖 漢. 地震數(shù)據(jù)的可視化技術(shù)研究[D]. 長(zhǎng)沙： 湖南大學(xué), 2007.

[12]SalvadorS, Chan P. Toward accurate dynamic time warping in linear time andSpace [J]. Intelligent Data Analysis, 2007, 11(5)： 561-580.

[13] Keogh E, KasettyS. On the need for timeSeries dataMining benchmarks： a survey and empirical demonstration [J]. DataMining and Knowledge Discovery, 2003, 7(4)： 349-371.

[14] 程文聰, 鄒 鵬, 賈 焰, 等. 基于DTW距離的偽周期數(shù)據(jù)流異常檢測(cè)[J]. 計(jì)算機(jī)研究與發(fā)展, 2010,47(5)： 893-902.

[15] 陸薛妹, 胡 軼, 方建安. 基于分段極值DTW距離的時(shí)間序列相似性度量[J]. 微計(jì)算機(jī)信息, 2007, 23(9-3)：204-206.

[16] 王 倩, 吳國平, 陳 琳. 特定人語音識(shí)別算法——DTW算法[J]. 軟件導(dǎo)刊, 2005, (20)：46-48.

Ground Penetrating Radar Data Visualization Based on Dynamic Time Warping

Duan Lili1, Yuan Da2, Nai Changxin2
(1.School of InformationScience and Engineering,Shandong Normal University, JinanShandong 250014, China; 2. Key Laboratory of Intelligent Information Processing in Universities ofShandong,Shandong Institute of Business and Technology, YantaiShandong 264005, China)

In this paper, we propose a visualizationMethod of ground penetrating radar data based on dynamic time warping distanceMetric, which uses the advantages of dynamic time warping algorithm on the time axisScaling, and combines with theSpecified clustering algorithm for ground penetrating radar data clustering and visualization analysis. The algorithm is used for theMeasured ground penetrating radar dataSets; the experimental resultsShow that compared with the traditional clustering algorithm, the proposed algorithm can get better clustering results.

k-means; dynamic time warping; ground penetrating radar data; visualization

TP 391

A

2095-302X(2015)02-0152-07

2014-10-08；定稿日期：2014-10-30

國家自然科學(xué)基金資助項(xiàng)目(61373079)

段麗麗(1989–)，女，山東濰坊人，碩士研究生。主要研究方向?yàn)榭梢暬c圖像處理。E-mail：duanlili.123@qq.com通訊作者：原 達(dá)(1968–)，男，遼寧建昌人，教授，博士。主要研究方向?yàn)榭梢暬c圖像處理。E-Mail：ydccec@126. com