張燕，張素榮

（中國(guó)地質(zhì)調(diào)查局天津地質(zhì)調(diào)查中心,天津 300170）

離散數(shù)據(jù)地球化學(xué)編圖中克里格法的應(yīng)用—以津巴布韋超低密度地球化學(xué)圖為例

張燕，張素榮

（中國(guó)地質(zhì)調(diào)查局天津地質(zhì)調(diào)查中心,天津 300170）

自由網(wǎng)離散數(shù)據(jù)的地球化學(xué)編圖以任意方式進(jìn)行網(wǎng)格化成圖存在明顯不確定性，會(huì)導(dǎo)致地球化學(xué)特征的改變或歪曲?？死锔穹ㄌ峁┝艘环N更科學(xué)合理的空間估計(jì)插值方法。筆者運(yùn)用該方法對(duì)津巴布韋超低密度地球化學(xué)調(diào)查數(shù)據(jù)進(jìn)行了空間分析并編制了地球化學(xué)圖件，研究和討論了自由網(wǎng)離散數(shù)據(jù)的地球化學(xué)編圖方法，評(píng)價(jià)了采用克里格法編制地球化學(xué)圖的效果。

編圖方法；克里格法；變差函數(shù)；津巴布韋；地球化學(xué)圖

地球化學(xué)等值線圖的應(yīng)用非常廣泛，通常是對(duì)地球化學(xué)填圖形成的數(shù)據(jù)進(jìn)行網(wǎng)格化插值來(lái)實(shí)現(xiàn)的。地球化學(xué)編圖所依賴(lài)的數(shù)據(jù)分布及其結(jié)構(gòu)特征對(duì)圖件編制有重要影響。目前對(duì)自由網(wǎng)地球化學(xué)數(shù)據(jù)的插值估計(jì)大多采用了距離倒數(shù)加權(quán)插值法,其在算法上比較簡(jiǎn)單,但是沒(méi)有考慮到數(shù)據(jù)點(diǎn)間的空間相關(guān)性和合理的半徑,計(jì)算和編圖結(jié)果一般是不確定的，甚至可能是錯(cuò)誤的。自由網(wǎng)地球化學(xué)數(shù)據(jù)的插值估計(jì)形成的曲面應(yīng)能最好地?cái)M合實(shí)測(cè)數(shù)據(jù)。地球化學(xué)調(diào)查數(shù)據(jù)的空間變化規(guī)律應(yīng)從數(shù)據(jù)本身分析，研究數(shù)據(jù)變化的每個(gè)變量的擬合變差函數(shù)的信息。地球化學(xué)調(diào)查數(shù)據(jù)的插值估計(jì)通?？梢圆捎每死锔穹ǎ@是一種有效的空間估計(jì)插值方法，在很大程度上依賴(lài)于數(shù)據(jù)模型的空間相關(guān)性。

本文以津巴布韋低超密度地球化學(xué)圖的編制為例，對(duì)離散數(shù)據(jù)的區(qū)域地球化學(xué)等值線圖編圖方法進(jìn)行了探索和研究，編制了一套自由網(wǎng)離散數(shù)據(jù)條件下的區(qū)域地球化學(xué)圖，并對(duì)克里格法的效果及其合理性進(jìn)行了討論。

1 克里格法的基本原理

克里格插值（Kriging）是建立在變異函數(shù)空間分析基礎(chǔ)上，對(duì)有限區(qū)域內(nèi)的區(qū)域化變量取值進(jìn)行無(wú)偏最優(yōu)估計(jì)的一種方法[1]。這是1951年南非采礦工程師D.G.K rige提出的對(duì)礦體品位和儲(chǔ)量的一種估值方法。1971年法國(guó)地統(tǒng)計(jì)學(xué)家Matheron把這種方法命名為Kriging，并在此基礎(chǔ)上提出了區(qū)域化變量理論，使傳統(tǒng)的地學(xué)方法與統(tǒng)計(jì)方法相結(jié)合，理論化、系統(tǒng)化，形成了完整的方法體系。

1.1 克里格插值的關(guān)鍵問(wèn)題

克里格法的基本原理是根據(jù)相鄰變量的數(shù)值，利用變差函數(shù)所揭示的區(qū)域化變量的內(nèi)在聯(lián)系來(lái)估計(jì)空間數(shù)值。按照空間場(chǎng)是否存在漂移可將克里格法（Kriging）分為普通Kriging法和泛Kriging法,其中以普通Kriging法在地球化學(xué)編圖中最為常用。

克里格插值是根據(jù)待插值點(diǎn)與其臨近實(shí)測(cè)數(shù)據(jù)點(diǎn)的空間位置，對(duì)待插值點(diǎn)的變量值進(jìn)行無(wú)偏最優(yōu)估計(jì)，通過(guò)生成一個(gè)關(guān)于區(qū)域化變量的克里格插值來(lái)表達(dá)研究區(qū)域的原始數(shù)值，其基本公式為式如下：

式中：Z（x0）表示未知樣點(diǎn)的值；Z（xi）表示未知樣點(diǎn)周?chē)囊阎獦颖军c(diǎn)的值；N為已知樣本點(diǎn)的個(gè)數(shù)；λi為第i個(gè)樣本點(diǎn)的權(quán)重，但權(quán)重一般并非距離倒數(shù)的權(quán)重，而是某種擬合數(shù)學(xué)模型。它的確定可以通過(guò)半方差的分析獲取，根據(jù)統(tǒng)計(jì)學(xué)上無(wú)偏和最優(yōu)的要求，利用拉格朗日極小化原理，可推導(dǎo)出權(quán)重值和半方差之間的公式[2]。

從求解估計(jì)值Z(X0)的過(guò)程可見(jiàn)，Kriging插值的關(guān)鍵是如何獲取變差函數(shù)的最佳計(jì)算公式。這涉及兩個(gè)主要問(wèn)題,變異函數(shù)理論模型的選取和模型參數(shù)的估計(jì)[3]。

1.2 半方差函數(shù)

所謂半方差函數(shù)γ(h)就是兩點(diǎn)間差值方差的一半，是用來(lái)描述地球化學(xué)特性空間連續(xù)變異的一個(gè)連續(xù)函數(shù)，反映了不同距離觀測(cè)值之間的地球化學(xué)特性的變化，也稱(chēng)為變異函數(shù)。它是地質(zhì)統(tǒng)計(jì)學(xué)中研究土壤、巖石或水系沉積物地球化學(xué)數(shù)據(jù)特征值變異性的關(guān)鍵函數(shù)。其表示形式如下：

式中：γ(h)為實(shí)驗(yàn)半方差函數(shù)；H=xi+h，h為步長(zhǎng)，即在一定方向上，距離為h的矢量；Z(xi+h)，Z(xi)為某變量的測(cè)定值；N(h)是以h為間距的所有觀測(cè)點(diǎn)的成對(duì)數(shù)目。

總體來(lái)講，距離愈近的已測(cè)點(diǎn)間和距離較遠(yuǎn)的已測(cè)點(diǎn)間相比，前者具有較小的平方差。每一樣點(diǎn)對(duì)被劃出后（分組后），就可以用一個(gè)模型來(lái)擬合它們?；谧钚《朔ㄅc最大似然法得到擬合最好的表面。理論上變異函數(shù)模型比較多,常用的有線性模型、指數(shù)模型、高斯模型和球狀模型等[3]。有幾個(gè)重要的參數(shù)可用來(lái)描述這些模型：塊金值、基臺(tái)值、自相關(guān)閾值（變程）。實(shí)際應(yīng)用的時(shí)候可以在高級(jí)參數(shù)里面對(duì)其進(jìn)行調(diào)整以減少誤差方差。

2 工作步驟

克里格插值方法需要完成以下主要步驟：

（1）分析數(shù)據(jù)的分布特征，揭示數(shù)據(jù)的空間規(guī)律。對(duì)數(shù)據(jù)的分布進(jìn)行分析，必要時(shí)進(jìn)行數(shù)據(jù)轉(zhuǎn)換，使數(shù)據(jù)符合正態(tài)分布的前提假設(shè)。正態(tài)分布轉(zhuǎn)換能夠改善插值結(jié)果。

（2）分析數(shù)據(jù)的空間變化趨勢(shì)。總體上,可以近似認(rèn)為變異函數(shù)存在幾何異向性。進(jìn)一步,需要將各向異性轉(zhuǎn)化為各向同性[4]。

（3）確定變異方差的函數(shù)類(lèi)型。生成變異函數(shù)和協(xié)方差函數(shù)，用于估算樣點(diǎn)值間的統(tǒng)計(jì)相關(guān)（空間自相關(guān)）。變異函數(shù)和協(xié)方差函數(shù)取決于自相關(guān)模型（擬合模型）。分別采用線性模型、指數(shù)模型、高斯模型和球狀模型進(jìn)行擬合,確定影響方差變異函數(shù)分析有四個(gè)基本參數(shù)：塊金（C0），基臺(tái)(C0＋C)，變程（a）和方差變異函數(shù)類(lèi)型（擬合模型）。這需要進(jìn)行多次模擬運(yùn)算，才能得到合適的結(jié)果。

（4）計(jì)算網(wǎng)格化的預(yù)測(cè)值。依據(jù)變程選擇合適的搜索半徑和臨近數(shù)據(jù)點(diǎn)數(shù)，選取合適的插值方法進(jìn)行網(wǎng)格化插值計(jì)算。

（5）結(jié)果分析和檢驗(yàn)。對(duì)擬合的連續(xù)曲面上的估計(jì)值與實(shí)測(cè)值的方差進(jìn)行檢驗(yàn)。一般而言，如果克里格插值結(jié)果具有最小的方差，而且交叉檢驗(yàn)中誤差的方差最小，那么該結(jié)果就是合適的。

在通常所研究的空間尺度情況下，半方差函數(shù)值都隨著樣點(diǎn)間距的增加而增大，并在樣點(diǎn)間距大于相關(guān)距離（a）時(shí)，該值趨于一個(gè)基本穩(wěn)定的常數(shù)(C0＋C)。

3 應(yīng)用實(shí)例-津巴布韋超低密度地球化學(xué)圖的編制

（1）樣品與數(shù)據(jù)

津巴布韋全國(guó)低密度地球化學(xué)調(diào)查的采樣是根據(jù)津巴布韋全國(guó)地形地貌和水系分布，并結(jié)合交通條件，采用水系沉積物樣品以超低的采樣密度在高級(jí)別的水系或河流中系統(tǒng)布置和采集的（圖1）。采樣力求控制面積均等和兼顧樣品分布的均勻性。采樣密度按1個(gè)樣品控制800～1000 km2的匯水面積，使每個(gè)采樣點(diǎn)都能有效地控制相應(yīng)的匯水面積。野外共采268件樣品，實(shí)際控制30余萬(wàn)平方千米的面積。選用多種分析方法的配套方案，對(duì)41種元素和化合物進(jìn)行了定量分析。獲得了每個(gè)采樣點(diǎn)的平面位置坐標(biāo)、高程和41種元素的含量數(shù)據(jù)[5]。

（2）數(shù)據(jù)統(tǒng)計(jì)和變差函數(shù)模型的選擇

津巴布韋超低密度區(qū)域地球化學(xué)調(diào)查數(shù)據(jù)的統(tǒng)計(jì)分析和變差函數(shù)模型的檢驗(yàn)采用地統(tǒng)計(jì)學(xué)軟件GS+7.0進(jìn)行分析。通過(guò)計(jì)算機(jī)軟件對(duì)數(shù)據(jù)的分布進(jìn)行分析，分析數(shù)據(jù)的空間變化趨勢(shì)，采用不同的變差函數(shù)模型進(jìn)行擬合計(jì)算。選擇決定系數(shù)（R）最大的模型，即選用的模型和參數(shù)對(duì)離散數(shù)據(jù)的擬合程度高，網(wǎng)格化計(jì)算出來(lái)的值與實(shí)際測(cè)量值的偏差最小。經(jīng)過(guò)計(jì)算機(jī)處理和計(jì)算，41種元素最終采用的網(wǎng)格化變差函數(shù)模型類(lèi)型為Ag、Au、Bi、Cao、Cd、Co、Cu、K2O、Li、Mo、Nb、P、Pb、Sb、Ti、W、Pt、Pt的變差函數(shù)為球狀模型；A l2O3、Be、Cr、MgO、Ni等的變差函數(shù)為高斯模型；Zn、As、B、Fe2O3、Sn、V的變差函數(shù)模型為線型；Ba、F、Hg、La、Mn、Na2O、SiO2、Sr、Th、U、Y、

Zr的變差函數(shù)模型為冪指數(shù)型。其中部分元素的函數(shù)模型的參數(shù)計(jì)算結(jié)果見(jiàn)表1。

（3）圖件編制

津巴布韋全國(guó)低密度地球化學(xué)圖使用的軟件版本為中地?cái)?shù)碼開(kāi)發(fā)的MapGIS6.7軟件，采用Kriging對(duì)離散數(shù)據(jù)進(jìn)行網(wǎng)格化插值成圖。根據(jù)變差函數(shù)模型的分析和檢驗(yàn)，成圖時(shí)對(duì)41種元素的變差函數(shù)模型分別采用高斯模型、球狀模型、線形模型、冪指數(shù)模型等變差函數(shù)擬合模型。插值網(wǎng)格化間距依據(jù)數(shù)據(jù)最小間距確定為8 km；搜索半徑依據(jù)各元素?cái)?shù)據(jù)模型的變程確定，最大搜索半徑100 km（圖2～5）。

4 結(jié)論

（1）在克里格插值中，權(quán)重不僅考慮了已知點(diǎn)與插值點(diǎn)間的距離，而且考慮了己知點(diǎn)的位置和屬性值整體的空間分布和格局，插值中的權(quán)重取決于已知點(diǎn)的擬合模型和已知點(diǎn)到插值點(diǎn)的距離以及插值點(diǎn)周?chē)囊阎獦狱c(diǎn)的空間關(guān)系等多重因素。

（2）克里格分析方法需要一定的數(shù)據(jù)量，計(jì)算中的樣本不夠可能影響效果?？死锔癫逯嫡`差隨塊金值的增大而增大；基臺(tái)值對(duì)結(jié)果影響不大，確定最佳變程值和最佳擬合函數(shù)是重要的，應(yīng)選擇決定系數(shù)

（R2）較高的函數(shù)。

圖1 津巴布韋超低密度區(qū)域地球化學(xué)調(diào)查采樣點(diǎn)位圖Fig.1 Sam p ling loca tions o f the u ltra-low density regiona l geochem ica l su rvey in Zim babw e

表1 津巴布韋超低密度區(qū)域地球化學(xué)調(diào)查數(shù)據(jù)變差函數(shù)模型Table 1 Va ria tion func tion m ode lo f u ltra-low density regiona lgeochem ica l su rvey da ta in Zim babw e

（3）津巴布韋超低密度區(qū)域地球化學(xué)圖件編制的效果是可以肯定的，總體上反映了區(qū)域地球化學(xué)背景特征的實(shí)際情況，清楚地指示了若干地球化學(xué)省、地球化學(xué)帶的空間展布規(guī)律，為區(qū)域地球化學(xué)研究提供了可信的基礎(chǔ)圖件。

圖2 津巴布韋Al2O3地球化學(xué)圖Fig.2 Geochem ica lm ap o f Alum inum Oxide in Zim babwe

圖3 津巴布韋K2O地球化學(xué)圖Fig.3 Geochem ica lm ap of Po tassium Oxide in Zimbabwe

圖4 津巴布韋B地球化學(xué)圖Fig.4 Geochem ica lmap o f Boron in Zimbabw e

圖5 津巴布韋Cr地球化學(xué)圖Fig.5 Geochem ica lmap o f Chrom ium in Zim babwe

[1]侯景儒，黃競(jìng)先.地質(zhì)統(tǒng)計(jì)學(xué)的理論與方法[M].北京:地質(zhì)出版社，1990.

[2]王艷妮，謝金梅，郭祥.ArcGIS中的地統(tǒng)計(jì)克里格插值法及其應(yīng)用[J].軟件導(dǎo)刊，2008，7(12):36-38.

[3]曾懷恩,黃聲享.基于Kriging方法的空間數(shù)據(jù)插值研究[J].測(cè)繪工程，2007，16(5):5-13.

[4]孫洪泉.地質(zhì)統(tǒng)計(jì)學(xué)及其應(yīng)用[M].北京:中國(guó)礦業(yè)大學(xué)出版社，1990.

[5]Gengxin Zhao，F(xiàn)uqing He,etal.Ultra-low density geochemicalmapping in Zimbabwe[J].Journal of Geochemical Exploration，2014.

Abstract:Theweightsof evidencemodel isan importantmethod for themetallogenetic prognosis.It isa statisticalmethod of geoscience based on two valued image.The metallogenetic prognosis based on GIS using the weights of evidencemodel can be done through overlaying the informations of geoscience related w ith them ineral resources.Themetallogenetic prognosis based on GIS using theweightsof evidencemodel is based on the data,and realized easily.Using datum in North Daxinganling,weights of evidence predictionmodel of gold show preferred precision for predictingm inerals resources.Information to create new preferable compositive variables was integrated using theARCGISsystem.

Key words:GIS;weightof evidencemethod;metallogenetic prediction;evidence layer

更正

本人及合作者在本刊2012年第2期發(fā)表的《華北克拉通南緣洛峪群－汝陽(yáng)群屬于中元古界長(zhǎng)城系—河南汝州洛峪口組層凝灰?guī)r鋯石LA-MC-ICPMSU-Pb年齡的直接約束》一文中，由于我們的測(cè)年材料為左景勛等（1997）文章所報(bào)道汝州陽(yáng)坡村附近洛峪口組內(nèi)的夾層狀層凝灰?guī)r，并誤以為我們所觀察剖面就是他們文章報(bào)道的“陽(yáng)坡剖面”，因此在拙文中直接引用了他們的剖面圖，并根據(jù)我們的觀察和區(qū)域資料，將其上覆層位的時(shí)代歸屬由原來(lái)的“寒武系”修改為“震旦系”。最近，經(jīng)與左景勛先生交流并再次觀察后確認(rèn)，我們這一修改其實(shí)有欠妥之處。實(shí)際上，拙文所報(bào)道的測(cè)年剖面是近年新開(kāi)劈的通往陽(yáng)坡村的鄉(xiāng)村公路，位于左景勛等（1997）“陽(yáng)坡剖面”的東側(cè)。雖然兩者相隔不足千米且完全平行，但由于兩條剖面的起點(diǎn)及終點(diǎn)位置畢竟不同等原因，拙文的剖面頂部及附近地層為“震旦系”，與左景勛等（1997）所報(bào)道“寒武系”應(yīng)有所不同。特此更正并向左先生及其合作者、向廣大讀者以及《地質(zhì)調(diào)查與研究》編輯部致歉！

蘇文博（中國(guó)地質(zhì)大學(xué)（北京）地學(xué)院）

2014年12月3日

Application of K rigingon Geochem icalM appingM ethod w ith DiscreteData:Illustrated by theU ltra-low Density Geochem icalM apping in Zimbabwe

ZHANG Yan,ZHANG Su-rong
(Tianjin Center,ChinaGeologicalSurvey,Tianjin 300170,China)

There are obvious defects and uncertainty on the geochem icalmapping w ith freemeshing and discrete data by grid interpolation.Itwill lead that thegeochem ical characteristics are changed or distorted.Krigingmethod provides amore scientific and reasonable interpolationmethod forestimating the data.We discuss the application of the Kriging method on the geochem icalmappingmethod w ith discrete data by spatial analysis and geochem icalmapping of ultra-low density geochem ical data in Zimbabwe,and assessed the effect of application about the Krigingmethod in thegeochem icalmapping.

mapping;Kriging;variation function;Zimbabwe;geochem icalmap

MineralResources Evaluation of Comprehensive In form ation by Weightsof EvidenceM ethod in North Daxingan ling A rea

TENG Fei,ZHANG Yan,ZHANG Su-rong,XING Yi
(Tianjin Center,ChinaGeologicalSurvey,Tianjin 300170,China)

P628

：A

：1672－4135（2014）04－0283-05

2014-09-02

中國(guó)地質(zhì)調(diào)查局承擔(dān)的中國(guó)與津巴布韋的合作項(xiàng)目（商援[2007]573號(hào)）

張燕（1962-），女，工程師，1988年畢業(yè)于西安地質(zhì)學(xué)院制圖專(zhuān)業(yè)，現(xiàn)從事地質(zhì)信息專(zhuān)業(yè)，Email:1342052428@qq.com。