遙感異常與地球化學指標間的映射關系模型探索

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

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

0 引言

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

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

1 研究方法及遙感數據預處理

1.1 研究方法

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

具體過程見圖1。

圖1 技術路線圖Fig.1 Technology roadmap

1.2 遙感數據預處理

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

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

2 蝕變信息提取

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

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

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

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

同理，研究區內相關黏土礦物在LandSat 8數據的Band 2 處的反射率呈現較高的反射峰，在Band 5處可見較小的吸收谷，而黏土礦物在Band 6處的反射率較高，在Band 7處黏土礦物均出現強烈的吸收特征。因此，選取Band 2、Band 5、Band 6、Band 7這4個波段進行主成份分析，并提取研究區內的羥基異常。主成份分析后得到的特征向量值見表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的加信息。根據含羥基類礦物的波譜特征，在LandSat 8遙感影像數據中的Band 2、Band 6波段強吸收特征，在Band 5、Band 7具有較高的反射率，通常還有Band 5、Band 7的貢獻值系數與Band 2的貢獻值系數相反，Band 2和Band 6符號相同的特征。最后選定PC 4作為主分量代表羥基異常，PC 4圖像中的高亮度信息就是表征含羥基類礦物的蝕變異常信息。

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

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

3.2 蝕變信息分級處理

為了進一步直觀地觀測異常信息的圖像，利用Surfer 軟件，在對相關主分量圖像進行3×3的高斯低通濾波處理和數字拉伸的基礎上，結合研究區實際情況，對于不同類型的蝕變信息建立相應的定量化的標準，用來進行閾值分割提取。本文采用主分量分析門限法劃分蝕變異常等級，即利用x+Kσ來劃分異常等級。x為均值，σ為標準差，K值一般取1～3。根據研究區內特點，異常均分為三級，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 化探數據處理與分析

3.1 化探數據的處理和統計分析

3.1.1 成礦元素的聚類分析

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

3.1.2 成礦元素的因子分析

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

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

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

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

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

表5 旋轉后的成份得分系數矩陣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 地球化學因子關系提取

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

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

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

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

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

分析提取出的兩種因子的性質，觀察其旋轉成份矩陣可以知道：F1主要代表了Cu、Co、Ni、Cr元素，反映研究去的巖石元素共生組合，F1因子中的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) 金川銅鎳礦區元素弱地球化學F1(Cu-Co-NI)因子異常分布圖 (b) 金川銅鎳礦區元素弱地球化學F2(Pb-Zn-Cr)因子異常分布圖

3.3 遙感地球化學統計學映射模型的建立

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

ρ=F(x)

(1)

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

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

可以求出：

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

(2)

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

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

(3)

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

4 應用研究

4.1 金川礦區測區地質地球化學概況

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

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

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

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

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

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

5 結論

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

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

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

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

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