鄭 緯,朱谷昌,,張遠飛,石菲菲,王冬寅,李 紅

(1.中南大學(xué)地學(xué)與環(huán)境工程學(xué)院,長沙410083;2.有色金屬礦產(chǎn)地質(zhì)調(diào)查中心,北京100012)

0 引言

礦化蝕變信息是地質(zhì)演化過程中的地質(zhì)異常事件,是礦床在生成、演化全過程中的成礦條件表現(xiàn)[1]。由于礦床與地質(zhì)體的各種異常屬性相關(guān)聯(lián),同時這些異常屬性又與相對均一的地質(zhì)體背景屬性所有偏離和差別,礦化蝕變信息就是包含在地質(zhì)體綜合信息中的各種屬性異常微弱信息,人們必須采用各種可能的方法去捕捉和分辨地質(zhì)體各種屬性的異常,才能更好地預(yù)測礦床。遙感地質(zhì)方法正是通過識別遙感信息中的礦化蝕變信息進行成礦預(yù)測的一種技術(shù)手段。

多層次分離提取技術(shù)就是以巖石和礦物的電磁波特征反射譜帶作為提取巖石蝕變帶信息的理論基礎(chǔ),通過數(shù)學(xué)建模和最佳變量的選擇,運用各種數(shù)學(xué)方法和圖像增強手段,多層次地從遙感信息中逐步剔除背景信息(即干擾信息),一次次地分離提取出礦化蝕變信息(即目標特征信息),最后得到包括目標特征信息的圖像[2]。

本文在“遙感蝕變信息多層次分離提取”的主技術(shù)框架下,通過遙感圖像特征診斷,把復(fù)雜的遙感蝕變信息提取轉(zhuǎn)化為遙感圖像的背景、干擾與蝕變信息3個主要對象的分析;然后,應(yīng)用光譜數(shù)據(jù)空間(即特征空間)的幾何結(jié)構(gòu)分析對背景、干擾與蝕變信息3個主要對象的聚類結(jié)構(gòu)進行研究與空間定位,特別是對干擾因素進行識別與劃分[3-5],在分析的基礎(chǔ)上再進行菲律賓宿務(wù)島地區(qū)礦化蝕變信息的檢測與提取。

1 宿務(wù)島礦區(qū)地質(zhì)背景與數(shù)據(jù)庫

工作區(qū)位于位于菲律賓宿務(wù)島北部宿務(wù)省宿務(wù)市 Kansi村,距宿務(wù)市北約 80 km,面積約16.92 km2。礦區(qū)地形為低山區(qū),炎熱潮濕,植被為雜草和樹木混合。

菲律賓群島屬于西太平洋新生代島弧火山巖帶,區(qū)內(nèi)地層以新生界為主,廣泛分布第三系中性火山巖,其基底為前侏羅系,分布于菲律賓群島中部和西部局部地區(qū)。

菲律賓島弧區(qū)火山活動早期以古新世及漸新世為主,主要分布在東部大斷裂帶的東緣;晚期以中新世及上新世為主,主要分布在大斷裂的兩側(cè)。侵入巖以閃長巖為主,還有花崗閃長巖、石英二長巖和正 長巖等。可分為2期:早期為35～60 Ma(漸新世),與塊狀硫化物和斑巖銅(鉬)礦有關(guān);晚期為1.5～2.5 Ma,與富金的斑巖銅礦有關(guān)。火山熱液型金銀礦與晚期中(基)性火山巖有關(guān)。超基性巖主要分布于大斷裂西側(cè),東側(cè)較少。

本文選用了1景 TM數(shù)據(jù),即113/53(時相為1992-06-29)。以 TM為代表的Landsat類多光譜遙感數(shù)據(jù),其地面覆蓋范圍寬,空間分辨率、光譜分辨率能滿足區(qū)域地質(zhì)調(diào)查及地質(zhì)信息提取的要求,數(shù)據(jù)形式易于增強處理,多波段優(yōu)化組合的圖像信息豐富,是中小比例尺區(qū)域遙感地質(zhì)調(diào)查的理想數(shù)據(jù)源。

2 礦化蝕變圖像特征分析與提取

2.1 礦化蝕變的多波段數(shù)據(jù)背景-異常子空間分析

根據(jù)多層次分離提取技術(shù)思路,首先對遙感圖像的背景-異常子空間進行劃分[6]。通過遙感多波段數(shù)據(jù)的背景端元數(shù)M值的估計方法,能夠基本上確定出背景-異常子空間。背景包括巖石、土壤等地物,同時也包含大的干擾(如植被區(qū)、云區(qū)、水系等)。小的干擾是通過光譜數(shù)據(jù)的空間結(jié)構(gòu)分析來識別。

背景端元數(shù)M值的估計是通過赤池信息準則(AIC準則)或最小描述長度準則(MDL準則)來實現(xiàn)的。

本文根據(jù) Crosta法先選擇 F(鐵化)因子的TM1,TM3,TM4,TM5等4個波段組合計算,得到標準差與方差分數(shù)(即貢獻率)(表1)、M值的估計函數(shù) Q(M)值(表 2)。

分析表1的方差分數(shù),前2個主成分的貢獻率已達到95.711%,表明主要背景(含大干擾)地物為前2個主成分;第3和第4主成分的合計貢獻率僅為4.29%,主要為一些小干擾與可能的鐵化蝕變異常信息。主成分4的貢獻率比較低,結(jié)合Q(M)值可以預(yù)計這里包含的除了噪聲之外,還有小干擾或弱信息。

M是主成分的本征數(shù),背景(含大干擾)的端元數(shù)目應(yīng)該等于M+1。分析表2中Q(M)值的變化會發(fā)現(xiàn),這些值具有2個臺階,第一個臺階值范圍是0.845～0.853,表明主要存在3類背景地物(含大干擾);第二個臺階值范圍是0.853～1.000,表明小干擾與可能的蝕變異常信息屬于臺階變化。這同前面的方差貢獻率分析是一致的。

表1 Crosta法TM1,TM3,TM4,TM5的標準差與貢獻率Table 1 Standard division and contribution rate of TM1,TM3,TM4 and TM5 calculated with Crosta method

表2 M值的估計函數(shù) Q(M)值Table 2 The estimated value of function Q(M)of M value

表3 Crosta法TM1,TM4,TM5,TM7的標準差與貢獻率Table 3 Standard division and contribution rate of TM1,TM4,TM5 and TM7 calculated with Crosta method

選擇 H(泥化)因子 TM1,TM4,TM5,TM7等4個波段組合計算,得到標準差與方差分數(shù)(即貢獻率)(表3)、M值的估計函數(shù) Q(M)值(表4)。

分析表3的方差分數(shù),前2個主成分的貢獻率已達到95.188%,表明主要背景(含大干擾)地物為前2個主成分;第3和第4主成分的合計貢獻率僅為4.812%,這2個主成分含有一些小干擾與可能的泥化蝕變異常信息。主成分4的貢獻率比較低,結(jié)合Q(M)值可以預(yù)計這里包含的除了噪聲之外,還有小干擾或弱信息。

表4 M值的估計函數(shù) Q(M)值Table 4 The estimated value of function Q(M)of M value

M是主成分的本征數(shù),背景(含大干擾)的端元數(shù)目應(yīng)該等于M+1。分析表4中Q(M)值的變化會發(fā)現(xiàn),這些值的變化具有2個臺階,第一個臺階值范圍是0.711～0.752,表明主要存在3類背景地物(含大干擾);第二個臺階值范圍是0.752～1.000,表明小干擾與可能的蝕變異常信息屬于這個臺階。這同前面的方差貢獻率分析是一致的。

結(jié)合遙感圖像分析,發(fā)現(xiàn)該地區(qū)植被覆蓋面很廣、干擾強烈,而且已構(gòu)成主背景地物,水體面積也較大,而其他巖石裸露區(qū)則是小背景地物。很明顯,一些小干擾與可能的蝕變異常信息會“湮沒”在這2類背景地物之中。

2.2 基于背景異常分析的礦化蝕變異常的特征波段選擇

在對遙感圖像數(shù)據(jù)進行背景-異常子空間分析后,通過回歸偏度曲線分析和地物光譜反射特征來選擇信息提取的特征波段。

回歸偏度曲線是用來描述線性回歸的二維散點圖的圖形不對稱性,回歸偏度較大,則預(yù)示著有不同于背景的異常或干擾信息存在[7]。

圖1是由剔除了水體影響后的7個波段的回歸偏度曲線。圖中顯示出 TM4,TM5,TM3,TM6和TM1等5個波段的偏度較大,說明它們載荷了較多的背景、干擾或蝕變異常信息。由地物波譜的物理機制可知,偏度最大的 TM4反映的是植被信息;第2位的 TM5既有植被信息、巖石信息,還可能載荷有泥化蝕變信息;第3位的TM3包含了鐵化蝕變信息;與 TM5形成反向偏度的 TM7應(yīng)該是由泥化蝕變信息引起的;同時注意到 TM3與 TM1的偏度反差不大,TM3與 TM4的偏度反差主要是由植被造成的,即 TM4高、TM3低。

由上述的分析得知,TM4,TM5,TM6,TM3和TM1是該地區(qū)重要的特征波段。所以波段選擇方案為:①剔除植被干擾的波段:TM4和 TM5;②巖石與裸露區(qū)識別的波段:TM5;③鐵化蝕變異常信息識別與提取波段:(偏度曲線反映不明顯,待光譜空間結(jié)構(gòu)分析后再定);④泥化蝕變異常信息識別與提取:TM5和 TM7。

2.3 礦化蝕變的波譜空間幾何結(jié)構(gòu)分析

在選定特征波段后,對遙感數(shù)據(jù)在波段組合中的光譜空間幾何結(jié)構(gòu)進行分析,從而確定不同蝕變信息提取的波段選擇。

遙感的多(高)光譜數(shù)據(jù)是一種多元數(shù)據(jù)集合,每1個像元代表1個波譜矢量。所以,圖像多元數(shù)據(jù)集合在高維空間中形成1個點陣,這個空間點陣具有一定的幾何結(jié)構(gòu)。

多(高)光譜數(shù)據(jù)集合屬高維空間中的點陣,高維空間的點陣幾何體形態(tài)是無法直接觀測的。但是,從低維空間入手,通過分析二維(2個波段)、三維(3個波段)數(shù)據(jù)的散點圖,可以探討高維情況下的數(shù)據(jù)點陣的結(jié)構(gòu)及在空間中的聚類分布特征。

以圖2為例來說明光譜數(shù)據(jù)的二維散點圖所包含的豐富信息:

(1)二維散點圖又稱二維直方圖,它是二維變量(波段)聯(lián)合概率密度分布的幾何表達,反映了二維變量數(shù)據(jù)空間的聚類結(jié)構(gòu)。

圖1 TM數(shù)據(jù)7個波段的回歸偏度曲線Fig.1 Partial regression curve of TM data of the 7 bands

圖2 光譜數(shù)據(jù)的二維散點圖Fig.2 Planar scatter plot of spectral data

(2)二維散點圖也是2個變量主成分分析的幾何解釋,其中 PC1為第1主成分軸,PC2為第2主成分軸。對于多波段(3個波段以上)主成分分析,比如Crosta定向主成分分析法,可以通過多個二維散點圖的觀測去重構(gòu)腦海中的高維光譜數(shù)據(jù)點陣幾何結(jié)構(gòu)。

(3)若約定圖形坐標系統(tǒng)的Bx(橫坐標)表示具有吸收峰的波段,比如 TM1或 TM7,而By(縱坐標)表示具有反射峰的波段,如 TM3或 TM5。根據(jù)幾何意義,坐標平面內(nèi)的點P(bx,by),位于對角線(圖2中綠色線)上方的數(shù)據(jù)點必定有by＞bx,即(by/bx)＞1.0;反之,對角線下方的點必定有by

圖3 TM3,TM4二維散點圖Fig.3 Planar scatter plot of band TM3 and TM4

(4)對于比值運算,例如,TM3/TM1,則等同于TM3為by,TM1為bx,且有 TM3/TM1=by/bx=tgθ,這里tgθ為坐標平面點P(bx,by)與原點P(0,0)聯(lián)線的斜率。

基于二維散點圖的上述重要性,所以,它將作為對圖像光譜數(shù)據(jù)空間幾何結(jié)構(gòu)特征分析的重要技術(shù)手段。

本區(qū)遙感數(shù)據(jù)的主要波段基于二維散點圖的空間幾何結(jié)構(gòu)分析如下。

圖3的背景主軸是沿著 TM4波段方向的,所以,該工作區(qū)的主要背景地物是植被。其他地物(含蝕變信息)僅構(gòu)成一個次橢圓。

圖4表明,散點信息全在對角線下方,從物理意義而言,表明 TM3與 TM1的比值未能反映出鐵化異常蝕變信息。

結(jié)合圖3與圖4,說明該地區(qū)的3價鐵(Fe3+)的吸收譜帶可能主要在0.90μm(TM4(0.76～0.90)),而不在 0.45μm(TM1(0.45～0.52))。TM3與 TM4能反映出鐵化異常蝕變信息。

圖5為 TM5與 TM7的二維散點圖,由于植被與泥化蝕變信息均為 TM5值高、TM7值低,所以主要包含植被信息的背景中心在對角線上方,說明植被的干擾還是比較嚴重的。

圖6是 TM5與 TM6的二維散點圖,利用 TM6與TM5的比值來提取硅化信息。

圖4 TM3,TM1二維散點圖Fig.4 Planar scatter plot of band TM3,TM1

以上二維散點圖所反映的主要地物在光譜空間的聚類結(jié)構(gòu),表明該地區(qū)提取蝕變異常信息主要受植被影響最大,其次是云區(qū)干擾。蝕變異常信息未能形成自身的“聚類”形態(tài),而是與干擾一起構(gòu)成“點群”,說明蝕變異常信息在空間上分布較分散,大片異常不多,是“小”而“碎”的異常較多。

經(jīng)分析認為,本區(qū)的植被、水體及云區(qū)是3種主要的干擾地物。

2.4 礦化蝕變圖像的識別與提取

通過對遙感圖像光譜數(shù)據(jù)的綜合分析,選擇基于Crosta法的F因子與 H因子分析方法,根據(jù)不同蝕變信息提取的波段選擇,以自行開發(fā)的軟件對工作區(qū)遙感圖像進行礦化蝕變信息提取。

(1)鐵化信息提取。根據(jù) Crosta法先選擇 F(鐵化)因子的 TM1,TM3,TM4,TM5等4個波段組合計算,得到特征向量矩陣(表5)。

經(jīng)分析比較,第3主成分(PCA3)的載荷 TM3和 TM5為正值,TM1和 TM4為負值,即 TM3與TM1,TM4載荷符號相反。所以,PCA3應(yīng)該屬于鐵化蝕變信息主成分。

表5 研究區(qū)TM 1,3,4,5波段的特征向量Table 5 Characteristic vector of band TM1,3,4,5 of the study area

在提取軟件中對濾波處理后的PCA3主成分影像進行密度分割,通過分析最優(yōu)分割段內(nèi)離差平方總和隨分割段數(shù)變化曲線發(fā)現(xiàn),當分割段數(shù)達到14后,曲線趨于平衡,因此14為合理分割段數(shù),分割段數(shù)14的段內(nèi)離差平方總和為28 108 033.71,最優(yōu)分割區(qū)間為:{255,238}{237,219}{218,200}{199,183}{182,162}{161,145}{144,127}{126,109}{108,93}{92,78}{77,62}{61,47}{46,1}{0,0}。

分析對比認為,第1、第2分割段(灰度級為219～255)為鐵化蝕變信息異常區(qū),對該灰度段作進一步的3段最優(yōu)密度分割,并分別賦以紅、黃、綠色,該灰度段以外的灰度級作為背景處理,賦為白色,得到鐵化蝕變信息異常圖(圖7)。

(2)泥化信息提取。選擇 H(泥化)因子的TM1,TM4,TM5,TM7等4個波段組合計算,得到特征向量矩陣(表6)。

由分析比較,第4主成分(PCA4)的載荷 TM5為負值,TM7為正值,即 TM5與 TM7的載荷符號相反。所以,PCA4應(yīng)該屬于泥化蝕變信息主成分。根據(jù)PCA4取反向的數(shù)值,TM5與 TM7提取的泥化蝕變信息圖像見圖8。具體軟件操作與鐵化提取類似,采用最優(yōu)彩色分割對泥化信息進行提取。

表6 研究區(qū)TM 1,4,5,7波段的特征向量Table 6 Characteristic vector of band TM1,4,5,6 of the study area

圖7 TM1,3,4,5的 PCA3主成分提取的鐵化信息圖像Fig.7 The extracted iron-based information Images of principal component PCA3 of TM1,3,4,5

(3)硅化信息提取。根據(jù)前文的分析,這里利用ETM6/ETM5來提取硅化信息。

2.5 礦化蝕變信息提取結(jié)果分析

從圖7、圖8、圖9中可以看出主要蝕變區(qū)集中于圖幅的南部,但在圖8與圖9中圖幅的東北部泥化、硅化均表現(xiàn)較弱,但卻有一致性。

(1)結(jié)合構(gòu)造解譯可見,圖幅南部的鐵化、泥化、硅化均有所反映,它正好與宿務(wù)島火山弧中的一個環(huán)形構(gòu)造帶相吻合,在環(huán)形構(gòu)造的中心即是托萊多斑巖銅礦(大型)的產(chǎn)出部位。托萊多斑巖銅礦具有良好的蝕變分帶:中央為黑云母的鉀化帶,其外側(cè)為硅化較強的石英絹云母化帶,再外則是泥化帶,最外側(cè)是青磐巖化帶。鉀化帶與青磐巖化帶的鐵質(zhì)都比較高。因為鉀化帶的黑云母、青磐巖化中的綠泥石和黃鐵礦均富含鐵質(zhì),所以會有較強的鐵化蝕變顯示。在石英絹云母化帶內(nèi),SiO2是主要的成分,必然有硅化的顯示,而泥化帶的主要成分是高嶺石、伊例石、葉臘石等,故而出現(xiàn)泥化反映。

(2)在圖幅的東北角也有一個較小的環(huán)形構(gòu)造和一組放射狀斷裂,其南側(cè)還有幾個次火山巖巖筒,這里的鐵化較弱,但泥化、硅化卻較明顯,雖然強度 不高,但鐵化、泥化、硅化也有吻合性,可能也是具有銅礦化的一種蝕變標志,可以列為找礦標志。

(3)東南角大片的泥化與硅化可能為濱海沙灘,那里可能有大量的石英和泥質(zhì)的海灘,他們并不是礦化的反映。

3 結(jié)語

將蝕變信息提取結(jié)果與地質(zhì)礦產(chǎn)資料相比較,發(fā)現(xiàn)區(qū)內(nèi)己知的礦床(點)和礦化蝕變帶均在所提取的蝕變信息異常內(nèi);通過野外實地調(diào)查表明,采用的礦化蝕變圖像識別與提取技術(shù)是有效和可靠的。

[1] 趙鵬大,陳永清,劉吉平,等.地質(zhì)異常成礦預(yù)測理論與實踐[M].武漢:中國地質(zhì)大學(xué)出版社,1999.

[2] 張遠飛,朱谷昌,吳德文.地質(zhì)礦產(chǎn)調(diào)查的遙感蝕變信息多層次分離提取技術(shù)與應(yīng)用[C]∥中國遙感應(yīng)用協(xié)會2007年學(xué)術(shù)年會論文集,2007.

[3] 張遠飛,吳德文,朱谷昌,等.遙感蝕變信息檢測中背景與干擾問題研究[J].國土資源遙感,2008(2):22-26.

[4] 張遠飛,楊自安,朱谷昌,等.遙感圖像蝕變信息檢測中的光譜數(shù)據(jù)空間結(jié)構(gòu)分析[J].遙感信息,2009(1):3-9.

[5] 張遠飛,吳健生.基于遙感圖像提取礦化蝕變信息[J].有色金屬礦產(chǎn)與勘查,1999,8(6):604-606.

[6] 張遠飛,楊自安,張普斌,等.高(多)光譜數(shù)據(jù)的背景-異常子空間模型研究[J].地球信息科學(xué)學(xué)報,2009,11(3):283-285.

[7] 張遠飛,吳德文,張艮中,等.高光譜數(shù)據(jù)的波段序結(jié)構(gòu)分析與應(yīng)用研究[J].國土資源遙感,2010(1):30-38.

[8] 張玉君,曾朝銘,陳薇.ETM+(TM)蝕變遙感異常提取方法研究與應(yīng)用——方法選擇和技術(shù)流程[J].國土資源遙感,2003(2):44-49.

[9] 王建平.基于遙感的河南盧氏西部地區(qū)蝕變信息提取與分析[J].地球信息科學(xué)學(xué)報,2007(6):111-115.

[10] 李智勇,郁文賢,匡綱要,等.基于高維幾何特性的高光譜異常檢測算法研究[J].遙感技術(shù)與應(yīng)用,2003(6):379-383.

[11] 李智勇,匡綱要,郁文賢,等.基于高光譜圖像主成分分量的小目標檢測算法研究[J].紅外與毫米波學(xué)報,2004(4):286-290.

[12] 方洪賓,李志中.遙感化探信息綜合分析在地質(zhì)找礦中的應(yīng)用研究[J].國土資源遙感,1998(4):33-36.

[13] 陳彩芬,舒寧.SAR影像與 TM影像的幾種融合處理方法[J].國土資源遙感,1999(4):53-57.

[14] 吳德文,朱谷昌,吳健生,等.青海芒崖地區(qū)巖石光譜特征分析及應(yīng)用[J].國土資源遙感,2001(4):28-34.

[15] 趙元洪,張福祥.波段比值的主成分復(fù)合在熱液蝕變信息提取中的應(yīng)用[J].國土資源遙感,1991(3):12-26.