Excel在地球化學數據處理中的高級應用

譚親平，夏 勇，謝卓君，閆 俊

（1．中國科學院地球化學研究所，貴陽 550002；2．中國科學院大學，北京 100049）

Excel在地球化學數據處理中的高級應用

譚親平1，2，夏 勇1*，謝卓君1，2，閆 俊1，2

（1．中國科學院地球化學研究所，貴陽 550002；2．中國科學院大學，北京 100049）

Excel強大的圖表繪制和數據計算能力，為地球化學數據的處理提供了便利。在研究地球化學數據處理原理的基礎之上，詳細解釋了三角圖解，頻率直方圖，概率格紙圖解法求異常下限，多重分形法計算異常下限，R型聚類分析和判別分析的具體計算步驟。這些方法有利于地質科技方面的研究。

三角圖解；異常下限；直方圖；聚類分析；判別分析

0 引言

在地質工作中，常常需要計算各種數據，繪制各種圖表［1－3］，目前尚無統一的軟件能夠滿足地質工作當中所有的數據處理的需求。Excel為大眾軟件，地質科技人員較為熟悉，并有一定的操作經驗。它的強大的圖表繪制和數據計算能力基本能滿足地質工作中的數據處理需求。作者以Microsoft Excel 2010為界面，以文獻和作者在工作中的數據為例，詳細解釋各種地質數據處理的原理和具體操作步驟，它適用于基層地質科技人員的工作性質和工作條件。本研究涉及了Excel的高級用法，需要了解相關的多元統計學知識，并能熟練操作Excel才能合理地進行運用。

1 三角圖解

Microsoft Excel已經提供了大量的圖表類型，但在地學中需要經常應用的三角圖解卻沒有提供。在Excel中，可以通過將三角坐標轉化為直角坐標的辦法來“迂回”實現。三角圖解有a、b、c三個軸，它們的值范圍都是0～100，且滿足a＋b＋c＝100。通過坐標變換可將三維坐標系（a、b、c）轉換為二維直角坐標系（X，Y），然后再利用轉換得來的直角坐標繪制散點圖即可實現三角圖解的繪制。

圖1中直角坐標的原點設為BC的中點，任意點D在三角圖解中的坐標設為（a，b，c）。根據幾何學原理，任意點D在直角坐標中的坐標是：X＝（b＋c）／2－b，Y＝sin60°*a。設A、B、C三個端點的坐標分別為：A（100，0，0）、B（0，100，0）、C（0，0，100），根據上面的公式計算出它們的直角坐標分別為：A（0，86．60254）、B（－50，0）、C（50，0）。在Excel中利用三個端點的直角坐標，每兩個點繪制帶直線的散點圖，就能得到三角圖解的外框。將需要處理的數據，根據同樣的公式，也轉換成直角坐標，然后添加為散點圖即可。

下面通過一個實際的例子來說明三角圖解的繪制過程。表1中的原始數據來自文獻［4］。以表1中的第2排數據為例，在D2中輸入函數：＝100／（A2＋B2＋C2），即可得到比例因子，再將每個變量乘以比例因子，可在E2：G2中得到歸一化后的值。歸一化后使三個變量之和等于100。在H2和I2中分別輸入函數：＝（F2＋G2）／2－F2和＝E2 *√3／2得到相對應的直角坐標，將得到的直角坐標添加到前面已生成的三角形外框中，在生成的圖中刪除直角坐標軸，插入三個文本框，在文本框中填入三個端元符號。圖2是應用本文的方法得到的圖，它與文獻［4］中的圖重合。該方法的關鍵為三角坐標向直角坐標的轉換原理。

表1 三角圖解數據表Tab．1 Data table of triangular chart

圖1 三角坐標轉直角坐標原理圖Fig．1 Schematic diagram of traingular coordinate transform into rectangular coordinate

2 頻率直方圖

頻率直方圖在Excel中的繪制過程，可利用實例來說明。假如有表2的A列中的一組數據需要繪制頻率直方圖，首先在B3和B7中分別輸入函數：＝MIN（A2：A13）和＝MAX（A2：A13），計算出該組數據的極值。根據極值和實際需要確定分組間隔，本次分組間隔為2。在C2中輸入極小值整數位，在C3中輸入：＝C2＋2，選中C2：C3，按住鼠標左鍵將其拖拽至極大值即可得到分組。然后應用“直方圖”工具就能繪制頻率直方圖。首次使用分析工具需加載計算工具：文件－選項－自定義功能區在“開發工具”選項前打勾－確定，然后選擇開發工具－加載項－選中“分析工具庫”和“規劃求解”加載項－確定。工具加載完成后：數據－數據分析－直方圖，就能出現頻率直方圖對話框（圖3）。

圖2 三角圖解應用實例圖Fig．2 Application example of triangular chart

表2 頻率直方圖數據表Tab．2 Data table of frequency histogram

圖3 頻率直方圖對話框Fig．3 Dialog box of frequency histogram

圖3中“輸入區域”需要輸入原始數據：＝A2∶A13，在“接收區域”需要輸入分組數據：＝C2∶C9，輸出選項中選中“輸出區域”，輸入任意空白區域，本文實例中輸入：＝D1∶E9，選擇“圖表輸出”，最后確定，就能得到D列數據和相對應的頻率直方圖。如需要將多組數據在一張圖中表示，只要將每組數據經過上面的步驟處理，然后在某一頻率直方圖中添加數據即可：右擊－選擇數據－添加。該方法的關鍵是分析工具庫中直方圖的使用技巧。

3 概率格紙圖解法求異常下限

概率格紙圖解法確定背景值和異常下限，是建立在元素在地質體中呈正態分布或對數正態分布的基礎上。應用這種方法時，統計元素各個分組區間內的累積頻率，并在概率格紙上繪出各個累積頻率分布點的連線，然后根據其在概率格紙上反映的正態分布（或對數正態分布）特點，確定背景值及異常下限。其在Excel中的具體做法和步驟如下。

圖4是在Excel中繪制的概率格紙，其中縱坐標是正態分布累積頻率的反函數值（實際標記為累積頻率），橫坐標是元素含量的對數值（數據已通過正態檢驗）。表3中A、B、C和D列是繪制概率格紙橫向網格的數據，在A和C列中假設一組累積頻率，在B和D列中計算各自對應的反函數值，比如在B2中輸入函數：＝NORMSINV（A2%）。函數NORMSINV（probability）返回標準正態分布累積函數的反函數值。

圖4 元素含量累積頻率分布圖Fig．4 Chart of cumulative frequency distribution of element content

表3 概率格紙法求異常下限數據表Tab．3 Data table of probability ruling paper method to calculate anomaly threshold

以縱坐標為0．1%的橫向網格線為例，介紹橫向網格的繪制方法。首先繪制帶直線的散點圖：插入－散點圖－帶直線和數據標記的散點圖－在圖中右擊－選擇數據－添加，出現“編輯數據系列”對話框，在“系列名稱”中輸入：＝A4（累積頻率值0．1%），在“X軸系列值”中輸入：＝G17∶G18（圖4中X軸的極大值和極小值），在“Y軸系列值”中輸入：＝B4∶B5，點擊“確定”后出現一條橫向直線。選中該直線右擊，選“添加數據標簽”，在直線端點處出現標簽數據，選中該標簽數據右擊，選擇“設置數據標簽格式”并出現對話框，在“標簽包括”里只選擇“系列名稱”，在“標簽位置”中選擇“靠左”，關閉后刪除右端點的數據標簽即可。其他的網格線利用此方法一一添加，并將X軸的極大值和極小值固定為G17∶G18中的數值，刪掉Y軸坐標即可。

概率格紙繪制好之后，將需要計算的數據添加到概率格紙中，就可以計算異常下限和背景值。表3中E、F、G、H、I列中的數據是數據的處理過程，其中分組和頻數的計算參考上面頻率直方圖的方法。H列中累積頻率的計算以H2和H3為例分別輸入函數：＝G2／1644*100，＝G3／1644*100＋H2。I列中利用函數NORMSINV（probability）返回累積頻率的反函數值。所有數據計算好之后，將數據添加到已繪制好的概率格紙中，X軸輸入：＝F2∶F14（分組），Y軸輸入：＝I2∶I14（實際應為反函數值，但Y軸標記為累積頻率），適當調整后就能得到圖4。圖4中連線與累積頻率為50%線的交點的橫坐標為背景值，連線上累積頻率為97．7%的點橫坐標即為異常下限。該方法的關鍵在于概率格紙的繪制，但繪制好之后可以多次使用，以后只需將新的數據添加在已繪制好的格紙上即可。

4 多重分形法計算異常下限

目前利用分形技術進行地球化學異常下限確定的方法主要有：含量－周長法、含量－面積法、含量－距離法、含量－頻數法等。這里采用含量－頻數法，設分形求和模型：N（Ci）＝kC－Di（i＞0），式中Ci為元素含量，又稱特征尺度，k為比例常數（k＞0），D為一般分維數，N（Ci）為當元素含量為Ci時所有大于等于Ci的元素含量的和數。分形求和模型兩邊分別取對數得到一元線性回歸模型：lg N（Ci）＝－Dlg（Ci）＋lg（k），用最小二乘法求出斜率D的量，即為分維數。Excel中多重分形法的計算過程如下。

表4中對A列中的數據進行分組并計算頻數，C列對B列中的分組數據求對數，如在C2中輸入：＝LOG10（B2），F列中為當元素含量為Ci時所有大于等于Ci的元素含量的和數，如在F5中輸入函數：＝E5＋F6，G列對F列中的數據求對數，如在G2中輸入：＝LOG10（F2）。然后繪制散點圖：X軸中輸入：＝C2∶C15，Y軸中輸入：＝G2∶G15，其散點大致分布在兩段直線上，同時在圖中可確定兩段直線的分界點。根據分界點，兩段直線在表4中的坐標數據分別為：X＝C2∶C8，Y＝G2：G8和X＝C9：C15，Y＝G9∶G15。確定好兩段直線的坐標數據之后，重新繪制散點圖，分別輸入兩段散點的X、Y坐標，生成兩段散點，然后分別選中一段散點，右擊彈出下拉菜單，選擇“添加趨勢線”，彈出“設置趨勢線格式”對話框，在“趨勢線選項／回歸分析類型”中選擇“線性”，同時在“顯示公式”前的方框中打勾，最后關閉即可。該方法的關鍵在于理解多重分形法的原理。

表4 多重分形法數據表Tab．4 Data table of multifractal method

圖5 元素含量－頻數雙對數曲線Fig．5 Chart of element content frequency double logarithmic curve

5 R型聚類分析

R型聚類分析是根據樣品的多種變量的測定數據進行數字分類，定量地確定變量之間的親疏程度［5］。進行數字分類，需要選擇合適的數量指標，以此衡量樣品之間的親疏程度。數量指標主要有距離系數，相似系數和相關系數。本研究用相關系數作相似性統計量的逐步計算形成法做R型聚類分析，將通過一個實際例子，說明在Excel中聚類分析的計算步驟。

由于各變量的單位，量級和數值變動范圍的差異很大，計算中往往突出了那些絕對值較高的變量。因此，在進行聚類分析之前需將各個變量換算成一致的相對值。常用的變換的方法有標準化和正規化，本例中選擇標準化。例如表5中B2∶G7為原始數據，在第8和第9行中分別計算均值和標準差，比如分別在B8、B9中輸入函數：＝AVERAGE（B2：B7）和＝STDEV（B2：B7），然后選中拖動即可全部算出。在B11∶G16中計算出標準化數據，比如在B11中輸入函數：＝（B2－B8）／B9。標準化之后的數據，均值為0，標準差為1。數據標準化之后即可計算相關矩陣：數據－數據分析－選中“相關系數”－確定，然后彈出“相關系數”對話框（圖6），在“輸入區域”選中B10∶G16，在“分組方式”中選擇逐列，選中“標志位于第一行”，在“輸出區域”中選擇A19，確定之后會得到表5中A9∶G25的相關矩陣。從相關矩陣中可以看出Cu和Co相關系數最大，因此首先將Cu和Co連為一組。逐步加權平均，即新的CuCo＝（Cu＋Co）／2，計算修正數據，將得到的新數據CuCo替換Cu和Co，并與其他數據一起，計算新的相關矩陣，直到所有元素均已分組完成（表6）。最后在繪圖軟件中根據上面的計算結果，繪制出譜系圖即可（圖7）。

表5 聚類分析數據表Tab．5 Data table of R cluster analysis

表6 R型聚類分析連結順序表Tab．6 Link sequence table of R cluster analysis

圖6 相關矩陣對話框Fig．6 Dialog box of correlation matrix

6 判別分析

若有兩個母體，按照某種準則，把它們的P種特征（變量）組合成一個綜合指標（關系式），稱為判別函數，使兩類母體的區分率達到最好［6］。判別函數中最簡單的形式是線性函數。有p個變量n個樣品的兩個母體的判別函數可以寫成：

代表每個變量的方差和，當k≠f時Skf代表兩兩變量之間的協方差和。方差／協方差的定義為同一／不同變量之間的偏差乘積的平均數。但是根據Skf表達式，Skf應為方差或協方差的n倍（n為樣品數）的和。實際胡以鏗［7］在《地球化學中的多元分析》中的計算實例中Sif也是方差或協方差的n倍的和。

圖7 R型聚類分析譜系圖Fig．7 Pedigree chart of R cluster analysis

作者以文獻［7］中的實例來說明在Excel中判別分析的計算過程。表7中第3到第9行為原始數據，共7個樣品（n＝7），在第10行計算每個變量的均值，比如在B10中輸入：＝AVERAGE（B3：B9），在第11行中計算兩組樣品同一變量的均值差，比如在B11中輸入：＝B10－E10，在第12行中計算每個變量方差的7倍，比如在B12中輸入：＝VAR．P（B3：B9）*7，在第13行中計算兩組樣品中同一變量方差的加和，比如在B13中輸入：＝B12＋E12，其中B13，C13，D13分別代表S11、S22、S33，同理在B14、B15中分別計算協方差的7倍以及協方差和，分別輸入：＝COVARIANCE．P（B3：B9，C3：C9）*7和＝B14＋E14，其中B15、C15、D15分別代表S12、S13、S23。經過上面方差和協方差的計算得到如下方程組：

將方程組的系數在第16、17、18行中按順序排列，在E16、E17、E18中分別輸入函數：＝B16*B19＋C16*C19＋D16*D19，＝B17*B19＋C17*C19＋D17*D19，＝B18*B19＋C18*C19＋D18* D19，此時E16、E17、E18中的值均顯示為零。然后利用“規劃求解”解線性方程組：數據－規劃求解－彈出“規劃求解參數”對話框（圖8），在“設置目標”中輸入：＝E16，在目標值中輸入：0．5571，在“更改可變單元格”中輸入：＝B19∶D19，在“遵守約束”中點擊“添加”，在彈出對話框中分別輸入：E17＝F17，E 18＝F18，最后把“使無約束變量為非負數”前的“勾號”去掉，確定即可得到B19∶D19中的三個解，同時E16、E17、E18中的值也變為非零。獲得方程組的解之后，就可以列出判別函數，判別顯著性檢驗以及對未知樣品的判別，本文略。本文計算結果與文獻［7］中實例計算中的微有差別，這是在計算過程中保留有效數字個數不同造成的（本文計算過程中保留9位有效數字，實際在表7中顯示4位）。該方法的難點為判別分析的原理，方差／協方差的計算以及利用“規劃求解”解方程組技巧。

圖8 規劃求解對話框Fig．8 Dialog box of solution of programming

7 結論

Excel為大眾軟件，除本研究中提到的應用外，微量／稀土配分圖，散點圖，相關性計算，正態檢驗，樣品化驗數據誤差的檢驗，一次趨勢面分析等等都可以用Excel來實現。因此熟練使用Excel基本能滿足地球化學的數據處理需求。地球化學數據處理多種多樣，有時需借助好幾種軟件甚至收費的軟件才能實現各種計算需求，每一種軟件又需要花費一定的時間去掌握和熟練。熟練使用Excel不失為較為理想地選擇，對地質人員較為適用。

表7 判別分析數據表Tab．7 Data table of discriminant analysis

參考文獻：

［1］ 蔣敬業，程建萍，祁士華，等．應用地球化學［M］．武漢：中國地質大學出版社，2006．

［2］ H．E．霍克斯，J．S．韋布．礦產勘查的地球化學［M］．謝學錦，譯，廊坊：地質科學院物探研究所，1974．

［3］ 伍宗華，古平．隱伏礦床的地球化學勘查［M］．北京：地質出版社，2000．

［4］ 李明欣，梁斌，王全偉，等．川西龍泉山白堊系泥質巖的元素地球化學特征［J］．高校地質學報，2013，19（2）：346－354．

［5］ 春乃芽．利用Excel實現R型聚類分析［J］．物探與化探，2007，31（4）：374－376．

［6］ 春乃芽．利用Excel實現判別分析［J］．物探化探計算技術，2007，29（6）：560－564．

［7］ 胡以鏗．地球化學中的多元分析［M］．武漢：武漢地質學院地球化學教研室，1984．

［8］ 羅先熔，文美蘭，歐陽菲，等．勘查地球化學［M］．北京：冶金工業出版社，2007．

Advanced application of Excel in geochemical data processing

TAN Qin-ping1，2，XIA Yong1*，XIE Zhuo-jun1，2，YAN Jun1，2
（1．State Key Laboratory of Ore Deposit Geochemistry，Institute of Geochemistry，Chinese Academy of Sciences，Guiyang 550002，China；2．University of Chinese Academy of Sciences，Beijing 100049，China）

The Excel functions of chart drawing and data computing facilitate geochemical data processing．With studying in the principle of geochemical data processing this paper explains the calculation steps of the triangular diagram，the frequency histogram，the probability ruling paper graphical method for anomaly threshold，the multifractal method to calculate anomaly threshold，R cluster analysis，and discriminant analysis．These methods are entirely applicable to geochemical data processes．

triangular diagram；anomaly threshold；histogram；cluster analysis；discriminant analysis

P 632

A

10．3969／j．issn．1001-1749．2014．05．21

1001-1749（2014）05-0626-08

2014-03-07 改回日期：2014-08-19

國家重點基礎研究發展計劃（973計劃）（2014CB440905）；礦床地球化學國家重點實驗室＂十二五＂項目群課題（SKLODG－ZY125－01）

譚親平（1986-），男，博士，從事構造地球化學研究，E-mail：565310821＠qq．com

*通訊作者：夏勇（1960-），男，博士生導師，從事礦床地球化學研究，E-mail：xiayong＠vip．gyig．ac．cn