朱家成 田善君

（1.江漢大學(xué)數(shù)學(xué)與計算機科學(xué)學(xué)院計算中心 武漢 430056）（2.中國地質(zhì)大學(xué)（武漢）計算機學(xué)院地質(zhì)信息技術(shù)研究所 武漢 430074）（3.山東交通學(xué)院交通土建工程學(xué)院 濟南 250000）

1 引言

變差函數(shù)是Motheron在1965年提出的一種區(qū)域變化估計方法。數(shù)學(xué)解釋為區(qū)域化變量增量平方的數(shù)學(xué)期望，變差函數(shù)的實際意義是，它反映了區(qū)域化變量在某個方向上某一距離范圍內(nèi)的變化程度，是地質(zhì)統(tǒng)計學(xué)克里格方法的基礎(chǔ)操作工具，能夠有效地描述區(qū)域化變量在空間上的結(jié)構(gòu)性與隨機性［1～3］，并可以在其他領(lǐng)域作為分析區(qū)域化變量隨機性和結(jié)構(gòu)性特征的有效工具。

在地質(zhì)統(tǒng)計學(xué)領(lǐng)域變差函數(shù)的實際應(yīng)用中，我們通常都將求解變差函數(shù)過程分解為實驗變差函數(shù)計算和理論變差函數(shù)擬合，目前主流的研究都是針對理論變差函數(shù)的擬合，因為理論擬合的參數(shù)好壞能決定一個變差函數(shù)應(yīng)用于地質(zhì)統(tǒng)計學(xué)的精確度。所以早期的地質(zhì)統(tǒng)計學(xué)研究領(lǐng)域主要是對變差函數(shù)理論擬合進行優(yōu)化，很少有涉及對實驗變差函數(shù)的計算過程進行優(yōu)化。而得到一個穩(wěn)健的實驗變差函數(shù)是進行理論擬合的基礎(chǔ)，且在當(dāng)前實際應(yīng)用中，局限變差函數(shù)應(yīng)用的一個主要問題就是數(shù)據(jù)量增大時，急劇下降的實驗變差函數(shù)的計算速度。真正在地質(zhì)統(tǒng)計學(xué)的應(yīng)用中，得到一個能滿足實際生產(chǎn)應(yīng)用的實驗變差函數(shù)結(jié)果，需要不斷進行不同方向上多個參數(shù)下的變差函數(shù)計算調(diào)試，需要來回多次進行實驗結(jié)果的比對，而每次計算的過程都是一次復(fù)雜的運算［1］。

在地質(zhì)統(tǒng)計學(xué)儲量估算領(lǐng)域，我們會對三維空間內(nèi)的鉆孔采樣點建模之后進行變差函數(shù)的計算，從而可以更加形象具體且真實的反映某個礦體的品位含量的三維空間的結(jié)構(gòu)性和隨機性，進而估算儲量。而真實的地礦數(shù)據(jù)較為龐大，在實際應(yīng)用過程中，阻礙儲量估算的難題已經(jīng)不再是對理論變差函數(shù)擬合的優(yōu)化來提高估算精確度，而是繁瑣的實驗變差函數(shù)計算嚴(yán)重影響了儲量估算的計算時間。一次完整的三維儲量估算過程，采用傳統(tǒng)的方法，需要多次的實驗變差函數(shù)調(diào)試，每次調(diào)試的過程都需要大量的計算，嚴(yán)重影響了整個計算過程，所以地質(zhì)統(tǒng)計學(xué)領(lǐng)域急需一個優(yōu)化的實驗變差函數(shù)計算算法來解決這個難題。

2 三維空間實驗變差函數(shù)

變差函數(shù)的基本定義是區(qū)域化變量增量平方的數(shù)學(xué)期望，也就是區(qū)域化變量增量的方差［3］。通常為了便于生產(chǎn)計算需要，一般取值半變差函數(shù)，本文在計算過程中，也遵循這一規(guī)則。總得來說可表述為在相距為h的兩個空間點x和x+h的權(quán)值Z（x）和Z（x+h）之間的方差，因此可得計算公式如下：

式（1）是在理想狀態(tài)下，有無窮的樣品對參與計算，最后求出一個平均值，其中，Z（x+h），Z（x）表示在x+h與x兩個不同位置時，樣品點的權(quán)值，而上式又稱為理論變差函數(shù)。而在實際應(yīng)用中，區(qū)域內(nèi)變量的數(shù)量是有限的，且取不同的增量值，會對樣品對數(shù)影響較大，增量值越大，樣品數(shù)越少；增量值變小，樣品數(shù)目會增多。而樣品對數(shù)越多，算出的結(jié)果才更準(zhǔn)確，有實際意義，但樣品數(shù)越多則會造成計算量的增加。而在地質(zhì)統(tǒng)計的實際應(yīng)用中，確定的空間位置點x0上只能有Z（x0）的一次觀測值，即樣品對的個數(shù)是有限的，故可以得到式（2）。

在式（2）中，N（h）表示相距為h的樣品對的數(shù)目，顯然可以看出，實驗變差函數(shù)是一個根據(jù)h的變化而變化的一條函數(shù)曲線，如圖1所示。

圖1 實驗變差函數(shù)曲線圖

2.1 不規(guī)則樣品對的選取

在真實的環(huán)境中，樣品數(shù)據(jù)往往都是不規(guī)則分布的，這就需要我們采用一定的處理方法才能獲取能夠參與運算的有效數(shù)據(jù)。

2.1.1 非三維空間的不規(guī)則數(shù)據(jù)處理

1）一維數(shù)據(jù)

當(dāng)數(shù)據(jù)處于一維分布，也就是待處理的數(shù)據(jù)都在一條直線上。如果我們還是按照固定的步長來計算的話，能夠滿足要求的樣品對會非常少，所以我們增加在距離上的約束，也就是兩個樣品點之間的距離在步長的約束范圍之內(nèi)存在一定的容差，滿足容差的樣品對都屬于合格樣品對。

距離容差可以保證在一個確定的方向上增加入選的樣品對，但容差也不能設(shè)置的太大，不然就會導(dǎo)致很多不符合計算要求的樣品對也被規(guī)劃進來。

2）二維數(shù)據(jù)

當(dāng)數(shù)據(jù)處于二維分布，也就是待處理的數(shù)據(jù)都在一個平面上。那么所選樣品對除了在距離上進行容差設(shè)計之外，還要在方向上進行一定的容差設(shè)計，也就是在方向上也要增加一定的角度范圍活動區(qū)間，屬于該范圍內(nèi)的樣品對都屬于符合要求的樣品。如圖2所示。

圖2 角度容差示意圖

假設(shè)點A分別與點BCD來進行角度判斷，有箭頭標(biāo)志的直線為基準(zhǔn)方向，那么樣品對（AB）、（AC）、（AD）都屬于符合基準(zhǔn)方向的樣品（假定夾角小于等于角度容差）。

2.1.2 三維空間不規(guī)則數(shù)據(jù)的處理

顯然，真實的地質(zhì)統(tǒng)計學(xué)過程中，數(shù)據(jù)都是在三維空間分布，在一維和二維的數(shù)據(jù)分析基礎(chǔ)上，把距離和角度的約束都添加進去，為了便于進行更好的數(shù)據(jù)結(jié)果分析比對與計算，以及算出的結(jié)果更有實際意義，可以參照自己一維，二維分布的容差涉及，在三維實際應(yīng)用中也對每次參與的樣品進行角度和距離的容差設(shè)計，便于更好地獲取能夠參與計算的樣品對，如圖3所示。

圖3 三維角度距離容差示意圖

2.2 三維穩(wěn)健實驗變差函數(shù)結(jié)果的調(diào)試

為了更好地獲得區(qū)域化變量在空間的結(jié)構(gòu)性和隨機性，以及為下一步的克里格插值做準(zhǔn)備，在進行三維空間變差函數(shù)分析時，通常會進行三個兩兩相互垂直的方向作為求得實驗變差函數(shù)結(jié)果的方向。這個三個方向的獲取，則是根據(jù)式（2）進行不同方向變差函數(shù)結(jié)果的比對而來。根據(jù)以上分析我們可以得到三維穩(wěn)健實驗變差函數(shù)結(jié)果的計算流程大致如下：

簡單來說計算流程主要可分為以下幾大步驟：

1）樣品數(shù)據(jù)模型的建立，即區(qū)域化變量的位置信息以及品位信息建模，即公式中任意樣品點A，B的位置xi與權(quán)值Z（xi）；

2）約束范圍模型的建立，確定角度約束，距離約束，以及其他條件約束，得到約束范圍內(nèi)的兩兩樣品對庫；

圖4 三維實驗變差函數(shù)流程圖

3）計算符合約束范圍的單個實驗變差函數(shù)值，得到當(dāng)前范圍的實驗變差函數(shù)結(jié)果；

4）變差函數(shù)結(jié)果的調(diào)試，專家進行當(dāng)前實驗變差函數(shù)結(jié)果的判定是否能達到后續(xù)工作的要求，進行不同的調(diào)試與重新計算的選擇；

5）穩(wěn)健實驗變差函數(shù)結(jié)果的生成，進行多次的實驗變差函數(shù)結(jié)果的調(diào)試之后，可以先后確定三個方向上的穩(wěn)健三維空間實驗變差函數(shù)結(jié)果。如圖所示黑色和灰色實線表示相互垂直的第一第二方向，虛線表示與前兩個方向分別垂直的第三方向。

圖5 三個相互垂直方向計算樣品圖

3 GPU加速的優(yōu)化算法

由以上步驟可以看出實驗變差函數(shù)的計算中，每一個Rz的計算是在確定一個固定的參數(shù)之后繼而實施同種類型的運算，具有模塊化和數(shù)據(jù)弱相關(guān)的特點，符合CUDA編程模型的“數(shù)據(jù)并行處理”機制，故可以采用GPU進行計算性能上的加速。該步驟主要體現(xiàn)在計算流程中的“計算半變差值”模塊里。

3.1 計算過程的優(yōu)化整合

基于GPU架構(gòu)的CUDA編程模型本質(zhì)上是一種流式計算，適用于數(shù)據(jù)并行處理的問題。數(shù)據(jù)并行是將數(shù)據(jù)劃分成子數(shù)據(jù)流，在一些相同執(zhí)行過程的程序模塊上對這些子數(shù)據(jù)流進行處理。這種方式的并行效果并不受問題求解步數(shù)的限制，只受制于相同處理模塊的數(shù)目和系統(tǒng)內(nèi)部的通信寬帶，并行度很高且具有良好的可擴展性，可以構(gòu)建包含大規(guī)模運算模塊的空間分析系統(tǒng)。

傳統(tǒng)的實驗變差函數(shù)計算效率低下的原因在于：傳統(tǒng)的算法在實現(xiàn)時沒有考慮到整個實驗變差函數(shù)值得計算過程所表現(xiàn)的內(nèi)在的處理并行性。

通過對式（1）變差值的計算過程，我們可以看出，每次重新約束規(guī)則進行一次新的計算，都是進行若干個半變差值的計算，而一個個樣品對進行變差值的計算又是相互獨立互不影響的過程，即多次不同參數(shù)下的同種類型的運算，故可在上述傳統(tǒng)三維穩(wěn)健實驗變差函數(shù)的計算流程中，將計算半變差值的過程，用CUDA編程模型進行運算過程的加速。

圖6 運算過程的替換

3.2 算法實現(xiàn)

3.2.1 建立模型

建立數(shù)據(jù)模型即創(chuàng)建一個需要進行實驗變差函數(shù)計算的三維模型，包括幾何信息模型與屬性信息模型。為了更便于操作以及空間展示，我們將每一個樣品都設(shè)計成三維空間圓點或者中心線的形式，屬性信息以隱藏的方式存儲，只有在調(diào)用的時候才展示出來。

Input：樣品的幾何信息，屬性信息。

Output：樣品數(shù)據(jù)模型。

1）創(chuàng)建樣品對象CPoint（X，Y，Z，Atr）；

//其中X，Y，Z表示樣品對象的空間坐標(biāo)，Atr表示該樣品對象的屬性權(quán)值；

2）依次循環(huán)，將幾何信息，屬性信息讀入創(chuàng)建的每一個對象；

3）進行數(shù)據(jù)模型信息的可視化輸出以及結(jié)構(gòu)輸出；

//可視化輸出表示對象的三維展示，結(jié)構(gòu)輸出表示記錄當(dāng)前的數(shù)據(jù)文本信息。

而建立約束模型，則需要在根據(jù)不同的專業(yè)領(lǐng)域，根據(jù)專家意見設(shè)定一個初步方向約束與距離約束，我們在一定范圍內(nèi)，設(shè)定方向約束與距離約束的調(diào)試范圍。

Input：約束規(guī)則的參數(shù)值，約束的表現(xiàn)形式。

Output：樣品的約束模型。

1）創(chuàng)建約束的表現(xiàn)形式CSearchCone，Cline；

//CSearchCone，Cline分別是一個圓錐體和直線段，為方向約束與距離約束的表現(xiàn)形式

2）根據(jù)樣品數(shù)據(jù)模型的參數(shù)確定約束規(guī)則的大致參數(shù)范圍；

3）將獲取的約束規(guī)則參數(shù)值賦予以上兩個約束對象；

4）確定約束變動的規(guī)則與范圍；

5）約束模型的可視化輸出與結(jié)構(gòu)輸出；

同樣下圖展示了約束模型在三維圖形顯示平臺中的展示形式與形態(tài)，圓錐體表示方向約束，中間的直線距離約束。

在進行約束規(guī)則更改之后，約束方向，或者約束長度都可以進行改變。

3.2.2 并行計算過程，建立樣品整體庫

并行計算的過程，實際上是計算樣品對的半變差值，樣品對空間距離以及樣品對的法向量。通過CUDA平臺，編程實現(xiàn)GPU加速。

1）根據(jù)當(dāng)前硬件條件確定CUDA中的Block的尺寸，再根據(jù)數(shù)據(jù)的規(guī)模確定Block的數(shù)量；

在本次實驗中，每個線程塊使用512個左右的線程進行并行計算，而Block的數(shù)量可以使用文件數(shù)據(jù)量與每個線程塊線程之比來計算。

2）根據(jù)CUDA全局內(nèi)存的容量和數(shù)據(jù)的規(guī)模，確定數(shù)據(jù)分塊的大小和數(shù)目，以便完成多批次計算。具體計算方法如下：

EleNum=MaxPitch/（Z*sizeof（float））

BlkNum=（m*n+EleNum-1）/EleNum

其中MaxPitch是從CUDA函數(shù)接口獲取的允許函數(shù)最大pitch的拷貝尺寸。

3）使用CUDA函數(shù)在顯存中開辟數(shù)據(jù)計算空間和并行計算結(jié)果空間，將層之間的數(shù)據(jù)按照分塊形式依次從內(nèi)存?zhèn)魅腼@存中，并設(shè)計核函數(shù)進行計算；

4）使用CUDA函數(shù)將計算結(jié)果從顯存轉(zhuǎn)移至內(nèi)存，保存至文件中。

而整體樣品庫的建立，則是依據(jù)CUDA函數(shù)的輸出結(jié)果，依照樣品庫Sap=（SapBasic，SapVar，Sap?Length，SapAngle）的結(jié)構(gòu)，依次將樣品庫所需要的信息錄入。

3.2.3 輸出實驗變差函數(shù)計算結(jié)果

在以上前期工作全部完成之后，接下來就是進行調(diào)試輸出，即根據(jù)當(dāng)前約束規(guī)則，快速得到該約束規(guī)則下的實驗變差函數(shù)計算結(jié)果。

1）將約束規(guī)則的變化轉(zhuǎn)換成參數(shù)大小的改變；

2）將改變后的約束參數(shù)與樣品整體庫中的Sap對象進行約束規(guī)則比對，當(dāng)條件符合，則記錄下該對象；

3）把所有符合約束規(guī)則的對象進行統(tǒng)計，按照實驗變差函數(shù)的坐標(biāo)曲線進行結(jié)果輸出；

4）判斷當(dāng)前的結(jié)果是否符合下一步計算要求，如果符合則完成運算，如果不符合，重復(fù)步驟1）。

圖7 輸出實驗變差函數(shù)

4 試驗結(jié)果

4.1 算法比對

上面主要介紹了GPU加速的實驗變差函數(shù)優(yōu)化算法，并在理論上與兩種不同的計算方法進行了對比。通過分析可以發(fā)現(xiàn)，該算法比傳統(tǒng)算法在快速得出大規(guī)模數(shù)據(jù)實驗預(yù)期結(jié)果的速度上有著顯著的提高。因此我們結(jié)合兩個算法的特點，進行了算法的優(yōu)缺點分析：

1）傳統(tǒng)算法，實質(zhì)上是按需所求的進行計算，根據(jù)當(dāng)前約束條件，得到約束條件內(nèi)的樣品對，每確認一個樣品對，進行一次計算。該計算方法，沒有復(fù)雜的調(diào)用過程，運算簡單清晰，在不需要大量調(diào)試，數(shù)據(jù)量比較小的環(huán)境里，該算法計算速度很快，但該算法所需環(huán)境在真實工程應(yīng)用中效率會大大降低；

2）GPU加速算法，在傳統(tǒng)算法計算的基礎(chǔ)上，考慮了計算半變差值的并行特點，在該計算過程采用GPU加速。這種方法在傳統(tǒng)的方法基礎(chǔ)上，對計算過程進行了加速，在原始數(shù)據(jù)空間分布均勻的環(huán)境里，不需要進行多次約束就可以快速得出穩(wěn)定的實驗變差函數(shù)結(jié)果。

根據(jù)上面的分析，可以對以上兩種算法的優(yōu)點和缺點進行總結(jié)如表1所示。

表1 兩種實驗變差函數(shù)計算算法有缺點比較

4.2 實驗比對

本次試驗所采用的數(shù)據(jù)，設(shè)有試驗?zāi)M數(shù)據(jù)，真實礦山數(shù)據(jù)等，包括點狀樣品數(shù)據(jù)模型和線狀樣品數(shù)據(jù)模型，根據(jù)不同的規(guī)模劃分為5個數(shù)量級。

圖8 五組數(shù)據(jù)

將以上五組實驗數(shù)據(jù)分別用本文優(yōu)化算法比較，因穩(wěn)健的實驗變差函數(shù)結(jié)果需要專家進行每次計算結(jié)果的分析比對，故將實驗結(jié)果分為兩部分：首次計算與調(diào)試計算，可得出實驗結(jié)果如下：

表2 首次計算三種算法結(jié)果比對

由表2可以初步得知，在首次計算過程中，GPU加速算法和傳統(tǒng)算法，在數(shù)據(jù)量增加的過程中都有較明顯的優(yōu)勢。

由圖9也可以看出隨著數(shù)據(jù)量的增加，兩種加速算法比傳統(tǒng)算法優(yōu)勢越來越明顯，且單純GPU加速算法比整體到局部算法還稍有優(yōu)勢，當(dāng)然一次完整的實驗變差函數(shù)的計算過程，絕不可能只是一次計算就能得出正確結(jié)果，故需要進行后期的調(diào)試計算。

圖9 算法首次運行時間對比圖（單位：ms）

表3 調(diào)試過程一次計算耗費時間比對

則可得相應(yīng)的對比圖如圖10所示。

圖10 調(diào)試一次耗費時間對比圖（單位：ms）

一次完整的過程所耗費的總時間為首次計算時間與調(diào)試時間的總合。

即TA=TF+（TD1+TD1+TD1+…+TDN），其中TA為耗費的總時間，TF為首次計算時間，TDi為單次調(diào)試時間，調(diào)試次數(shù)N會隨著數(shù)據(jù)的改變而改變，一般情形N≥15。

表4 15次調(diào)試三種算法時間對比

從上表3、4可以看出，優(yōu)化過后的算法，在實驗變差函數(shù)的計算過程中有著顯著的優(yōu)勢，并且，當(dāng)數(shù)據(jù)量增大，優(yōu)勢會更加突出。

5 結(jié)語

綜上所述，本文提出了一個三維空間實驗變差函數(shù)優(yōu)化算法，并應(yīng)用到真實的地質(zhì)統(tǒng)計學(xué)應(yīng)用軟件中去，可以解決真實環(huán)境中實驗變差函數(shù)計算效率低下的問題，且通過多種數(shù)據(jù)下的實驗驗證了實驗的目的。

當(dāng)然本次算法也有不足之處，通過對約束分析，我們可以看到，每次進行一次新的約束，都可能存在著與前一次約束相同的樣品對重復(fù)計算，所以本文的算法雖然通過GPU加速，能夠提升計算效率，但仍因為重復(fù)的樣品對，導(dǎo)致了大量的時間與空間的浪費，所以如何能夠優(yōu)化這個過程，是解決這個問題的關(guān)鍵。一次性把所有樣品對的實驗變差值結(jié)果通過GPU加速算法計算出來，再進行實驗變差函數(shù)計算的約束分析之后，判定哪些計算結(jié)果是符合約束條件的樣品對是一大難題，因此以后的工作室設(shè)計一個便于后期進行約束篩選的樣品對庫來實現(xiàn)實驗變差函數(shù)的進一步優(yōu)化。