地質雷達數據中克里金插值采樣數據選擇算法

王振武， 卜異亞， 馬鍵強

(中國礦業大學(北京) 機電與信息工程學院，北京 100083)

?

地質雷達數據中克里金插值采樣數據選擇算法

王振武， 卜異亞， 馬鍵強

(中國礦業大學(北京) 機電與信息工程學院，北京 100083)

采樣數據選擇的合理性是克里金算法插值結果是否準確的前提。針對如何從地質雷達數據中選擇恰當的克里金插值采樣數據的問題，本文提出了兩種克里金插值采樣數據選擇算法，即動態三維球體覆蓋(D3DGC)選點算法和雙次反距離選點(DRD)算法。兩種算法在基于粒子群優化的普通克里金算法下驗證了采樣數據選擇的合理性，實驗結果分析表明，D3DGC算法在最佳選點個數下誤差率比傳統選點方法降低了6.7%，在其他選點個數下誤差率(<10%)也明顯低于傳統選點方法，DRD算法在參數e和c為4以及最佳選點個數的條件下誤差率也達到了6.2%，在其他情況下誤差率也低于傳統選點方法。

克里金算法；插值；地質雷達；采樣數據選擇算法；動態三維球體覆蓋選點算法；雙次反距離選點算法；粒子群優化；誤差率

地質雷達(ground penetrating radar，GPR)技術是一種基于電磁波反射原理，用于淺層地質構造探測和工程治理檢測的地球物理方法[1]。為了深入勘查或檢測探測區域的地質構造或地下結構，需要對測線數據做插值處理來獲取探測區域的整體數據。在常見的空間插值方法中以克里金(Kriging)插值方法為代表的地質統計方法提供了揭示空間變量結構相關性與隨機性的理論，是地質領域應用最廣泛的插值方法之一[2]。在對克里金算法的研究中，有些學者根據不同的應用情況，將不同類型的克里金算法組合使用。牛文杰等將同位置協同克里金(Collocated CoKriging)和貝葉斯克里金方法相結合，推導出一種新的貝葉斯同位置協同克里金估值方法，綜合了上述兩種方法的優點[3]。李曉軍等針對不連續地層，提出了一種采用指示克里金和普通克里金相結合的方法，即先用指示克里金法估計地層分布范圍，再用普通克里金方法估計分布范圍內的地層厚度，更多的學者將克里金方法和其他方法結合使用[4]。王輝贊等針對變異函數難以刻畫實際數據分布的不足，采用最小二乘支持向量機從實際資料場擬合重構變異函數，表現出較好的針對性和客觀性[5]。賈雨等將粒子群算法和克里金算法相結合，在粒子群優化過程中，通過高斯變異、樣本點權重系數設定、搜索范圍約束等方式提高了插值精度[6]。顏華等針對三維溫度場重建問題，首先采用最小二乘法獲得一個由少量像素描述的溫度場，然后運用克里金內插和外推運算獲得整個被測溫度場的細致描述，這種方法可實現復雜的三維溫度場的高精度重建[7]。徐馳等將高光譜反演與協同克里金方法相結合，將土壤表層含水率反演結果作為協同克里金插值的協變量來輔助進行協同克里金插值，提高了預測效率[8]。ZHANG等先采用克里金算法進行大型球磨機的功能適配，然后利用遺傳算法優化了齒輪傳送的可靠性問題[9]。FANG等將序列克里金算法和人工蜂群算法結合來優化點焊接頭以改善其疲勞壽命，該方法用克里金模型來近似表達點焊接頭設計模型的目標函數，然后采用人工蜂群算法進行迭代優化[10]。

在上述的研究中，對克里金插值方法的研究多集中在對不同類型克里金算法的組合[3-4]以及與其他算法理論的結合[5-11]中，鮮有學者研究克里金插值方法中采樣數據的選擇問題。一般情況下，按照空間相關性原理，地質數據應該選取距離插值點較近的數據點作為采樣數據。文獻[12]針對某儲層數據，研究了空間相關分析的主要影響因素與克里金估計結構的關系，指出影響插值結構的因素和參數包括采樣點的位置分布特征、變差函數模型的塊金效應、變差函數模型的變程、變差函數模型的類型、空間結構的各向異性以及空間搜索范圍等。文獻[11]進一步指出，距離估計點較近的采樣數據對估計結構的影響較大，且大多數情況下應該選用球狀模型或指數模型。但文獻[11]也只是給出了一般性的選點原則，對具體數據而言，應該如何選擇采樣數據并沒有進行深入分析。文獻[2]針對傳統的克里金方法采用歐式距離測度描述空間結構相關性的缺點，基于最小方差準則，提出了綜合曼氏、歐式和切氏三種距離測度的多測度加權克里金方法，并通過實驗表明該方法提高了插值精度和可靠性，但是該方法也僅僅是對距離測度公式做了擴展，沒有深入討論采樣數據如何選擇的問題，而選擇合適的采樣數據是提高克里金插值準確性的前提。論文針對這一問題提出了動態三維球體覆蓋(dynamic three-dimensional globe cover，D3DGC)選點算法和雙次反距離(double reverse distance，DRD)選點算法，并基于粒子群優化的普通克里金算法驗證了兩種算法的有效性。

1 克里金插值采樣數據選擇算法描述

1.1 基本概念

本文在已獲取地質雷達測線數據的基礎上研究克里金插值的采樣數據選擇問題。為了對所提出的算法進行描述，首先給出相關概念。

定義1 地質數據點δ=〈ζ，γ〉：δ是一個二元組，其中ζ=〈x，y，z〉表示三維坐標值，γ表示地質點的屬性值，如物質對電磁波的反射值等。

(1)

式中：δij表示測線數據上的第i行j列的地質數據點，且1≤i≤m，1≤j≤n。為表述方便，設一共p條測線數據，第l條測線表示為l，將測線上第h行的地質數據點集合定義為δh，h∈[1，m]，l∈[1，p]。

定義3 測線數據集Γ：Γ表示對同一探測區域的同一次探測活動中獲取的測線數據總和，即

Γ={1，2，…，lp}

(2)

定義4 測線數據凸包Δ：Δ表示當前測線數據集所描述的最大三維探測區域。

Δ={δ|δ.ζ.x∈[Xmin，Xmax]，δ.ζ.y∈[Ymin，Ymax]，
δ.ζ.z∈[Zmin，Zmax]}

Xmin=min{k.δij.ζ.x}，Xmax=max{k.δij.ζ.x}
Ymin=min{k.δij.ζ.y}，Ymax=max{k.δij.ζ.y}
Zmin=min{k.δij.ζ.z}，Zmax=max{k.δij.ζ.z}

(3)

式中：k∈[1，p]，i∈[1，mk]，j∈[1，nk]，Xmin和Xmax分別表示x坐標的最大和最小值，min{ }和max{ }表示取集合的最小值和最大值操作。

定義5 待插值數據點集ψ：ψ=Δ-Γ，克里金插值的目的就是從Γ中選擇合適的采樣點對ψ進行插值。

定義6 凸包水平分割層η：η表示Δ按照測線平行方向的層次劃分，即

(4)

Ψ={φ1，φ2，…,φu}

(5)

式中：φr表示第r層的待插值數據點的集合， 1≤r≤u且存在關系測線l中的δi滿足δi?l.δi，l∈[1，p]，i∈[1，m]。

定義7 待插值點的采樣數據點集φr[ε]→π[μ]：假設選擇μ個采樣點進行克里金插值，則φr中第ε個插值點的μ個采樣點的集合表示為φr[ε]→π[μ]。

定義8 距離計算公式τ：論文采用歐幾里得距離作為距離測度公式。設ψ中的某待插值點為θ(a，b，c)，測線上某點l中δ的ζ值為(x，y，z)，則距離τ定義為

(6)

1.2 采樣數據選擇算法原理

地質數據具有空間相關性和隨機性的特點，針對不同數據的特點，本文提出了D3DGC選點算法和DRD選點算法。

1.2.1 D3DGC 算法

根據地質數據的相關性原理，一般而言選取距離插值點較近的數據點作為采樣點，針對地質雷達數據特點，論文提出了D3DGC選點算法，算法原理如圖1所示。

圖1 D3DGC算法原理圖Fig.1 Theory graph of D3DGC

定義9 球與直線交點方程：設待插值點為θ(a，b，c)， 球的半徑為R，l上的δh與球的交點φ(x，y，z)計算公式為

(7)

設測線垂方向采樣點步長為d，每個待插值點選取μ個采樣點進行插值?；谌缟隙x，D3DGC算法的流程如圖2所示。

圖2 D3DGC算法流程圖Fig.2 The flow graph of D3DGC

對于每一個待插值點θ(a，b，c)，球的半徑以d為單位遞增，計算l的δh與球體交截后包含在球體內的地質數據點個數kl，這里直到條件滿足時為止。通過逐層計算ηi上的φi，直到所有待插值點的采樣點選擇完畢為止。

如果θ(a，b，c)的插值采樣點只從ηi中的δi選擇，其中δi屬于l，1≤≤p則D3DGC算法就簡化為二維模式，如圖3所示，稱之為動態圓覆蓋(dynamic circle cover，DCC)選點算法。

準時制（JIT)，可概括為“任何時間、任何地點、任何事情”，都可以及時解決和處理，以達到零時間、零距離的高效作業模式。審批平臺可以引入此理論在審批業務流程的設計和模型的架構當中，以實現審批的高效、便捷。具體JIT模式的優勢見圖1。

1.2.2 DRD算法

圖3 DCC算法原理圖Fig.3 Theory graph of DCC algorithm

定義10 點到直線的距離?。設φi中當前的待插值點為θ(a，b，c)，l上的δh可由式(7)的方程2描述，若將方程2簡記為則θ(a，b，c)到l上的δh的距離可定義為

(8)

設待插值點選取μ個樣本點進行插值，.δh的點數為n(為插值準確性和效率考慮，通常μ?n)，按照式(8)計算距離后，分配到第l條測線上的采樣點數由式(9)計算：

(9)

(10)

圖4 DRD算法原理圖Fig.4 Theory graph of DRD

圖5 DRD算法流程圖Fig.5 The flow graph of DRD

2 算法實驗分析

由于主觀設置的變差函數模型參數對克里金算法的插值結果影響較大，為了客觀準確地驗證兩種采樣點選擇算法，論文采用了粒子群優化(particle swarm optimization，PSO)算法對普通克里金算法的變差函數模型參數進行優化。

PSO算法將優化問題的候選解假設為一個無體積和質量的粒子，在d維空間中粒子通過自身和其他粒子的經驗來調整自己的位置。設有n個粒子組成的粒子群，其中第i個粒子在d維空間中的速度為vi=[vi1vi2…vid]T，位置為xi=[xi1xi2…xid]T，粒子i的當前個體最優位置為pid，整個粒子群的當前全局最優位置為gd，粒子i通過式(11)的迭代來更新自己的速度和位置。

(11)

定義11 適應度函數F:本文給定的PSO算法的適應度函數表示為

(12)

式中：F(δo) 為第o個粒子(地質數據點)的適應度函數值，V*(δo)表示粒子的克里金估計量，V(δo)表示例子的實際值。論文采用普通Kriging算法進行實驗結果的驗證，變差函數采用球狀模型，拱高f、塊金值g、變程s，參數采用PSO算法進行優化，粒子的速度和位置用向量〈f,g,s〉表示，粒子群的最大迭代次數α=1 000，適應度函數閾值β=10-6。PSO優化的克里金(PSO-Kriging)算法的流程如圖6所示。

圖6 PSO-Kriging算法流程圖Fig.6 The flow graph of PSO-Kriging

本文基于地質雷達探測數據來驗證所提出的兩種采樣點選擇算法的有效性。地質雷達儀為中國礦業大學(北京)自主研發的MTGR-4F車載地質雷達，雷達參數為：探測深度5m，時間窗口100ns，采樣點數512，天線主頻200MHz。實驗數據探測路段為北京市朝陽區東三環中路，探測目的有3個：地底的分層結構、有無地下管線以及有無空洞。實驗數據取值范圍為[-23 679，15 421]；測線條數為3，采樣道數分別為18 632、18 088和18 838，實驗所用數據的采集工作照如圖7所示。

圖7 實驗數據的獲取Fig.7 The acquirement of experiment data

實驗硬件和軟件環境信息如下：CPU為Pentium(R)Dual-CoreT4300 2.10GHz，內存為RAM2GB，操作系統為Windows7 32位，實驗平臺為VS.NET2012，開發語言為C++。

傳統的采樣點選擇方式是對同一ηi，1≤i≤u上所有l的δi數據按照式(6)進行排序，選擇距離最近的某些地質數據點作為克里金插值的采樣點。本文將所提出的兩種算法和這種傳統的選點方法進行了比較，為準確起見，誤差率統計的為N0=2 000個點的平均誤差率，誤差率ρ由式(13)計算：

(13)

實驗探測數據的最優采樣數據點數與采樣道數、采樣點數等具體的地質雷達儀參數有關?？紤]到插值算法的運算效率和插值結果精度等因素，本實驗選取的數據點個數μ∈[160，300]，依次按10個點的數量遞增。如圖8所示，在不同的選點個數下，D3DGC算法的誤差率明顯低于傳統選點算法。當選點個數為210時，兩種算法的誤差率都為最小，D3DGC算法的誤差率僅為1.3%，而傳統選點算法的誤差率卻為8%，而且當選點個數在180～270時，D3DGC算法的誤差率都小于10%，而傳統選點算法只有在200～210誤差才小于10%，且都大于8%，在其他范圍內的最大誤差達到了31.1%。與傳統選點算法相比，D3DGC算法效果好的原因在于，該算法在三維空間中考慮相關性，按照球狀半徑遞增的原理逐步尋找采樣數據點，比傳統選點算法更為準確合理。

圖8 D3DGC算法實驗結果Fig.8 The experiment results of D3DGC

表1 DRD算法實驗結果

3 結論

1)D3DGC算法是在三維空間中按照球狀半徑遞增原理選取的采樣點，因此其克里金插值誤差率明顯地比傳統選點方法好；2)針對數據量大且相鄰數據值較為接近的地質雷達數據，DRD算法綜合考慮了地質數據的空間相關性和隨機性，在參數e和c取值為4的條件下也取得了不錯的效果。

D3DGC和DRD算法均已在道路塌陷的檢測問題中得到了應用，使得克里金插值的精度得到了提高。對于DRD算法，針對不同的地質雷達數據，如何實現參數e和c的自適應取值是將來需要進一步研究的問題。

[1]李晉平，邵丕彥，谷牧.地質雷達在鐵路隧道工程質量檢測中的應用[J].中國鐵道科學， 2006， 27(2)： 56-59.

LI Jinping， SHAO Piyan， GU Mu. Application and analysis of GPR in railway tunnel engineering quality inspection[J]. China railway science， 2006，27(2)： 56-59.

[2]劉志平，何秀鳳，張淑輝.多測度加權克里金法在高邊坡變形穩定性分析中的應用[J].水利學報， 2009， 40(6)： 709-715.

LIU Zhiping， HE Xiufeng， ZHANG Shuhui. Multi-distance measures weighted Kriging method for deformation stability analysis of steep slopes[J]. Shuili xuebao， 2009，40(6)： 709-715.

[3]牛文杰，朱大培，陳其明.貝葉斯同位置協同克里金(Bayesian Collocated CoKriging)估值法的研究[J]計算機輔助設計與圖形學學報， 2002， 14(4)： 343-347.

NIU Wenjie， ZHU Dapei， CHEN Qiming. Bayesian collocated CoKriging[J]. Journal of computer-aided design & computer graphics， 2002， 14(4)： 343-347.

[4]李曉軍，張振遠.基于指示和普通克里金的不連續地層厚度估計方法[J].巖土力學， 2014， 35(10)： 2881-2887.

LI Xiaojun， ZHANG Zhenyuan. A combined indicator-ordinary Kriging method for estimating thicknesses of discontinuous geological strata[J]. Rock and soil mechanics， 2014， 35(10)： 2881-2887.

[5]王輝贊，張韌，劉巍，等.支持向量機優化的克里金插值算法及其海洋資料對比試驗[J].大氣科學學報， 2011， 34(5)： 567-573.

WANG Huizan， ZHANG Ren， LIU Wei， et al. Kriging interpolation method optimized by support vector machine and its application in oceanic data[J]. Transactions of atmospheric sciences， 2011， 34(5)： 567-573.

[6]賈雨，鄧世武，姚興苗，等.基于約束粒子群優化的克里金插值算法[J].成都理工大學學報：自然科學版， 2015， 42(1)： 104-109.

JIA Yu， DENG Shiwu， YAO Xingmiao， et al. Kriging interpolation algorithm based on constraint particle swarm optimization[J]. Journal of Chengdu University of Technology：science & technology edition， 2015， 42(1)： 104-109.

[7]顏華，李欣，王善輝.基于最小二乘法和克里金插值的三維溫度場重建[J].沈陽工業大學學報， 2014， 36(3)： 303-307.

YAN Hua， LI Xin， WANG Shanhui. Reconstruction of three-dimensional temperature field based on least-square method and kriging interpolation[J]. Journal of Shenyang University of Technology， 2014， 36(3)： 303-307.

[8]徐馳，曾文治，黃介生，等.基于高光譜與系統克里金的土壤耕作層含水率反演[J].農業工程學報， 2014， 30(13)： 94-103.

XU Chi， ZENG Wenzhi， HUANG Jiesheng， et al. Modelling soil water content in plow layer using co-kriging interpolation based on hyperspectral data [J]. Transactions of the Chinese society of agricultural engineering， 2014，30(13)： 94-103.

[9]ZHANG Guanyu， WANG Guoqiang， LI Xuefei， et al. Global optimization of reliability design for large ball mill gear transmission based on the Kriging model and genetic algorithm[J]. Mechanism and machine theory， 2013， 69:321-336.

[10]FANG Jianguang， GAO Yunkai， SUN Guangyong， et al. Optimization of spot-welded joints combined artificial bee Colony algorithm with sequential sriging optimization[J]. Advances in mechanical engineering， 2014: 1-10.

[11]DENG Mingguo， LI Wenchang， LI Bo， et al. Application of log kriging on estimated reserves of the 10-9 ore body of lutangba in the gejiu tin deposits[J]. Journal of China university of mining & technology， 2007， 17(2)： 286-289.

[12]劉永社，印興耀，賀維勝.空間相關分析因素對儲層建模中克里金估計結果的影響[J].石油大學學報：自然科學版， 2004， 28(2)： 24-29.

LIU Yongshe， YIN Xingyao， HE Weisheng. Effects of the factors for spatial correlation analysis on Kriging estimation in reservoir modeling[J]. Journal of the University of Petroleum： science & technology edition， 2004， 28(2)： 24-29.

本文引用格式：

王振武， 卜異亞， 馬鍵強. 地質雷達數據中克里金插值采樣數據選擇算法[J]. 哈爾濱工程大學學報， 2017， 38(5)： 784-790.

WANG Zhenwu，BU Yiya，MA Jianqiang. Sampling data selection algorithm for Kriging interpolation based on ground penetrating radar data[J]. Journal of Harbin Engineering University， 2017， 38(5)： 784-790.

Sampling data selection algorithm for Kriging interpolation based on ground penetrating radar data

WANG Zhenwu，BU Yiya，MA Jianqiang

(School of Mechanical Electronic & Information Engineering， China University of Mining & Technology (Beijing)， Beijing 100083， China)

The rationality of the data sampling selection method is a prerequisite for judging whether the Kriging interpolation results are accurate. Considering the problem of how to properly choose sampling data from the probing data of ground penetrating radar， in this paper， we propose two sampling selection algorithms for Kriging interpolation: the dynamic three-dimensional globe cover (D3DGC) algorithm and the double reverse distance (DRD) algorithm. We verified the rationality of these two methods with respect to the ordinary Kriging algorithm based on particle swarm optimization. The experimental results show the D3DGC algorithm to reduce the error rate to 6.7% compared with the traditional method for an optimal number of sampling points. With respect to other sampling numbers， D3DGC also demonstrated obvious advantages (error rate<10%). The DRD algorithm error rate was 6.2% for the optimal number of sampling points when the parameter e and c values were 4， respectively. For other sampling numbers， the DRD algorithm also obtained better results than the traditional method.

Kriging algorithm; interpolation; ground penetrating radar; sampling selection algorithm; dynamic three-dimensional globe cover; double reverse distance; particle swarm optimization(PSO); error

2016-01-07.

日期：2017-04-26.

國家高技術研究發展計劃重大專項(2012AA12A308)；核設施退役及放射性廢物治理科研項目(FZ1402-08)；國家自然科學基金項目(61302157)；北京市高等學校青年英才計劃項目(YETP0939).

王振武(1978-)， 男， 副教授，博士.

王振武，E-mail: wangzhenwu@126.com.

10.11990/jheu.201601022

TP39

A

1006-7043(2017)05-0784-07

網絡出版地址：http://www.cnki.net/kcms/detail/23.1390.u.20170426.1152.046.html