菲律賓宿務島地區礦化蝕變圖像識別與提取分析

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

(1.中南大學地學與環境工程學院,長沙410083;2.有色金屬礦產地質調查中心,北京100012)

0 引言

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

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

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

1 宿務島礦區地質背景與數據庫

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

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

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

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

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

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

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

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

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

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

M是主成分的本征數,背景(含大干擾)的端元數目應該等于M+1。分析表2中Q(M)值的變化會發現,這些值具有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值的估計函數 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個波段組合計算,得到標準差與方差分數(即貢獻率)(表3)、M值的估計函數 Q(M)值(表4)。

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

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

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

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

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

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

回歸偏度曲線是用來描述線性回歸的二維散點圖的圖形不對稱性,回歸偏度較大,則預示著有不同于背景的異?；蚋蓴_信息存在[7]。

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

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

2.3 礦化蝕變的波譜空間幾何結構分析

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

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

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

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

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

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

圖2 光譜數據的二維散點圖Fig.2 Planar scatter plot of spectral data

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

(3)若約定圖形坐標系統的Bx(橫坐標)表示具有吸收峰的波段,比如 TM1或 TM7,而By(縱坐標)表示具有反射峰的波段,如 TM3或 TM5。根據幾何意義,坐標平面內的點P(bx,by),位于對角線(圖2中綠色線)上方的數據點必定有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)聯線的斜率。

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

本區遙感數據的主要波段基于二維散點圖的空間幾何結構分析如下。

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

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

結合圖3與圖4,說明該地區的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

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

經分析認為,本區的植被、水體及云區是3種主要的干擾地物。

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

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

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

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

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

在提取軟件中對濾波處理后的PCA3主成分影像進行密度分割,通過分析最優分割段內離差平方總和隨分割段數變化曲線發現,當分割段數達到14后,曲線趨于平衡,因此14為合理分割段數,分割段數14的段內離差平方總和為28 108 033.71,最優分割區間為:{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)為鐵化蝕變信息異常區,對該灰度段作進一步的3段最優密度分割,并分別賦以紅、黃、綠色,該灰度段以外的灰度級作為背景處理,賦為白色,得到鐵化蝕變信息異常圖(圖7)。

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

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

表6 研究區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)硅化信息提取。根據前文的分析,這里利用ETM6/ETM5來提取硅化信息。

2.5 礦化蝕變信息提取結果分析

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

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

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

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

3 結語

將蝕變信息提取結果與地質礦產資料相比較,發現區內己知的礦床(點)和礦化蝕變帶均在所提取的蝕變信息異常內;通過野外實地調查表明,采用的礦化蝕變圖像識別與提取技術是有效和可靠的。

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

[2] 張遠飛,朱谷昌,吳德文.地質礦產調查的遙感蝕變信息多層次分離提取技術與應用[C]∥中國遙感應用協會2007年學術年會論文集,2007.

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

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

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

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

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

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

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

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

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

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

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

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

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