陳三明，馬云來，韋 龍，董 霖，宋世鵬

(桂林理工大學 地球科學學院，廣西 桂林 541006)

0 引言

前人在這方面的工作成果頗為豐富，代表性研究工作有陳濤[1]基于TM數(shù)據(jù)利用掩膜去除干擾因素、采用主成份分析、密度分割方法結(jié)合地質(zhì)信息提取研究區(qū)礦化異常信息；蘇一鳴[2]等采用去干擾異常主分量門限化技術(shù)，提取羥基、鐵染異常；王宏[3]等利用高分一號、ETM+影像數(shù)據(jù)對區(qū)內(nèi)巖性、構(gòu)造進行解譯，通過主成份分析法和比值法結(jié)合ETM+和ASTER影像數(shù)據(jù)對區(qū)內(nèi)礦化蝕變信息進行提??；唐堯[4]等利用遙感數(shù)據(jù)與蝕變礦物的譜帶特征提取礦化蝕變信息，對ETM的不同波段組合進行增強處理后確定圍巖蝕變的范圍和強度；汪子義[5]等利用混合調(diào)制匹配濾波方法開展基于OLI數(shù)據(jù)的蝕變礦物組合提取研究，并與Hyperion 高光譜數(shù)據(jù)提取結(jié)果進行對比驗證，分析遠景區(qū)的蝕變礦物組合分帶特征；賀金鑫[6]等基于LandSat 8數(shù)據(jù)，以主成份分析法和比值法結(jié)合的方法增強蝕變信息，等等。但由于遙感信息的不確定性和反演多解性，單純的利用遙感礦化蝕變信息找礦實際效果并不太理想。

筆者認為，如果將宏觀層面上的遙感礦化蝕變信息與微觀層面具有明確找礦指示意義的地球化學示礦指標在中觀層面上建立模型聯(lián)系起來，通過對比印證，就能夠大大減少單一運用遙感信息帶來的不確定性和反演的多解性，從而就能增加遙感找礦的可靠性。

1 研究方法及遙感數(shù)據(jù)預處理

1.1 研究方法

選取LandSat 8 數(shù)據(jù)、6種元素(Cu、Pb、Zn、Co、Ni和Cr)的化探數(shù)據(jù)為基本的信息源，先對遙感影像的多光譜波段(Band 1～Band 7)進行預處理，包括輻射校正(包括輻射定標、大氣校正)、坐標精校正、高程校正，再和配準后的全色波段(Band 8)融合，然后對融合后的影像利用已知的工作區(qū)數(shù)據(jù)建立MASK掩膜，并選取相應(yīng)的波段，使用PCA主成份分析法、比值法結(jié)合閾值分割法得到增強的蝕變異常信息，初步劃定異常區(qū)域。最后把得到的增強后的弱化蝕變信息同化探數(shù)據(jù)分析出的弱化異常信息對比分析，在研究區(qū)地質(zhì)資料及氣象水文等資料相結(jié)合的基礎(chǔ)上，通過建立地球化學指標與LandSat 8的特征波段凈反射率之間的統(tǒng)計學映射關(guān)系，建立示礦綜合地球化學指標與LandSat 8特征波段凈反射率組合關(guān)系定量模型，以此為基礎(chǔ)，對未知區(qū)域開展遙感地球化學找礦預測，最終圈定異常區(qū)域和預測區(qū)域。

具體過程見圖1。

圖1 技術(shù)路線圖Fig.1 Technology roadmap

1.2 遙感數(shù)據(jù)預處理

這里的數(shù)據(jù)預處理指的是對遙感數(shù)據(jù)LandSat 8的處理，化探數(shù)據(jù)的處理會在下面的章節(jié)中另外闡述。數(shù)據(jù)預處理包括遙感原始數(shù)據(jù)的輻射定標、大氣校正、坐標精校正、高程校正和圖像融合。輻射定標是將傳感器記錄的電壓或數(shù)字量化值轉(zhuǎn)化為輻射亮度(Radiance)或者轉(zhuǎn)化為地表(表觀)反射率(Reflectance)和亮度溫度(Brightness Temperature)的過程，方法有絕對定標和相對定標兩種。本次工作直接利用ENVI 5.3的輻射定標工具(Radiometric Calibration)來完成對原始數(shù)據(jù)的定標操作。大氣校正的目的是消除大氣和光照等因素對地物反射的影響，獲取地物反射率、輻射、地表溫度等真實物理模型參數(shù)，用來消除大氣中水蒸氣、氧氣、二氧化碳、甲烷和臭氧等對地物反射的影響，消除大氣分子和氣溶膠散射的影響。大多數(shù)情況下，大氣校正同時也是反演地物真實反射率的過程，方法包括基于大氣參數(shù)的模型校正方法以及基于遙感數(shù)據(jù)統(tǒng)計特征規(guī)律的回歸分析方法。坐標精校正是賦予遙感數(shù)據(jù)像元的空間定位屬性，根據(jù)研究區(qū)的分帶本次工作統(tǒng)一使用西安80坐標系。高程校正是為了減少地面高程差對提取結(jié)果的影響，在高差地貌其伏較大的山地丘陵地區(qū)尤為明顯。圖像融合是指使用特定的數(shù)學算法，將兩幅或多幅遙感影像數(shù)據(jù)，在設(shè)定好的地理坐標系中融合合成一幅新的影像的過程。本次工作采用多光譜影像數(shù)據(jù)(Band 1～Band 7)和全色波段(Band 8)融合，既保留了豐富的光譜信息，同時提高了的空間分辨率，為接下來的信息提取提高了精度。圖2為預處理前后遙感影像的對比圖。

圖2 數(shù)據(jù)預處理前后遙感圖像對比圖Fig.2 Comparison of remote sensing images before and after data preprocessing

2 蝕變信息提取

2.1 主成份分析法與比值法相結(jié)合的蝕變信息提取

2.1.1 鐵染和羥基異常信息提取

研究區(qū)內(nèi)與鐵離子有關(guān)的礦物主要是赤鐵礦、針鐵礦等，根據(jù)以上礦物的波譜特征，在0.45～0.52 μm和0.875 μm附近的波曲線具有2個吸收谷，其正好對應(yīng)LandSat 8數(shù)據(jù)的Band 2和Band 5波段。0.630～0.680 μm和1.500～1.700 μm附近可見較高的反射峰，分別對應(yīng)LandSat 8數(shù)據(jù)的Band 4和Band 6波段[6]。因此可以選取Band 2、Band 4、Band 5、Band 6波段進行主成份分析提取鐵染異常。但是，經(jīng)分析后發(fā)現(xiàn)這4個波段所提提取的特征向量值并不符合鐵離子礦物的波譜特征。根據(jù)梁昊[7]、馬威[8]等人的資料進一步研究發(fā)現(xiàn)，利用比值法即Band 6/Band 5，可以進一步增強鐵化蝕變圖像的亮度反差。由此可見，選用band 6/Band 5代替Band 6，即用Band 2、Band 4、Band 5、Band 6/Band 5波段組合進行主成份分析，能夠進一步增強鐵染蝕變信息。經(jīng)過主成份分析后計算所得特征向量值見表1。

負值代表吸收，正值代表反射，絕對值的大小代表能力的強弱。

由表1可見，PC 1主要反映的是Band 4、5的加信息；PC 2主要反映的是Band 2的加信息和Band 5的減信息；PC 3主要反映的是Band 3的加信息和Band 2的減信息；PC 4則幾乎全部反映Band 6/Band 5的加信息。

表1 鐵染蝕變異常特征向量Table 1 Feature vector of iron-stained alteration anomaly

根據(jù)鐵染蝕變異常的提取原理，波段特征信息在LandSat 8遙感影像數(shù)據(jù)中的Band 2、Band 5波段具有強吸收特征與Band 4或Band 6/Band 5具有較高的反射率，通常還具有Band 2、5的貢獻值系數(shù)與Band 4的貢獻值系數(shù)相反，Band 4和Band 6符號相同的特點。因此選擇PC 3作為包含鐵染蝕變異常的主分量進行分析。

同理，研究區(qū)內(nèi)相關(guān)黏土礦物在LandSat 8數(shù)據(jù)的Band 2 處的反射率呈現(xiàn)較高的反射峰，在Band 5處可見較小的吸收谷，而黏土礦物在Band 6處的反射率較高，在Band 7處黏土礦物均出現(xiàn)強烈的吸收特征。因此，選取Band 2、Band 5、Band 6、Band 7這4個波段進行主成份分析，并提取研究區(qū)內(nèi)的羥基異常。主成份分析后得到的特征向量值見表2。

表2 羥基異常特征向量Table 2 Characteristic vector of hydroxyl anomaly

由表2可見，PC 1主要反映的是Band 5、Band 6、Band 7的的加信息；PC 2主要反映的是Band 2、Band 5的加信息和Band 7的減信息；PC 3主要反映的是Band 2的加信息；PC 4主要反映的是Band 6的減信息和Band 7的加信息。根據(jù)含羥基類礦物的波譜特征，在LandSat 8遙感影像數(shù)據(jù)中的Band 2、Band 6波段強吸收特征，在Band 5、Band 7具有較高的反射率，通常還有Band 5、Band 7的貢獻值系數(shù)與Band 2的貢獻值系數(shù)相反，Band 2和Band 6符號相同的特征。最后選定PC 4作為主分量代表羥基異常，PC 4圖像中的高亮度信息就是表征含羥基類礦物的蝕變異常信息。

3.1.2 交互式閾值拉伸及低通濾波處理

通過分析特征向量后，利用交互式拉伸工具增強弱化的鐵染、羥基蝕變信息，能更直觀的觀測到異常。通常取統(tǒng)計數(shù)據(jù)的均值加二至三倍的均方差數(shù)值拉伸異常，然而這樣得到的遙感異常出現(xiàn)所謂的“焦鹽現(xiàn)象”，這不利于異常的空間性規(guī)律識別。此時再將異常的數(shù)據(jù)在ENVI中的數(shù)學形態(tài)學模塊里，進行低通濾波處理，得到模糊影像，但遙感異常的空間性得到了顯著的增強，一些碎點異常得到了有效的抑制。

3.2 蝕變信息分級處理

為了進一步直觀地觀測異常信息的圖像，利用Surfer 軟件，在對相關(guān)主分量圖像進行3×3的高斯低通濾波處理和數(shù)字拉伸的基礎(chǔ)上，結(jié)合研究區(qū)實際情況，對于不同類型的蝕變信息建立相應(yīng)的定量化的標準，用來進行閾值分割提取。本文采用主分量分析門限法劃分蝕變異常等級，即利用x+Kσ來劃分異常等級。x為均值，σ為標準差，K值一般取1～3。根據(jù)研究區(qū)內(nèi)特點，異常均分為三級，K值分別為1、1.5、2，即以σ、1.5σ、2σ為閾值劃分為一級、二級和三級異常。圖3、圖4分別為加載底圖后的鐵染蝕變異常和羥基異常。

圖3 加載底圖后的鐵染蝕變異常分布圖Fig.3 Distribution of iron-stained alteration anomaly after loading of the base map

3 化探數(shù)據(jù)處理與分析

3.1 化探數(shù)據(jù)的處理和統(tǒng)計分析

3.1.1 成礦元素的聚類分析

在實驗室內(nèi)先對野外樣品進行分析測試，然后利用SPSS對地球化學數(shù)據(jù)進行統(tǒng)計分析，根據(jù)分析結(jié)果，最終選定Co、Ni、Cr、Cu、Pb、Zn六種元素作為數(shù)據(jù)樣本。通常為了描述數(shù)據(jù)衡量各組數(shù)據(jù)之間的近似源關(guān)系，通過不同的計算方法將需要研究的數(shù)據(jù)進行分類，將多個類別數(shù)據(jù)歸結(jié)到少數(shù)幾個大類的不同簇中，一般使用聚類分析的方法。地球化學元素本身在地表漫長的演化過程中，經(jīng)過各種地質(zhì)作用以及地球化學作用下形成如今的地球化學環(huán)境，其組合具有一定的解釋成礦規(guī)律的能力。處理時，先把六種元素的地球化學數(shù)值做甚高值去除后標準化，然后聚類分析，得到樹狀圖。根據(jù)六元素聚類分析樹狀聯(lián)接圖可以看出，這些元素之間存在的關(guān)聯(lián)性可知大體可分為兩大類。其中Pb、Co聚為一類，Cu、Ni、Zn、Cr聚為另一類。聚類分析的結(jié)果不是特別的理想，表明我們需要進行進一步的地球化學指標的提取和組合。

3.1.2 成礦元素的因子分析

因子分析是較常用的一種統(tǒng)計學分析方法，一般伴隨使用聚類分析作為參考。元素地球化學數(shù)據(jù)因子分析結(jié)果見表3。

從相關(guān)性矩陣表3可以求得其特征值和累積百分比，見表4。從主因子的累積百分比看到，前兩個特征值所代表的方差已經(jīng)占總方差的83%。

根據(jù)元素統(tǒng)計情況及因子分析情況結(jié)合碎石圖，可知一共有3個特征值超過1的主成份，這3個成份的貢獻率按順序分別是72.838%，10.470%，7.833%，可以看出最后因子的貢獻率較低不足10%，因此舍棄只保留前2個主成份因子。由表4得知，前2個主成份因子的累積貢獻率達到了83.308%，即代表這2個主成份因子可以解釋表達了原數(shù)據(jù)的83.308%的主要內(nèi)容。

表4 總方差解釋Table 4 Explanation of total variance

結(jié)合因子分析旋轉(zhuǎn)成份矩陣表(表5)分析，兩個因子F1、F2已經(jīng)提取了Cu、Pb、Zn、Co、Ni、Cr，元素的大部分信息。根據(jù)旋轉(zhuǎn)成分矩陣，計算得出2個因子公式如下：

表5 旋轉(zhuǎn)后的成份得分系數(shù)矩陣Table 5 Component score coefficient matrix after rotation

F1=0.393×Cu-0.316×Pb+0.023×Zn+0.024×Co+0.485×Ni+0.038×Cr

F2=-0.323×Cu+0.573×Pb+0.330×Zn-0.198×Co-0.262×Ni+0.684×Cr

3.2 地球化學因子關(guān)系提取

指示元素的組合關(guān)系一般包括定性組合和定量(半定量)組合[9]。定量(半定量)的組合關(guān)系主要用來表征各個元素間數(shù)學特征值及其空間分布規(guī)律狀況。按組合性質(zhì)分為如下幾種表現(xiàn)形式：

1)簡單含量關(guān)系：主要根據(jù)相關(guān)性質(zhì)，如劃分限定的閾值、利用某些指示元素之間的含量特征形式來指示含礦異常和非礦異常。

2)比值：指單元素之間的比值關(guān)系，因其在礦床地球化學異常三度空間中隨空間部位動態(tài)變化形成關(guān)系式效應(yīng)，故比值的關(guān)系映射動態(tài)變化具有一定的規(guī)律性。

3)多元素組合比值：分為壘乘異常及比值和壘加異常及比值，例如：(a×b)/(c×d)和(a+b)/(c+d)。這種組合關(guān)系一般是相對的，沒有比較嚴格的組合標準?？梢詤⒖荚氐幕顒有皂樞騺韯澐纸M合關(guān)系。

4)指示元素之間含量的關(guān)系：主要表現(xiàn)在元素對或多元素的組合含量之間的共消長關(guān)系和互消長關(guān)系。利用指示元素含量間相關(guān)性來確定地球化學異常時，將其中某些元素的相關(guān)系數(shù)用作指標。

分析提取出的兩種因子的性質(zhì)，觀察其旋轉(zhuǎn)成份矩陣可以知道：F1主要代表了Cu、Co、Ni、Cr元素，反映研究去的巖石元素共生組合，F(xiàn)1因子中的Cu、Co、Cr、Zn、Pb在巖漿型銅鎳硫化物礦床中對成礦有明顯指示作用[7]，且Ni本身為良好指示元素。F2因子主要提取了Pb、Zn、Cr元素的信息，代表著成礦需要的重要元素組合。通過多次的組合計算分析，最終選定地球化學指標因子F1(Cu-Co-Ni)與因子F2(Pb-Zn-Cr)，把F1/F2的值作為找礦地球化學敏感指標，這樣就建立了遙感地球化學指標間的模型的因變量(圖5)。

圖5 F1、F2因子異常分布Fig.5 Anomaly distribution of F1 and F2 factors(a) 金川銅鎳礦區(qū)元素弱地球化學F1(Cu-Co-NI)因子異常分布圖 (b) 金川銅鎳礦區(qū)元素弱地球化學F2(Pb-Zn-Cr)因子異常分布圖

3.3 遙感地球化學統(tǒng)計學映射模型的建立

統(tǒng)計分析地球化學元素時，在樣本微量元素與光譜反射率統(tǒng)計分析的元素選擇中，要先對微量元素的結(jié)果進行多元統(tǒng)計，分析其與主成礦元素的相關(guān)性；然后篩選出研究區(qū)中和主成礦元素成礦關(guān)系密切的元素組合，對選出的元素組合進行標準化，提取出最能代表主成礦元素成礦的元素綜合異常數(shù)據(jù)(用x表示)，建立地物光譜反射率與微量元素地球化學指標之間的關(guān)系函數(shù)。對上述篩選出的樣本光譜反射率、元素綜合異常數(shù)據(jù)進行統(tǒng)計分析，得出其相關(guān)函數(shù)，見公式(1)：

ρ=F(x)

(1)

其中：ρ為樣本光譜反射率；F為樣本光譜反射率與微量元素地球化學指標之間的統(tǒng)計關(guān)系(函數(shù))?？梢郧蟪鲞b感影像數(shù)據(jù)灰度值與地物光譜反射率之間的方程為

L=ρ(E0·cosθ)/π·d2=DN·gain+bias，

可以求出：

ρ=(DN·gain+bias)·π·d2/E0·cosθ

(2)

其中：L為地物在大氣頂部的輻射能量值；ρ為樣品的光譜反射率；DN為樣本的灰度均值；d為日地天文單位距離；E0為大氣頂部的太陽輻照度；θ為成像時的太陽天頂角；gain和bias分別為圖像的增益量與偏移量，可從遙感數(shù)據(jù)的頭文件中讀取。

由公式(1)、(2)組成方程組，可以推導出DN與X之間的對應(yīng)關(guān)系：

(3)

由此通過F(x)函數(shù)，建立了遙感影像灰度值與地球化學異常之間的判別函數(shù)關(guān)系，也即提取出的遙感地球化學異常模型。最后將遙感數(shù)據(jù)各波段的灰度值代入式(3)，就可以得到代表綜合元素異常值的一組數(shù)據(jù)，其值的高低代表所求遙感地球化學異常的強弱，x的高值區(qū)就是通過遙感數(shù)據(jù)獲取的遙感地球化學異常區(qū)。

4 應(yīng)用研究

4.1 金川礦區(qū)測區(qū)地質(zhì)地球化學概況

金川礦床位于甘肅西部龍首山脈東段北坡，是世界第三大銅鎳硫化物礦床。大地構(gòu)造位置上位于華北地臺阿拉善地塊西南緣的龍首山隆起帶內(nèi)，北依準噶爾晚古生代褶皺系，南鄰北祁連褶皺帶，東部為華北地臺，西隔祁連早古生代褶皺系與柴達木地塊相望。

金川礦區(qū)南延測區(qū)出露從太古代到古生代的一系列地層，老地層的原始沉積建造已被徹底改造。受區(qū)域性構(gòu)造作用影響，各巖性帶呈近NW-SE向帶狀平行展布。測區(qū)出現(xiàn)多期強烈的基性—超基性巖漿活動為特征，孕育著超大型銅鎳礦床的白家嘴子超基性雜巖體位于普查區(qū)北面，輝綠巖體主要呈脈狀或巖枝產(chǎn)于各時代的地層中。

測區(qū)地球化學統(tǒng)計表明，Cu、Pb、Zn等金屬元素在各時代蝕變火山巖中含量最高，在變質(zhì)碎屑巖及含碳變質(zhì)巖中較高，在泥質(zhì)及碳酸鹽巖中含量普遍較低。區(qū)內(nèi)Cu、Pb、Zn等金屬元素的富集成礦與火山活動及區(qū)域變質(zhì)作用密切相關(guān)。印支期侵入巖Pb元素含量普遍較高，為地殼克拉克值的2倍以上；海西晚期鉀長花崗巖Pb含量均較高，為地殼克拉克值的2倍；二長花崗巖Pb、Zn、As、Ni、Co、Mo含量均較高，超過了地殼克拉克值。海西中期各類巖體的Pb元素含量普遍偏高，其他元素含量則相對偏低。超基性巖Ni、Co元素含量是克拉克值的200倍以上，其他元素含量則較低。加里東期斜長花崗巖中的Pb、As含量普遍大于地殼克拉克值，Cu、Zn只在局部地段含量較高：石英閃長巖As、Ni、Co含量都高于地殼克拉克值數(shù)倍以上。不同地區(qū)不同期次的中性、酸性侵入巖中，Cu、Pb、Zn元素含量普遍較高，并呈正消長關(guān)系[10]。

4.2 遙感地球化學異常指標的提取結(jié)果分析及異常評價

根據(jù)LandSat 8遙感地球化學異常綜合模型進行信息提取的空間展布特征，依據(jù)異常大小進行密度分割，本區(qū)可劃分為2個成礦遠景預測區(qū)，分別編號為1、2號(圖6)。第1個預測大區(qū)(一級預測區(qū))之中分2個小預測區(qū)，分別標記為1-1、1-2、，第2個預測大區(qū)(二級預測區(qū))只有1個別標記為2-1。其中，第1預測區(qū)1-1為已開采的礦床，從預測的結(jié)果來看說明本文建立的基于灰色關(guān)聯(lián)度的高光譜遙感地球化學綜合預測模型達到了較為理想的結(jié)果。1-2預測區(qū)則與NW向的基性超基性巖脈的空間分布較為一致，且形成NW向的條狀軸向分布格局。地層巖性主要為上元古界震旦系燒火筒群青灰色微晶白云巖、含礫絹云千枚巖等。構(gòu)造以斷裂為主，F(xiàn)2、F3和F4斷裂總體走向NW，異常位于3組斷裂圍限的震旦系中。巖漿巖主要為基性的輝綠巖脈，少量超基性輝橄巖脈，特點是巖脈出露面積較大，形態(tài)與走向與全區(qū)一致，輝綠巖脈一般長1500～1800 m，寬30～50 m，其余小巖脈呈零星分布。第2預測區(qū)2-1為團狀分布格局，似近SN向展布，與礦區(qū)的整體構(gòu)造走向近乎垂直，需要進一步綜合其他勘查手段作進一步的分析研究。

圖6 金川銅鎳礦區(qū)遙感地球化學異常預測Fig.6 Remote sensing geochemical anomaly prediction in Jinchuan copper-nickel mining area

5 結(jié)論

1)從本次實驗的結(jié)果來看，如果僅從遙感異常上提取信息容易受不確定性因素的影響；如果僅有化探數(shù)據(jù)異常的提取，也難于發(fā)現(xiàn)異常的空間展布規(guī)律。所以融合遙感異常及地球化學指標異常兩種方法，可以更為準確、高效地提取礦致異常信息。

2)對特定的地質(zhì)研究體，以LandSat 8 OLI的多光譜及全色波段的遙感影像特征，通過統(tǒng)計分析，建立地球化學成礦因子與遙感影像閾值分級相結(jié)合的遙感地球化學指標模型，利用遙感地球化學指標勾繪出遙感地球化學異常圖件，可以對未知區(qū)進行初步的預測。

3)遙感地球化學模型的建立必須依托于實際的化探工作與已知礦點，因此不同的地質(zhì)地球化學背景的預測應(yīng)有不同的模型，并且預測模型應(yīng)有一定的外延范疇，它在同一巖性區(qū)域有較好的預測效果。但遇到構(gòu)造較強烈的地段，表現(xiàn)出的規(guī)律性不強，建議對遙感構(gòu)造地球化學進行專門的研究。

4)基于反射波譜建立的遙感地球化學預測模型的機理還需要進行深入的研究。

5)遙感地球化學綜合異常提取的結(jié)果一定要做后期的異常評價，需要與實際的數(shù)據(jù)對比驗證，結(jié)合地質(zhì)、地球物理勘探資料才能更好的對未知區(qū)進行科學地預測。