基于空間插值算法確定離子型稀土礦體的底板

謝文強， 龍平， 胡世麗

（江西理工大學土木與測繪工程學院，江西 贛州 341000）

離子型稀土礦是我國重要的礦產資源，素有“工業(yè)維生素”和“工業(yè)黃金”之美譽，對現(xiàn)代高技術產業(yè)發(fā)展具有重要的作用[1]。離子型稀土礦的成礦條件極其復雜，花崗巖、火山巖、凝灰?guī)r和混合巖等含稀土母巖，在溫暖濕潤的環(huán)境下，經物理、生物和化學等作用，母巖中的礦物相稀土轉化為離子相形態(tài)，隨著水流遷移，被黏土礦物吸附富集，形成離子型稀土礦[2]。目前，離子型稀土礦主要采用原地浸礦工藝進行開采[3]。原地浸礦是將浸礦劑溶液通過注液孔注入礦體，通過離子交換作用，將吸附在礦物表面的稀土離子解吸，形成稀土母液，稀土母液通過導流孔或巷道等收液工程收集，經過除雜、沉淀后實現(xiàn)稀土回收[4-6]。然而，該工藝在開采過程中常因礦體底板形狀不明，導致注液、收液工程設計不合理，常出現(xiàn)收液率低，造成資源浪費、收液工程施工遇石返工耽誤工期和反吸附嚴重等問題[7-9]。因此，準確確定離子型稀土礦體底板形狀，對合理設計注液與收液工程、稀土資源回收和環(huán)境保護都具有重要意義。

地質勘探工作是當前探明礦體底板形狀的有效手段[10]，現(xiàn)場采用贛南鉆進行勘探，通過鉆孔打穿地表直至見水或見石處，以確定該勘探孔孔底的性質。然而，實際相鄰勘探孔的間距較大，無法反映礦體底板局部區(qū)域的形狀特征。因此，如何將有限勘探孔數(shù)據實現(xiàn)由點向面的轉化，是確定礦體底板的有效方法。高原等[11]在徑向基函數(shù)插值法的基礎上，建立顧及地形特征的定權方法，有效反映了局部地形細節(jié)特征。申靜等[12]將普通克里金法應用于海底地形的可視化表達。ZHANG等[13]提出的基于地形粗糙度的克里金法，相比傳統(tǒng)的插值方法能顯著提高插值效率與精度。李曉軍等[14]對煤層底板進行貝葉斯地統(tǒng)計模型插值構造，認為該方法在地質層面高程預測方面具有一定的優(yōu)越性。蔣偉達等[15]對埕島海域地形演變進行多種插值方法分析，發(fā)現(xiàn)反距離權重法能更真實地反映該地區(qū)的水深地形。CHEN等[16]采用5種插值方法對不同粗糙程度的地形進行了評估，發(fā)現(xiàn)反距離權重法輸出表面粗糙，難以體現(xiàn)地形細節(jié)特征。由此可見，空間插值方法在地質領域已被廣泛運用，但基于空間插值算法確定離子型稀土礦體底板形狀的研究較少。

本文采用常用的反距離權重法、徑向基函數(shù)插值法、局部多項式法、普通克里金法和經驗貝葉斯克里金法共5種空間插值方法，對離子型稀土礦體底板形狀進行空間插值，利用獨立數(shù)據集比較各方法的精度，分析稀土礦體底板的空間特征，為合理設計離子型稀土礦注液與收液工程提供指導。

1 實驗部分

1.1 研究區(qū)域概況

試驗場地位于福建省某離子型稀土礦山，該區(qū)域氣候溫和，雨水充沛，屬于亞熱帶季風性氣候，年均氣溫約18 ℃，年均降水量約1 700 mm[17]。該區(qū)域礦床母巖多以花崗巖為主，形成了儲量較大、品位較高的離子型稀土礦，是我國南方重要的離子型稀土礦產資源富集區(qū)。

1.2 數(shù)據采集

勘探網度設置為約20 m×20 m，勘探孔直徑設置為100 mm，采用贛南鉆進行勘探工作。以潛水位或半風化層作為勘探孔的終孔依據，當勘探孔終孔至潛水位時標示孔底信息為見水，當勘探孔終孔至半風化層時標示孔底信息為見石，測量勘探孔深度。若某一勘探孔孔深明顯小于周邊的勘探孔，該位置可能遇到孤石，在附近再選一個位置進行勘探。該礦塊共有440個見石勘探孔和39個見水勘探孔。采用GPS-RTK測量勘探孔孔口坐標，所有勘探孔的空間位置如圖1所示。以見石勘探孔為研究對象，孔口高程減去相應勘探孔深度，即見石勘探孔孔底的高程。

圖1 福建省某離子型稀土礦勘探孔分布Fig.1 Distribution map of exploratory holes of an ionic rare earth ore in Fujian province

1.3 研究方法

采用泰森多邊形聚類方法[18]對440個見石勘探孔孔底坐標數(shù)據進行離群檢驗，去除離群數(shù)據?；谀m指數(shù)[19]和K-S檢驗[20]對處理后的數(shù)據進行空間自相關性分析和正態(tài)性檢驗，判斷是否具有正態(tài)性，能否進行地統(tǒng)計學方法分析，避免計算結果出現(xiàn)比例效應[18]。莫蘭指數(shù)的取值范圍為[-1,1]，越趨近于1，說明數(shù)據的空間自相關性越好。通過K-S檢驗的P值判斷數(shù)據是否服從正態(tài)分布，若P＞ 0.05，則服從正態(tài)分布；若P＜ 0.05，則不服從正態(tài)分布，需對數(shù)據進行轉換，使其符合或基本符合正態(tài)分布。

分別隨機取55%、65%、75%、85%和95% 5組數(shù)據作為訓練點（數(shù)據點數(shù)按四舍五入取整），以其他數(shù)據作為檢驗點，通過反距離權重法[21]、徑向基函數(shù)插值法[22]、局部多項式插值法[23]、普通克里金法[24]（半變異函數(shù)為穩(wěn)定函數(shù)）和經驗貝葉斯克里金法[14]共5種空間插值方法分析訓練點，預測檢驗點處的見石高程，以標準差（σ，式（1））、平均絕對誤差（MAE，式（2））和均方根誤差（RMSE，式（3））作為評價指標，量化檢驗點處高程預測值與實際值的誤差。為了減少訓練點選擇的隨機性對誤差分析的影響，隨機12次訓練點，取12次檢驗點誤差指標參數(shù)（σ、MAE、RMSE）的平均值作為每種方法的誤差，根據誤差值評價5種方法。

式（1）—式（3）中：N為樣本數(shù)；Zcheck,i為實測高程；為預測高程為實測高程的平均值。

2 結果與討論

2.1 數(shù)據點檢驗

440個見石勘探孔的聚類泰森多邊形如圖2所示，若某個多邊形單元的級區(qū)間與其相鄰單元的級區(qū)間都有所不同，則該單元的勘探孔為離群點，將離群單元用黑色表示，也繪制于圖2中，本礦塊的見石勘探孔共有15個離群點。將15個離群點去除，實際用于分析的見石勘探孔共有425個，離群處理后的見石勘探孔的聚類泰森多邊形如圖2（b）所示，由圖2（b）可知，經過一次離群處理的數(shù)據，已經無新的離群點。

圖2 離群值去除前（a）和去除后（b）的泰森多邊形示意Fig.2 Schematic diagram of Tyson polygon before （a） and after （b） outlier removal

425個見石勘探孔的數(shù)據的莫蘭指數(shù)為0.93，說明數(shù)據具有較高的空間自相關性，可應用地統(tǒng)計學方法進行空間插值分析。采用K-S檢驗分析勘探孔數(shù)據，得到P值為0.01（＜ 0.05），表明所分析的425個見石勘探孔數(shù)據不符合正態(tài)分布。運用對數(shù)變換對數(shù)據進行處理，在采用K-S檢驗進行分析，得到P值為0.06（＞ 0.05）。對數(shù)變換前，見石勘探孔高程頻率分布直方圖的偏度系數(shù)[25]為0.16；對數(shù)變換后，偏度系數(shù)為0.05。由此可見，采用對數(shù)轉換處理后，可以提高數(shù)據的正態(tài)性，見石勘探孔高程數(shù)據可以視為正態(tài)分布（見圖3 ），避免采用地統(tǒng)計學方法分析時的比例效應，保證了空間插值結果的精度。

圖3 采用對數(shù)處理前（a）和后（b）見石勘探孔的頻率分布直方圖Fig.3 Frequency distribution histogram of stone exploratory holes before （a） and after （b） logarithmic processing

2.2 插值方法比較

從圖4可知，5種插值方法的標準差（σ）、平均絕對誤差（MAE）和均方根誤差（RMSE）驗證指標結果并不相同。隨著檢驗點個數(shù)的增加，不同插值方法的精度均有所降低，總體上σ、MAE和RMSE隨著檢驗點個數(shù)的增加而逐漸提升。在檢驗點個數(shù)相同的情況下，雖然5種插值方法在不同檢驗點個數(shù)的預測誤差范圍內均有所差異，但總體上徑向基函數(shù)插值法的預測誤差范圍最小，其σ和RMSE均低于6.1，MAE低于4.9。由此可見，徑向基函數(shù)插值法在插值精度上優(yōu)于其他插值方法，預測效果更佳。

圖4 檢驗點誤差指標參數(shù)結果：（a） 標準差σ；（b） 平均絕對誤差MAE；（c） 均方根誤差RMSEFig.4 Parameter results of error index of checkpoints：（a） σ；（b） MAE；（c） RMSE

不同插值方法得到的稀土礦體底板空間形狀格局大體相似，總體呈西北至東南方向逐漸降低的格局，但局部區(qū)域的細節(jié)特征及曲線平滑程度存在一定差異（見圖5）。通過對比各插值方法的檢驗點誤差分布（見圖6），5種插值方法對區(qū)域 Ⅰ 內的見石勘探孔高程預測誤差范圍較大，主要是因為1號見石勘探孔（Z1=326.72 m）與相對鄰近2號和3號（Z2=338.18 m，Z3=341.69 m）見石勘探孔實測高程相差頗遠，極值點嚴重影響插值方法對局部區(qū)域的插值精度。各插值方法中，反距離權重法出現(xiàn)的頻率相對更多，這是因為在插值過程中，該方法忽視各訓練點之間的空間自相關性，僅以待測點附近訓練點的距離大小，作為權重系數(shù)進行賦值，因而插值生成的擬合曲面極易受到極值點或點集群的影響，在局部區(qū)域內呈現(xiàn)孤立的“牛眼狀”分布（見圖6（a）），導致局部區(qū)域插值精度降低。

圖5 基于5種插值方法生成的稀土礦體底板等值線：（a）反距離權重法；（b） 徑向基函數(shù)插值法；（c） 局部多項式插值法；（d） 普通克里金法；（e） 經驗貝葉斯克里金法Fig.5 Contour map of bottom plate of rare earth ore body based on five interpolation methods：（a）inverse distance weighting method；（b） radial basis function method；（c） local polynomial method；（d） ordinary Kriging method；（e） empirical Bayesian Kriging method

圖6 檢驗點的誤差分布：（a）反距離權重法；（b） 徑向基函數(shù)插值法；（c） 局部多項式插值法；（d） 普通克里金法；（e） 經驗貝葉斯克里金法Fig. 6 Error distribution of inspection point：（a）inverse distance weighting method；（b） radial basis function method；（c） local polynomial method；（d） ordinary Kriging method；（e） empirical Bayesian Kriging method

高程變化明顯的區(qū)域也是各插值產生偏差的原因。在區(qū)域 Ⅱ 中，4號、5號和6號見石勘探孔高程（Z4=368.36 m，Z5=370.83 m，Z6=378.05 m）由東南至西北方向逐漸升高，該區(qū)域各插值方法的預測誤差較其他平穩(wěn)區(qū)域相差較大，表明各插值方法擬合形成的曲面均受到高程變化而影響精度。經驗貝葉斯克里金法雖然通過利用訓練點數(shù)據能模擬多個半變異函數(shù)進行插值分析，以此提高插值精度，但由圖5（e）可知，該方法生成的預測表面過于平滑，難以反映局部高程區(qū)域的細節(jié)變化，而相比于插值精度最高的徑向基函數(shù)插值法，該方法具備對訓練點數(shù)據擬合程度高、生成預測表面平滑能力強等優(yōu)點，能凸顯礦體底板高程的變化幅度，對高程變化波動較大的數(shù)據，插值精度通常較高，能較好地反映底板局部變化特征與趨勢。然而，徑向基函數(shù)并非總體最優(yōu)，相對局部訓練點空白區(qū)域 Ⅲ 而言，該區(qū)域徑向基函數(shù)插值法相比反距離權重法誤差范圍更大，表明在處理無訓練點區(qū)域，徑向基函數(shù)插值法插值精度相對較低。

此外，四周邊緣區(qū)域缺少訓練點數(shù)據作為參考，對稀土礦體底板預測空間分布也會造成較大的影響。從區(qū)域 Ⅳ 可知，8號和9號見石勘探孔分布在礦區(qū)的右邊緣處，相比于中間區(qū)域，該區(qū)域訓練點進行插值時，插值精度會受到周邊訓練點數(shù)目的影響。通過對比各插值方法的誤差分布，局部多項式插值法和普通克里金法因邊緣區(qū)域缺少足夠的訓練點數(shù)據，在區(qū)域的邊緣化問題中，易產生更大的誤差，造成周邊區(qū)域插值精度降低，這與費坤等[21]的研究相吻合。

由此可見，雖然在局部訓練點空白區(qū)域，徑向基函數(shù)插值法的插值精度低于反距離權重法，但總體在插值精度或是插值效果上都優(yōu)于其他方法，采用徑向基函數(shù)插值法是確定離子型稀土礦體底板空間形狀較好的方法。

3 結 論

本文基于有限的見石勘探孔高程數(shù)據，采用空間插值算法確定離子型稀土礦體底板形狀，并通過標準差、平均絕對誤差和均方根誤差分析各種方法的精度，在此過程中得到以下結論：

1）去除離群點和對數(shù)處理后，見石勘探孔數(shù)據的莫蘭指數(shù)和P值分別為0.93和0.06，說明對數(shù)處理后數(shù)據具有較好的空間自相關性和正態(tài)分布特性，能用地統(tǒng)計學插值方法進行分析。

2）隨著檢驗點的增多，檢驗點的σ、RMSE和MAE均呈增加趨勢，當檢驗點占比超過35%時，σ、RMSE和MAE量化參數(shù)趨于穩(wěn)定。

3）徑向基函數(shù)插值法的σ和RMSE均低于6.1以下，MAE低于4.9，總體插值精度和插值效果都優(yōu)于其他插值方法，能更好地呈現(xiàn)離子型稀土礦體底板總體形狀和局部變化的細節(jié)特征。