吳德文,張遠飛,袁繼明,楊自安,張建國
(有色金屬礦產地質調查中心,北京 100012)
高光譜數據以其圖譜特征在巖石礦物識別方面具有獨特的優勢,特別是航空高光譜數據,兼具高空間分辨率和高光譜分辨率,可以實現對礦物的精細識別和定量提取。隨著高光譜遙感技術的發展,航空高光譜礦物識別與填圖技術在地質礦產勘查中逐漸發揮了重要作用,已成為當前國內外遙感地質應用領域研究的一個熱點,并取得了顯著的應用成效(甘甫平等,2003;葉發旺等,2018)。但這些研究和應用主要集中在典型礦物的光譜特征、光譜識別與蝕變礦物填圖方面,而對巖相-巖性的識別,乃至全巖信息提取和地質填圖則較少涉及。實際上,巖相-巖性、構造的光譜識別是智能化地質填圖的基礎,對區域地質背景和成礦地質條件研究,以及快速進行找礦預測與礦產資源潛力評價同樣具有重要意義。本團隊長期開展高(多)光譜遙感技術在金屬礦產勘查中的應用研究,注重于技術方法總結以及新技術方法和數據模型研發,在多源地學信息挖掘方面,逐步形成了一些新的技術思路和技術方法,并在自主開發的軟件系統中得以實現。尤其是構建起高光譜巖相、礦物填圖的技術框架和方法體系,在地質礦產調查中取得了良好的應用效果,本文以新疆東疆地區白干湖幅(K46E012020)1∶5萬航空高光譜巖相、礦物填圖為例加以介紹。
有色金屬礦產地質調查中心于2014—2015年組織實施中國地質調查局“新疆重點地區航空高光譜調查與找礦預測技術研究”項目時,利用HyMap機載成像儀,在東戈壁—野馬泉西地區飛行獲取了約4000 km2的高光譜數據,本次選取了其中的白干湖幅1∶5萬圖幅范圍的數據開展研究工作。
該次飛行獲取的高光譜數據的空間分辨率為4.5 m,光譜分辨率在可見光近紅外(VNIR)波長范圍約為15 nm,短波紅外(SWIR)波長范圍約為18 nm。在獲取的128個波段數據中,有效波段為125個。產品數據經過了輻射定標、大氣校正(光譜重建)、幾何校正和地理編碼等處理,以及各航帶圖像的匹配和無縫拼接。主要問題是圖像局部存在拼接錯位和色調不匹配(鑲嵌線兩側色差明顯),但采用合適的處理方法,可以在一定程度上消除這些問題帶來的影響,基本滿足本次填圖應用的要求。
遙感圖像的數據質量和信息特征(包括光譜信息量、光譜值的分布形態、背景區域、異常分布等)可通過直方圖及統計特征值(均值、眾數、中值、方差、偏度、峰度等)反映出來。遙感圖像各個波段的灰度直方圖和統計特征值,在多波段圖像分析中可以逐一進行觀察和對比分析,但對于多達一百個以上波段的高光譜圖像,難以做到這一點。為了總體把握與分析高光譜圖像數據特征,可采用波段序列直方圖的描述與表達方式:橫坐標(X軸)是波段序列(即波段序號由小到大排列),縱坐標(Y軸)是輻射值(相當于灰度值),第三維(Z軸)為輻射值出現的頻數,采用彩色分級的方式予以可視化表達。這樣的波段序列直方圖能夠很好地總覽與分析高光譜圖像所有波段的直方圖及統計特征信息。
白干湖幅HyMap數據波段序列直方圖表明(圖1):①HyMap數據各波段的直方圖均呈近于高斯分布(正態分布)的正常形態;②第62波段(1391 nm)和63波段(1406 nm)出現輻射值范圍的突變,是因為這2個波段位于大氣水帶吸收波長上,反映的基本是大氣水汽的特征,極少包含地面信息,故本次被棄用;③序列直方圖頻帶(不同顏色表示)反映出了光譜曲線的整體變化趨勢與形態特征,如波段序列兩端的波段(波段1~9、124~125)的輻射值相對較低,輻射值范圍相對較窄,統計得到的標準差也相對較小,說明這些波段的信息量相對較少,其它波段的信號和噪聲分布正常;④各波段直方圖均存在一定偏度變化(綠色、青色至藍色的各個色級分布區域),分析認為大的偏度往往是干擾因素引起的,而小的偏度可以看作異常信息作用的結果。

圖1 白干湖幅HyMap數據波段序列直方圖
進行遙感蝕變信息提取時,為了減少盲目性,首先需要了解工作區內各類蝕變的存在及其在圖像上的光譜響應,通常情況下,可以參考已有地質資料,把已知地段作為試驗區進行對比分析來加以驗證。但是,在缺少礦產資料或者已知地段不能準確定位的情況下,難以進行蝕變信息的判別和分析。針對這個問題,本團隊提出了波段序結構分析的解決方案,主要技術方法為回歸偏度序結構分析和相關序結構分析,通過遙感數據本身的特征統計信息來評價工作區“蝕變異常”的潛在性和蝕變強度(張遠飛等,2009,2010)。
2.2.1 回歸偏度序結構分析
上面已經提到,單波段直方圖的偏度反映了不同于背景的異常或干擾信息的存在;同樣地,二維散點圖上背景橢圓范圍之外的點群或回歸軸兩側的圖形不對稱現象也反映了異常信息的存在,在統計意義上可用兩個波段的回歸偏度來表征。對高光譜數據進行兩兩波段的回歸偏度統計,用圖形表示就形成回歸偏度曲線圖,可反映波段之間的相似程度漸變的過程,其起伏與大小決定了高光譜波段序列結構的變化規律,并通過波形突變揭示蝕變特征譜帶。
統計獲得白干湖幅HyMap數據全波段回歸偏度曲線(圖2~3),依據曲線的相似度和波形變化規律大體可劃分為如下波段(譜帶)組:1—9波段(471—590 nm)、10—45波段(605—1109 nm)、46—61波段(1124—1335 nm)、64波段(1421 nm)、65—93 波 段(1435—1800 nm)、94—95 波 段(1941—1960 nm)、96—105 波 段(1979—2147 nm)、106—113 波段(2164—2287 nm)、114—118波段(2304—2371 nm)、119—123 波 段(2387—2452 nm)和124—125波段(2467—2483 nm)。圖2和圖3分別為10—45波段組和96—105波段組的回歸偏度曲線圖,這兩組曲線圖較好地反映了區內存在的蝕變礦物的光譜響應。

圖3 白干湖幅HyMap數據96-105波段回歸偏度曲線圖
圖2、3表明,檢測到的特征譜帶主要有:1.42 μm附近的水帶或Al-OH弱吸收譜帶、1.94 μm附近的水帶、2.22 μm附近的Al-OH吸收譜帶、2.25 μm附近的Fe-OH吸收譜帶、2.34 μm附近的Mg-OH吸收譜帶、2.35 μm附近的碳酸根吸收譜帶和2.44 μm附近的伴隨Mg-OH吸收譜帶等。一般來說,回歸偏度值越大,在曲線上的波形變化越顯著,反映的蝕變礦物信息越強。根據上述曲線分析,結合已有地質資料判斷,該圖幅內普遍存在的蝕變類型有綠泥石化和絹云母化,其次有綠簾石化、高嶺石化、方解石化、黃鉀鐵礬化等。
2.2.2 相關序結構分析
相關(似)系數能夠反映兩個波段之間的相關或相似程度。所有相鄰波段之間的相關性是比較高的,高光譜序列波段中兩兩波段之間的相關性通常情況下是按序列方式漸變的,當出現突變時,就表明有“異常”光譜地物存在,這種“異常”光譜地物可能是“干擾”地物或蝕變礦物。與回歸偏度統計類似,對高光譜數據進行兩兩波段的相關系數統計,用圖形表示就得到序列相似分析曲線,其波形突變部位反映了蝕變特征譜帶。
統計獲得白干湖幅HyMap數據全波段相似分析曲線,根據曲線的相似度和波形變化規律分析,獲得了與回歸偏度序結構分析基本相同的檢測結果。
2.3.1 高光譜地質信息的分類提取
常規的高光譜蝕變信息識別與提取方法主要是基于光譜相似度的整體光譜匹配法和基于診斷性光譜吸收譜帶的局部光譜識別法,這些方法對提取某些特定蝕變礦物信息是有效的,但對不同巖性的巖石信息提取卻難以實現,因為巖石一般為不同礦物的集合體,并受諸多地質因素的影響,其光譜復雜而多變,基本上無標準光譜與之匹配。此外,它們對高光譜數據的校正有較高的要求,即重建的光譜應同標準光譜相近或相似,當無法重建光譜或重建的光譜出現嚴重的“變形”時,都會導致無法應用光譜匹配方法。至于多維數據的分類方法,監督分類很難根據已知樣本數據確定出最佳的“分類模式”,而非監督分類的最佳分類數和分類效果也難以估計和達到。本團隊研制的最優動態聚類法可以有效地解決這些問題,實現全巖信息提取,該方法是一種改進與優化的動態集群分類(屬非監督分類),分類算法能自動估算出最佳分類數目和初始類別中心,實現聚類核最佳逼近,讓分類結果達到全局最佳,可稱之為高維聚類核最佳逼近動態集群分類,簡稱最優動態聚類(張遠飛等,2011)。分類過程中全部波段參與,保證了圖像光譜信息的完整性,達到有效識別不同的地物類型。
利用白干湖幅HyMap數據的全部有效波段參與最優動態聚類分類。運算過程中,首先通過各波段直方圖最優分割獲取初始類中心,根據初始聚類的類-頻數曲線(圖4)確定最佳分類數,然后進行動態循環聚類至類中心收斂(最佳逼近)。
從圖4可以看到,按頻數由大到小排序后,前5類已聚類了絕大部分像元,對應于較大范圍分布的地質體;排后的基本為小類,對應于小范圍分布的地質體或局部異常,20類以后只存在極少像元,因此取20為全區大類分類數,將其余小類進行歸并。分類完成后提取各類單元的圖像光譜曲線,根據各類單元中以云母類為主的Al-OH礦物(統稱為絹云母化)和以綠泥石為主的Mg-OH礦物(統稱為綠泥石化)的光譜響應強弱變化進行賦色(表1),得到分類結果圖像(圖5)。

圖4 白干湖幅HyMap數據最優動態聚類的初始類-頻數曲線
2.3.2 分類信息的光譜識別
對照研究區原有地質圖可以看到,高光譜分類圖能夠較好地區分不同的地質體(包括巖性層、巖體、巖脈、松散堆積物等),并明細表現出構造(帶)、蝕變帶等地質現象。各類地質體的巖性特征可依據其光譜特征進行判斷,在軟件中通過實時聯動的圖像與光譜分析來實現,稱之為“圖譜一體化識別”。從高光譜數據中提取各分類單元的圖像光譜曲線(類均值曲線),結合典型巖石的地面光譜測試數據,通過光譜特征參數提取與譜帶特征分析進行巖性識別。圖6為白干湖幅各分類單元的圖像光譜曲線,各類單元的圖像光譜特征及地物類型判別見表2。

表2 白干湖幅分類單元圖像光譜特征簡表

續表

圖6 白干湖幅各分類單元光譜曲線圖
在上述高光譜數據處理分析與信息獲取工作的基礎上,借助高光譜信息地質解譯編制巖相、礦物分布地質圖。航空高光譜信息地質解譯以HyMap數據的全波段分類圖像為基礎圖像,輔以假彩色合成圖像(第120、60、1波段RGB合成)。解譯基本原則是:在地層單元內圈出不同的巖相(建造)單位,按巖性單元圈定侵入巖。地層組群的劃分及巖體時代歸屬以已有區域地質資料為依據,解譯單元根據不同巖石(組合)或建造類型在圖像上的可分性來確定。各解譯單元的解譯標志依據其影像特征建立,以分類圖的顏色或顏色組合(代表地物類別)為主要標志,輔以合成圖像的色調、色彩、影紋和空間分布等特征標志;構造解譯標志主要是線性影像特征。通過野外地質調查和驗證工作,對高光譜信息的地質解譯標志進行補充、修正和完善。
根據各分類單元的蝕變類型和強弱變化(或礦物種類和相對含量)圈定蝕變礦物分布區(帶),以不同的顏色加以區分。白干湖地區的蝕變礦物主要為以絹云母為主的Al-OH礦物、以綠泥石為主的Mg-OH礦物及碳酸鹽礦物等。據此劃分的蝕變區(帶)包括強絹云母化帶、絹云母化帶、弱絹云母化帶、弱絹云母化+綠泥石化帶、強綠泥石化+絹云母化帶、強綠泥石化帶、綠泥石化帶、弱綠泥石化帶、綠泥石化+絹云母化帶、強硅化帶、無蝕變或弱蝕變區、碳酸鹽巖區等。
成果圖的圖件框架結構和圖面表達形式按照相關制圖要求編制。圖7為白干湖幅原區域地質圖與高光譜地質圖的局部對比,從對比圖可以看到,航空高光譜填圖的成果圖信息量更加豐富,對地質體、構造線的勾繪更加精細、完整、全面、客觀。

圖7 白干湖幅原區調地質圖(左)與高光譜地質圖(右)對比(局部)
高光譜遙感技術是一種新的高效、綠色勘查方法技術,通過地質找礦綜合應用,可以快速查明找礦目標區內礦物信息空間分布與組合特征,快速進行巖相-巖性-構造識別,結合成礦理論快速實現找礦定位預測與找礦預測區查證,并能快速對找礦靶區的資源環境進行評價,為后續礦業開發及基地規劃提供相關地質環境信息。在整個勘查過程的綜合調查與評價中,可以大大節省勘查投入。
高光譜遙感地質找礦綜合應用已成為國外礦產勘查重要技術支撐之一。在國內還存在不足,一是受高光譜數據獲取限制,尚以試驗性應用為主,工程化應用還未實質性地開展;二是以數據處理和信息提取技術方法研究為主,與地質找礦和礦產資源勘查評價結合程度低。這些不足加之研究成果的碎片化,使得高光譜遙感技術在礦產勘查中的優勢作用還沒有充分顯現出來。近年來,隨著高光譜數據獲取技術的不斷發展,特別是我國高分五號高光譜衛星的投入使用,改變了因高光譜數據采集不易或數據缺乏對地質礦產勘查應用的制約,必將促進國內高光譜遙感技術的深入研究和推廣應用。
本文所述是一種技術研究和探索,目的是為了拓展和深化高光譜遙感技術的地質礦產勘查應用。該套技術方法在全巖地質信息提取和智能化地質填圖應用方面獲得了顯著效果,通過進一步完善和提高,將具有重要的推廣應用價值。