宋 坤 王恩德 付建飛 姚玉增 孫 健

(1.東北大學(xué)資源與土木工程學(xué)院,遼寧 沈陽 110819;2.中國海洋大學(xué)海底科學(xué)與探測技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室,山東 青島 266100)

遼寧鞍本地區(qū)弓長嶺鐵礦是我國唯一由鞍山式貧鐵礦經(jīng)后期熱液改造形成的大型磁鐵富礦,也是唯一具有工業(yè)開采價值的富鐵礦[1]。國內(nèi)外諸多學(xué)者在弓長嶺鐵礦成因方面進(jìn)行了深入研究,其找礦方法以地質(zhì)、地球物理及地球化學(xué)方法為主,對于遙感方法的應(yīng)用較少。自20世紀(jì)70年代以來,遙感技術(shù)以其快速、無損的優(yōu)勢被廣泛應(yīng)用在地質(zhì)領(lǐng)域[2]。相對于較傳統(tǒng)的地質(zhì)找礦方法在礦產(chǎn)勘探工作中的高成本、低效率等不足來說,遙感技術(shù)可在較大范圍內(nèi)提取研究區(qū)蝕變礦化信息,為礦產(chǎn)資源潛力評估提供幫助,為找礦勘探工作提供方向[3]。本研究以弓長嶺鐵礦為例,選取Landsat8OLI遙感數(shù)據(jù),以羥基礦物和鐵染礦物在特定波段的吸收或反射特征為理論依據(jù),應(yīng)用主成分分析法和波段比值法對遙感異常蝕變信息進(jìn)行提取。通過總結(jié)并分析蝕變信息與鐵礦體之間的關(guān)系,為鞍本地區(qū)及其他同類型礦床的找礦工作提供依據(jù)。

1 研究區(qū)概況

1.1 研究區(qū)地質(zhì)背景

研究區(qū)位于弓長嶺及其周邊區(qū)域,大致以弓長嶺為中心區(qū)域,向西至遼陽市東側(cè),向東至南芬鐵礦區(qū),呈70 km×50 km的矩形,面積約3 500 km2(圖 1)。弓長嶺鐵礦資源十分豐富,具有儲量大、分布廣等特點(diǎn)。研究區(qū)位于華北克拉通北緣東段的遼東地區(qū),具有下部太古宙基底與上部古生界沉積蓋層的雙層結(jié)構(gòu)?；讕r性主要由鞍山群茨溝組、大峪溝組和櫻桃園組變質(zhì)巖組成[4]。我國賦存于太古宙變質(zhì)巖系的鐵礦最早發(fā)現(xiàn)于這一帶,稱之為鞍山式鐵礦,以貧礦多、少有富礦、儲量大為主要特點(diǎn)。據(jù)統(tǒng)計,遼寧省鐵礦資源儲量占全國總量的1/4左右,且鞍山式鐵礦占遼寧省鐵礦總量的96%。

圖1 研究區(qū)區(qū)域地質(zhì)特征Fig.1 Regional geological characteristics of the study area

1.2 礦床地質(zhì)特征

鞍本地區(qū)以弓長嶺礦區(qū)規(guī)模最大,也最為典型。弓長嶺礦區(qū)主要為沉積變質(zhì)型鐵礦,磁鐵富礦體一部分分布在各層條帶狀鐵礦床中,還有部分含鐵礦體產(chǎn)在蝕變圍巖或混合巖中[5]。這表明富鐵礦體可能是條帶狀磁鐵礦體經(jīng)熱液改造形成[6]。礦集區(qū)呈NW—SE向展布,弓長嶺一、二、三礦區(qū)以及獨(dú)木礦區(qū)為鐵礦床的主要分布區(qū)域(圖2)。鞍山群茨溝組是該地區(qū)的主要賦礦地層,二礦區(qū)茨溝組地層出露最為完整。弓長嶺背斜為控礦構(gòu)造[7]。礦區(qū)內(nèi)礦石類型主要為磁鐵石英巖,其中二礦區(qū)內(nèi)含有富鐵礦石(圖2)。出露的巖漿巖主要為晚太古代混合花崗巖,弓長嶺礦區(qū)呈孤島狀分布在混合花崗巖中[8]。

圖2 弓長嶺礦區(qū)地質(zhì)特征Fig.2 Geological characteristics of Gongchangling mining area

2 遙感數(shù)據(jù)處理

2.1 數(shù)據(jù)源

Landsat 8衛(wèi)星于2013年成功發(fā)射,衛(wèi)星上搭載了OLI陸地成像儀和TIRS熱紅外傳感器,其空間分辨率和光譜分辨率等方面與之前的Landsat系列衛(wèi)星基本保持一致。Landsat 8與Landsat 7衛(wèi)星波段對應(yīng)關(guān)系參見表1。Landsat 8遙感影像數(shù)據(jù)分辨率較高,且影像涉及范圍廣,在遙感異常信息提取中被廣泛使用;同時,OLI和ETM+數(shù)據(jù)波段分配大體一致,且OLI波段范圍比ETM+更窄,可更好地表現(xiàn)出蝕變礦物波的譜特征[9]。

表1 Landsat 8與Landsat 7衛(wèi)星波段參數(shù)Table 1 Satellite bands parameters of Landsat 8 and Landsat 7

通過對TM和ETM+數(shù)據(jù)波段特征及其應(yīng)用的分析研究,總結(jié)出OLI傳感器中不同波段的光譜特征及其適用領(lǐng)域,即Band 5、Band 6和Band 7 3個波段與蝕變異常提取相關(guān),其中,對碳酸鹽以及黏土礦物等有明顯吸收特征的光譜波段為Band 7。

研究區(qū)Landsat 8 OLI圖像的成像時間為2015-04-20,衛(wèi)星軌道號為119-31,經(jīng)反復(fù)比對分析,該時段圖像地物輪廓清晰、層次感強(qiáng),圖上無云層,地面無積雪,植被覆蓋率低,地表巖石、土壤裸露較好,成像質(zhì)量最佳,有效滿足了本次研究的要求(圖3)。研究區(qū)圖像范圍為東經(jīng) 123°04′40″～123°58′01″,北緯40°53′14″～41°25′54″。

圖3 研究區(qū)原始遙感影像Fig.3 Original remote sensing image of the study area

2.2 遙感影像預(yù)處理

在遙感成像過程中,遙感影像因受其搭載平臺、傳感器誤差以及大氣環(huán)境、太陽高度角、地形等不同因素影響,得到的測量值與地物本身的光譜反射率存在誤差。為了消除或修正輻射畸變,得到真實(shí)的地物反射率,須對原始影像數(shù)據(jù)進(jìn)行輻射校正,通常包括傳感器校正和大氣校正兩類[10]。

對傳感器誤差的校正是通過輻射定標(biāo)進(jìn)行,將傳感器中的電壓或數(shù)字量化值(Digital Number,DN)換算為大氣外層的表面反射率。大氣校正是通過將大氣頂層的輻射亮度值轉(zhuǎn)換為地表反射率消除或修正誤差,從而獲取地物的真實(shí)反射率(圖4)。上述處理均使用ENVI4.7軟件的輻射定標(biāo)和FLAASH大氣校正模塊,參照遙感數(shù)據(jù)的參數(shù)設(shè)置完成校正。

圖4 大氣校正前后遙感圖像及其波譜曲線Fig.4 Remote sensing images and their pop curves before and after atmospheric correction

3 研究方法

不同的礦物、巖石及礦化蝕變信息都具有其特有的光譜特征[11]。依據(jù)蝕變礦物的不同波譜特征,可以有效提取遙感蝕變異常信息[12]。通常情況下,蝕變異常信息提取的主要相關(guān)礦物離子包括鐵離子、羥基或碳酸根等離子(表2)。

表2 對巖石反射光譜特征起主導(dǎo)作用的離子和基團(tuán)的光譜特征[11]Table 2 Spectral characteristics of ions and groups that dominate the spectral characteristics of rock reflections

分析表2可知:

(1)含鐵(Fe2+、Fe3+)離子。以含F(xiàn)e3+的褐鐵礦、赤鐵礦等蝕變礦物為主,其中含F(xiàn)e2+的蝕變礦物相對較少。如圖5所示,Fe3+礦物在0.40～0.55μm和0.85～0.95μm波段表現(xiàn)為明顯的吸收特征,對應(yīng)OLI 2(0.45～0.51μm)波段、OLI 3(0.53～0.59μm)波段和OLI 5(0.85～0.88μm)波段的強(qiáng)吸收特征,以及在OLI 4(0.64～0.67μm)波段的強(qiáng)反射特征。

圖5 含鐵類蝕變礦物反射光譜曲線Fig.5 Reflectance spectral curves of iron dyeing altered minerals

(2)含羥基(OH-)或碳酸根(CO32-)離子。該類離子蝕變礦物主要為黏土礦物,如圖6所示,在1.4 μm及2.2～2.4μm波段存在吸收特征,對應(yīng)OLI 7(2.11～2.29μm)波段產(chǎn)生低反射值,在 OLI 6(1.57～1.65μm)波段產(chǎn)生高反射值。含碳酸根礦物主要存在5個特征吸收譜帶,分別在2.35μm、2.55μm波段處表現(xiàn)出較強(qiáng)吸收,在1.9、2.0、2.16 μm波段處表現(xiàn)出相對較弱吸收,對應(yīng)OLI 7(2.11～2.29μm)波段形成吸收特征。

圖6 含羥基類蝕變礦物反射光譜曲線Fig.6 Reflectance spectral curves of hydroxyl altered minerals

3.1 波段運(yùn)算法

波段運(yùn)算法包括波段的加減運(yùn)算與比值運(yùn)算。波段加法通過將多幅圖像亮度值相加求平均,達(dá)到圖像視覺增強(qiáng)效果、邊緣增強(qiáng)效果。圖像減法運(yùn)算則可以去背景,突出研究對象。蝕變圍巖與周圍巖石之間的亮度值存在差異,可通過波段加減法運(yùn)算增強(qiáng)礦化蝕變信息[13]。

為了增強(qiáng)蝕變信息,還可以運(yùn)用波段比值法。根據(jù)每個象元在兩個不同波段上的亮度比值形成新的比值圖像[14],可消除地物反射或避免光照引起的陰影干擾,突出重要信息。

比值法通常包括簡單比值法、差和比值法以及截取比值法等。結(jié)合弓長嶺礦區(qū)實(shí)際情況,突出異常信息,本研究運(yùn)用簡單比值法對圖像進(jìn)行處理[15]。簡單比值法運(yùn)算公式為

式中,DNm(x,y),DNn(x,y)分別為像元 (x,y)在m和n波段上的亮度值;Rmn(x,y)為輸出的亮度比值。

OLI遙感影像通?？蓪﹁F氧化物、氫氧化物類,羥基礦物類,碳酸鹽、泥化類3類蝕變礦物進(jìn)行識別[16]。在遙感異常信息提取過程中增強(qiáng)蝕變信息的主要OLI波段比值為:①Band 4/Band 2,主要增強(qiáng)鐵氧化物、氫氧化物類蝕變;②Band 6/Band 4,主要增強(qiáng)硅化類蝕變;③Band 6/Band 7,主要增強(qiáng)碳酸鹽化、泥化類蝕變。

本研究利用ENVI 4.7軟件中的“band math”工具進(jìn)行上述相關(guān)波段進(jìn)行運(yùn)算,結(jié)果如圖7所示。

圖7 波段比值法圖像處理結(jié)果Fig.7 Image processing results by band ratio method

3.2 主成分分析法

主成分分析法(Principal Components Analysis,PCA)是較常用的蝕變信息提取方法。通過對反應(yīng)蝕變信息的多個波段進(jìn)行特征統(tǒng)計的多維正交線性變換(K-L變換),生成1組新的組分圖像[17]。主要目的是集中多波段中有用的信息到盡可能少的新組分圖像中,確保地物信息之間互不干擾再去相關(guān)性[18]。同時,選用計算結(jié)果中的新波段作為多波段原始圖像,可實(shí)現(xiàn)增強(qiáng)重要信息的目的[19]。

主成分分析法在遙感蝕變異常信息提取中應(yīng)用廣泛,被稱為CROSTA法。根據(jù)前人研究成果,對于OL I數(shù)據(jù),常用OLI 2、OLI 5、OLI 6、OLI 7這4個波段提取羥基蝕變異常,用OLI 2、OLI 4、OLI 5、OLI 6這4個波段提取鐵染遙感蝕變異常。

鐵染異常主分量的判斷準(zhǔn)則為含鐵染異常礦物OLI 4與OLI 2及OLI 5的系數(shù)符號相反、OLI 4與OLI 6的系數(shù)符號相同[20]。同時,還可以根據(jù)不同波段在不同主分量上的載荷情況來輔助確定異常主成分(該主分量在OLI 4波段載荷較大,在OLI 2波段載荷較小)。本研究將與上述準(zhǔn)則對應(yīng)的主分量稱為鐵染異常主分量。

羥基異常主分量的判斷準(zhǔn)則為含羥基異常礦物的OLI 6系數(shù)應(yīng)與OLI 7及OLI 5的系數(shù)符號相反,OLI 2一般與OLI 6系數(shù)符號相同[21-22],且該主分量在OLI 6、OLI 7波段載荷較大。本研究將與上述準(zhǔn)則對應(yīng)的主分量稱為羥基異常主分量。

3.2.1 鐵染蝕變異常信息提取

根據(jù)上述PCA方法提取原理,對研究區(qū)圖像的OLI 2、OLI4、OLI 5、OLI6波段進(jìn)行主成分分析,得到4個主分量,并對其進(jìn)行統(tǒng)計分析,結(jié)果見表3。

表3 OLI 2、4、5、6波段提取鐵染異常特征向量矩陣及各主分量Table 3 Eigenvector matrix and eigenvalues of each principal component of iron dyeing anomaly with OLI 2,4,5 and 6

表3主要體現(xiàn)了4個波段的亮度信息,其特征向量值在4個波段中均為負(fù)值,不符合判斷準(zhǔn)則。同時,根據(jù)上述鐵染異常主分量判斷準(zhǔn)則,發(fā)現(xiàn)并沒有完全符合條件的主成分。但是第4主成分(PC4)最接近判斷準(zhǔn)則,其在 OLI 4上的載荷系數(shù)(0.609 854)與在OLI 2上的載荷系數(shù)(-0.782 576)符號相反,而在OLI 5上的載荷系數(shù)(0.002 123)雖然是正,但數(shù)值很小接近于0?？紤]到獲取圖像和預(yù)處理時可能會帶入一些誤差,因此可以將PC4視為滿足判斷準(zhǔn)則的主分量。從含鐵礦物光譜特征來看,含鐵蝕變巖石在OLI 2強(qiáng)吸收、OLI 4強(qiáng)反射,因此含鐵染信息的主分量在OLI 2和OLI4顯示為強(qiáng)吸收和強(qiáng)反射且在這兩個波段上具有較大的載荷,可以看出PC4主分量滿足判斷準(zhǔn)則,代表了鐵染蝕變信息。

通過分析PC4灰度圖像發(fā)現(xiàn),與實(shí)際礦點(diǎn)對比發(fā)現(xiàn)異常區(qū)域?yàn)榈椭?顯示為暗),因此對PC4數(shù)據(jù)取相反數(shù),使得鐵染異常區(qū)顯示為高值(亮),取反后的-PC4灰度圖像如圖8所示。

圖8 鐵染異常主分量DN值取反圖像Fig.8 Principal component DN value inverse image of iron dyeing anomaly

3.2.2 羥基蝕變異常信息提取

根據(jù)上述PCA方法提取原理,對研究區(qū)圖像的OLI 2、OLI5、OLI6、OLI7波段進(jìn)行主成分分析,得到4個主分量,并對其進(jìn)行了統(tǒng)計分析,結(jié)果見表4。

表4 OLI 2、5、6、7 提取羥基異常特征向量矩陣及各主分量Table 4 Eigenvector matrix and eigenvalues of each principal component of hydroxyl anomaly with OLI 2,5,6 and 7

由表4可知:根據(jù)上述羥基異常主分量判斷準(zhǔn)則,該主分量在OLI 6、OLI 7波段載荷較大。同時,可以發(fā)現(xiàn) PC4主分量符合判定準(zhǔn)則,其 OLI 6(0.617 392)系數(shù)與OLI 7(-0.543 045)及OLI 5(-0.331 771)系數(shù)符號相反,OLI2(0.462447)與OLI 6(0.617 392)系數(shù)符號相同。PC4特征向量載荷系數(shù)絕對值較大的是OLI 6和OLI 7,分別為0.617 392和0.543 045,二者符號相反,而其他波段載荷較小,表明PC4主成分的信息主要來自O(shè)LI 6和OLI 7,與黏土礦物和碳酸鹽礦物存在有關(guān),且與OLI 6高反射和OLI7強(qiáng)吸收相符,突出了泥化蝕變信息。因此,本研究選擇PC4作為羥基信息。結(jié)合研究區(qū)實(shí)際情況發(fā)現(xiàn)灰度圖像中亮的區(qū)域?yàn)榱u基異常區(qū)(圖 9)。

圖9 羥基異常主分量圖像Fig.9 Principal component image of hydroxyl anomaly

4 弓長嶺礦區(qū)蝕變異常信息提取

4.1 鐵染異常

首先對OLI 2、4、5、6 主成分分析得到的鐵染異常主分量-PC4數(shù)據(jù)進(jìn)行數(shù)學(xué)統(tǒng)計(表5),利用門限法分級,多次試驗(yàn)確定最優(yōu)分割(表6),得到如圖10所示的鐵異常分級結(jié)果。

表5 鐵染異常主分量-PC4數(shù)據(jù)數(shù)學(xué)統(tǒng)計信息Table 5 Mathematical statistics information of the principal component-PC4 data of iron dyeing anomaly

表6 鐵染異常主分量異常分級Table 6 Classification of principal component of iron dyeing anomaly

圖10 弓長嶺地區(qū)鐵異常圖Fig.10 Iron anomaly map of Gongchangling Area

4.2 羥基異常

對OLI 2、5、6、7主成分分析得到的羥基異常主分量PC4數(shù)據(jù)進(jìn)行數(shù)學(xué)統(tǒng)計(表7),利用門限法分級,多次試驗(yàn)確定最優(yōu)分割(表8),得到如圖11所示的羥基異常分級結(jié)果。

表7 羥基異常主分量PC4數(shù)據(jù)數(shù)學(xué)統(tǒng)計信息Table 7 Mathematical statistics information of the principal component PC4 data of hydroxyl anomaly

表8 羥基異常主分量異常分級Table 8 Classification of principal component of hydroxyl anomaly

圖11 弓長嶺地區(qū)羥基異常圖Fig.11 Hydroxyl anomaly map of Gongchangling Area

針對弓長嶺礦區(qū)范圍進(jìn)行詳細(xì)的蝕變異常分級,觀察其在小范圍內(nèi)(礦區(qū))是否能表達(dá)出更多、更有規(guī)律的信息。將主成分分析得到的鐵染和羥基異常主分量進(jìn)行裁剪,只保留弓長嶺礦區(qū)圖像。對裁剪得到的圖像數(shù)據(jù)重新進(jìn)行數(shù)學(xué)統(tǒng)計,然后利用門限法重新分級[23]。蝕變信息特征中大部分礦區(qū)以外的背景值、假異常等干擾因素可通過處理去除,結(jié)果如圖12和圖13所示。

圖12 弓長嶺礦區(qū)鐵異常圖Fig.12 Iron anomaly map of Gongchangling Mning Area

圖13 弓長嶺礦區(qū)羥基異常圖Fig.13 Hydroxyl anomaly map of Gongchangling Mining Area

4.3 分析驗(yàn)證

主成分分析法對蝕變信息的提取效果比較理想,研究區(qū)內(nèi)主要礦區(qū)都有異常顯示,能明顯的與非異常區(qū)分。鐵和羥基的蝕變異常區(qū)域重疊較好,也吻合了研究區(qū)內(nèi)礦化區(qū)域的鐵和羥基含量均較高的特征。為了進(jìn)一步驗(yàn)證遙感蝕變異常信息與區(qū)域礦化信息的吻合程度,本研究將鐵和羥基蝕變異常進(jìn)行疊加處理,并將研究區(qū)內(nèi)已知礦點(diǎn)投影到異常疊加圖上,分析蝕變異常與礦點(diǎn)的吻合程度,結(jié)果如圖14所示。

圖14 蝕變異常疊加圖和部分已知礦點(diǎn)投影Fig.14 Alteration anomaly superposition map and partial known ore occurrences projection map

由圖14可知:投影的7個超大型礦床和6個大型礦床中有5個超大型礦床和4個大型礦床位于強(qiáng)蝕變異常區(qū)范圍內(nèi),吻合程度較好,2個大型礦床處于輕微蝕變異常區(qū)范圍內(nèi),可能是礦床圍巖蝕變作用較弱或?yàn)殡[伏礦床所致。2個超大型礦床不在明顯的蝕變異常區(qū)內(nèi),吻合程度較差,分別為大臺溝礦床和思山嶺礦床。大臺溝礦床2011年詳查顯示埋深達(dá)到1 100～1 400m,含鐵巖系分布在太古宙鞍山群,上覆遼河群及近水平的青白口系、震旦系等沉積蓋層[24]。思山嶺鐵礦埋深達(dá)400 m以上,賦礦巖系為鞍山群茨溝組,上覆沉積蓋層為青白口系和震旦系巖層[25]。兩個礦床的變質(zhì)蝕變信息均被沉積蓋層覆蓋,影響蝕變信息提取結(jié)果。除了以上兩個超大埋深的隱伏礦床外,其余礦床與蝕變異常區(qū)域基本吻合,證明主成分分析法能有效地提取礦化蝕變異常信息,在同類型礦床中具有較好的應(yīng)用前景。

5 結(jié) 論

(1)以Landsat 8遙感數(shù)據(jù)為基礎(chǔ),通過提取弓長嶺礦區(qū)蝕變礦物信息,在同類型BIF鐵礦床中建立了一種快速準(zhǔn)確的遙感找礦勘查方法。

(2)采用主成分分析法和波段比值法分別對研究區(qū)以及弓長嶺礦區(qū)進(jìn)行蝕變異常信息提取,在研究區(qū)內(nèi)主成分分析法可以有效地抑制干擾信息,提取鐵和羥基蝕變信息效果較好。弓長嶺礦區(qū)鐵異常主要分布在礦區(qū)中心且向南延展。羥基異常范圍較鐵染異常面積大,且礦區(qū)南部有較多弱異常分布,其中有部分信息為地表覆蓋下的假異常。

(3)提取的蝕變異常信息與研究區(qū)內(nèi)已知礦點(diǎn)對比發(fā)現(xiàn),異常區(qū)與礦點(diǎn)吻合較好。蝕變信息集中區(qū)域?yàn)閷φ业V有指示意義的區(qū)域,反映出本研究方法對于鞍本地區(qū)及其他同類型礦床研究有一定的參考價值。但該方法在提取地表覆蓋層較厚的隱伏礦床蝕變信息過程中,對假異常的識別和蝕變?nèi)跣畔⒌脑鰪?qiáng)仍存在不足。下一步應(yīng)結(jié)合鞍本地區(qū)的其他礦床信息進(jìn)行補(bǔ)充和完善,進(jìn)一步提升該方法的普適性。