帥 爽,孟 丹,穆 超,謝 菲,胡海潮,邵澤興
(湖北省國土測繪院,湖北武漢 430010)
TM和ASTER中等分辨率遙感數據相繼應用于蝕變巖信息提取、巖性識別和影像地層劃分,在區域地質調查、礦產勘查中發揮了重大作用,遙感技術的應用已逐步走向成熟[1-5]。受其空間分辨率影響,影像像元混合嚴重,其光譜信息、結構信息不能滿足更為復雜的巖石礦物信息提取。隨著遙感技術的發展,遙感探測頻譜范圍不斷拓寬、分辨率不斷提高,可利用的數據源也不斷增多,如Worldview-2遙感數據作為新的數據源,其全色空間分辨率達0.46 m,同時多光譜具有8個波段。Landsat-8數據相對于Landsat-7新增一個深藍波段、一個短波紅外波段和一個熱紅外波段,并且波段寬度變窄。如何挖掘多源遙感數據各自的優點,將其協同應用于遙感地質領域是一項意義重大的課題,也是現代遙感找礦技術發展的必然趨勢。
在使用多種數據進行巖性識別和劃分方面,前人也做了一些研究。余海闊等為了充分利用ASTER和LandSat ETM+數據特征,將兩類數據共21個波段進行協同識別巖性,結果對巖性分類的精度比任何以其中單一數據得到效果更好[6]。但這兩種數據的大部分光譜覆蓋范圍均為重合,未達到理想的光譜協同效果。Rowan等在美國內華達州銅礦區蝕變巖信息提取中,釆用ASTER、AVIRIS和TM數據進行試驗,并評價了三種數據的礦物信息識別能力[7]。D.RoS等協同利用TM、SPOT數據,數字高程模型,地球化學數據,對愛琴海盆地島弧中的低溫熱液型金礦進行研究,用5/7、3/1等比值及主成分分析等數學運算進行蝕變帶識別[8]。這些方法并沒有深入探討多源數據各自的優勢,從光譜和空間探測能力上對數據進行擴展和提升。
本文以協同思想為指導,分析Landsat-8數據和Worldview-2數據光譜信息、空間結構信息在巖性劃分方面的作用,將不同覆蓋范圍的波段進行疊加協同,達到光譜分辨率上的協同;并將Landsat-8多光譜數據與Worldview-2數據融合,使其空間分辨率提高到2 m,最終得到兼顧光譜信息和高空間分辨率的協同數據,利用協同遙感數據可高精度地對區域巖性進行劃分。
Worldview-2和Landsat-8數據在光譜分辨率、空間分辨率、時間分辨率及輻射分辨率方面均有所差異。從地質體的屬性結構上看,影像的時間分辨率及輻射分辨率對遙感巖性的識別影響較小。因此僅對兩種數據在光譜、空間上的探測能力進行對比分析。
光譜分辨率指傳感器探測器件接收電磁波輻射時所能區分的最小波長范圍,也指傳感器在其波長范圍內所能劃分的波段量度,波段的波長范圍越小,波段越多,則波譜分辨率越高。Worldview-2、Landsat-8數據波譜范圍對比圖1可看出,在可見光—近紅外波譜范圍內,Worldview-2數據波段范圍連續分布,基本上實現該波譜范圍內波譜全覆蓋,平均波段寬度50 nm,光譜分辨率高。Landsat-8除紅綠藍波段外,還有一個深藍波段和一個近紅外波段,平均波段寬度也為50 nm,但在此波段范圍內波段覆蓋范圍沒有Worldview-2數據大,且波段不連續。
在中紅外—短波紅外波譜范圍內,Worldview-2數據沒有波段設置,Landsat-8數據在短波紅外波段設置有3個波段,平均波段寬度110 nm,光譜分辨率要高于Worldview-2數據。但Landsat-8第9波段圖像噪聲非常大,主要用于識別卷云,輻射定標后的值集中在0.011 697~0.023 395。因此,在下面的數據處理中將第9波段剔除。

圖1 Worldview-2、Landsat-8數據波譜范圍對比Fig.1 Comparison of spectrum wavelength range of Worldview-2、Landsat-8 data
空間分辨率是指圖像中可辨認的臨界物體空間幾何長度的最小極限,即對細微結構的分辨率,也可理解為影像像元的尺寸。Worldview-2數據全色空間分辨率為0.46 m,多光譜空間分辨率為2.0 m,能較好地探測地表細節信息,對地層中的節理、巖層層理及巖性單元之間的接觸關系都具有較好地表現能力。Landsat-8數據全色空間分辨率為15 m,可見光、近紅外和短波紅外分辨率均為30 m。對于在地表有一定出露規模,均質性較好的地質體,能充分發揮其識別能力;但對較大比例尺的礦物填圖和構造識別,受到混合像元影響,增大巖礦光譜研究的不確定性。因此,在地表空間結構信息挖掘上,Landsat-8數據難以監測到地表重要的細微結構信息。遙感對巖性的識別除了利用巖性單元光譜特性外,巖石的結構信息也是其重要信息。從Worldview-2數據和Landsat-8數據對空間信息的探測能力的對比上看,Worldview-2具有明顯的優勢。
綜上,Worldview-2數據與Landsat-8數據在空間探測能力和光譜探測能力上各有優勢。單獨利用某一種數據難以達到準確地進行巖石地層分類的要求。理論上將Worldview-2數據的空間探測能力優勢與Landsat-8短波紅外光譜探測能力的優勢協同,可提高遙感對巖性單元的識別能力。
研究使用的數據是新疆維吾爾自治區阿克陶縣2010年5月23日采集的Worldview-2數據和2013年6月3日采集的Landsat-8數據。
首先需要對Worldview-2多光譜數據和Landsat-8多光譜數據進行輻射定標,將影像DN值轉化輻射亮度值(Radiance),Worldview-2和Landsat-8數據輻射定標公式分布如下:

式中:Lλ為Landsat-8數據單波段的輻射亮度值;ML為單波段乘數轉換因子;AL為單波段加法轉換因子;Qcal為影像 DN 值。Gain、offset、ML、AL 從影像頭文件中讀取。
然后利用ENVI4.8平臺上的FLAASH模塊對輻射定標后的Worldview-2數據和Landsat-8數據進行大氣校正。FLAASH模塊采用MODTRAN4+輻射傳輸模型對圖像逐像元地糾正大氣中的二氧化碳、水汽及氣溶膠等吸收、散射影響。具體校正公式如式(3)。

式中:L*表示衛星傳感器上接受到的單一像元的輻射亮度值;ρ為像元在地表的反射率;ρe表示單一像元及周邊的混合像元的地表反射率;S表示大氣的球面反射率;La*為大氣輻射進入傳感器通道內的輻射亮度值;A、B為大氣條件系數。在大多數校正模型中,一般地取ρ=ρe,即是忽略“鄰近像元”效應。Landsat-8數據校正結果如圖2。

圖2 Landsat-8數據大氣校正效果Fig.2 Atmospheric correction effect of Landsat-8 data
從圖2中可以看出,大氣校正后,圖像的清晰度明顯增強,細部紋理也清晰可見。對比校正前后植被的波譜曲線,發現大氣校正后植被的波譜曲線趨于正常,大氣校正效果良好。
將經過大氣校正的Landsat-8多光譜數據與其全色數據采用Gram-Schmidt變換融合法進行融合,得到空間分辨率為15 m的數據。該方法對原始多光譜信息能較好地保留。最后將以大氣校正后的Worldview-2反射率數據為參照對Landsat-8反射率數據進行幾何校正,保證Landsat數據與Worldview-2數據在空間上配準,平均誤差保證在三個像元內。
協同學是20世紀70年代發展起來的一門新興的交叉科學,它研究一個開放系統和內部系統之間的協同合作,形成宏觀有序結構的機理和規律。1971年聯邦德國斯圖加特大學教授哈肯教授提出了協同的概念。協同理論認為,自然界是由各個層次上不同的系統所構成的統一體,在各個層次、各個系統之間既存在有相互作用和對立影響,又存在著相互制約和合化。由于系統之間的相互聯系和制約,構成了一個互為依存協調一致的具有一定發展規律的世界。對于自然界中的每個子系統是如此,對于整個由子系統構成的自然界同樣如此。自然界中各系統、層次之間的聯系促使了實物的發展和進步[9]。協同就成為連接各子系統的紐帶。因此,協同各子系統之間的內在聯系,比單一系統所起的作用更大,即表現為協同效應。
協同理論應用于遙感地質領域,是將不同衛星平臺上搭載的不同傳感器作為統一體中的不同層次,同一傳感器數據的光譜分辨率、空間分辨率又構成一個系統,它們相互聯系。由于遙感成像中瞬時視場[10]的關系,同一傳感器難以同時獲得高光譜分辨率和高空間分辨率。所以數據傳感器的光譜分辨率和空間分辨率又存在相互作用和對立的關系。
遙感巖性識別中,高空間分辨率遙感數據能較好地探測地表細節信息,對于不同類型巖石的節理、巖石層理及巖性單元之間的接觸關系都具有較好地表現能力;中等分辨率多光譜數據的短波紅外波段數據對于巖石礦物的光譜差異相對于可見光波段而言有較好地表現[11-12],對于大尺度巖性劃分、蝕變礦物的提取有一定優勢(圖3)。
因此多源遙感數據對于巖石礦物信息識別能力的發展和進步,是各個傳感器光譜探測能力和空間探測能力相互影響、相互制約“協同”推動的[13]。所以用何種方式實現多源遙感數據的空間分辨率優勢和光譜分辨率優勢的協同是遙感地質發展的趨勢。

圖3 中等分辨率多光譜數據與高空間分辨率數據巖性識別效果對比Fig.3 Comparison of multispectral data of medium resolution and data lithology of high spatial resolution
2.2.1 結構協同方案
通過對Worldview-2、Landsat-8數據空間探測能力上的分析,希望把Worldview-2數據在空間探測能力上的優勢“嫁接”到Landsat-8數據上,從而提高影像上巖性判斷的精度。首先將Landsat-8數據1~7多波段數據與其全色波段融合成15 m分辨率數據,然后再與Worldview-2多光譜數據中的某單波段數據進行融合。最終獲得具備Landsat-8數據1~7波段多光譜信息和Worldview-2數據2 m空間分辨率的空間結構協同數據。需要指出的是遙感數據融合中,如果兩種待融合數據空間分辨率比>20∶1,融合數據將損失大量信息。所以本次研究中沒有選擇0.5 m空間分辨率的Worldview-2全色數據進行融合。
2.2.2 方案優化
為使協同數據能較好地保持Landsat-8數據1~7波段多光譜信息對結構協同方案進行優化,采用具有高光譜保真度的主成分變換融合法和Gram-Schmidt變換融合法進行試融合,并對融合結果計算與原始多光譜數據的相關系數和光譜扭曲度,一般情況下,相關系數越大,光譜扭曲度越小,融合結果的光譜保真度越好(表1)。

表1 融合結果相關系數與光譜扭曲度對比表Table 1 Contrast table of correlatior coefficient of fusion results and twisting degree of spectrum
由于傳統像素級融合一對多的模式中(全色與多光譜融合),全色波段光譜響應范圍與多光譜波段光譜響應范圍不一致,傳統融合方法又沒有從物理上改變光譜不匹配的現狀,因此出現了光譜失真的問題。由于Worldview-2多光譜數據的光譜范圍不包括短波紅外波段,而短波紅外波段對于巖性劃分又至關重要,所以融合結果出現一定的光譜失真在所難免。通過計算試融合結果的相關系數和光譜扭曲度,最終選擇了光譜保真度相對較好的Gram-Schmidt變換融合法,使用了Worldview-2多光譜數據中第8波段作為高空間分辨率波段Landsat-8自融合數據進行融合。這是因為第8波段相對其他7個波段,光譜范圍最寬,相對能減小光譜響應范圍不一致帶來的光譜失真問題。
經過對融合結果的分析,這種光譜失真并不影響影像上巖性的劃分。最終結構協同效果如圖4。
圖4中可以看出,融合數據兼具了Landsat-8數據的光譜信息和Worldview-2數據的空間結構信息,達到了數據結構協同的效果。
2.3.1 光譜協同方案
通過對Worldview-2、Landsat-8數據的光譜覆蓋范圍及空間探測能力的分析,提出了將兩種數據進行協同識別巖性的構想。在光譜覆蓋范圍上,Worldview-2、Landsat-8數據光譜優勢互補,可進行光譜協同。具體方案是將Worldview-2可見光—近紅外1~8波段與Landsat-8數據短波紅外6、7波段進行疊加獲得10個波段光譜協同數據(圖5)。
2.3.2 協同方案優化
由于Landsat-8數據與Worldview-2數據傳感器及成像時間不同,經過大氣校正和高斯拉伸變換后同一地物在相同波譜范圍內的反射率值不同。如圖6,同一地物反射率,Worldview-2數據(紅色曲線)較Landsat-8數據(黑色曲線)整體偏高。

圖4 結構協同前后數據顯示效果圖Fig.4 Effect drawing of data display of structure synergy

圖5 Worldview-2、Landsat-8數據光譜協同方案圖Fig.5 Synergistic proposal map of Worldview-2,Landsat-8

圖6 Worldview-2數據、Landsat-8數據上同一地物的反射率差異Fig.6 Reflectivity difference of ground objects of WordView-2,Landsat-8 data
為了使兩種數據的反射率值匹配一致,需要確定兩種數據進行匹配的因子“S”[14]。從Worldview-2數據與Landsat-8數據的波段覆蓋范圍上分析,Worldview-2數據Band2(藍波段)波譜范圍為0.450~0.510 μm,Landsat-8數據Band2(藍波段)波譜范圍為0.450~0.515 μm。可以通過計算 Worldview-2數據 Band2均值與Landsat-8數據Band2均值的比值關系來獲得兩種數據之間的反射率值匹配因子。

計算出Worldview-2數據Band2均值為1 805.733 1,Landsat-8數據Band2均值為1 287.469 2,得到匹配系數為1.402 5。將Landsat-8自融合影像6、7波段數據乘以匹配因子1.402 5后再與Worldview-2多光譜1~8波段數據進行疊加協同,以達到光譜協同方案的優化。
通過對具體結構協同方案和光譜協同方案的實驗、研究和優化,最終確定了Worldview-2數據與Landsat-8數據協同方案,獲得具有10波段光譜信息和2 m空間分辨率的協同數據(圖7)。
分別使用Worldview-2、Landsat-8影像以及協同影像進行了目視巖性劃分,并以三組數據上烏魯克恰特組(K1w)、英吉莎群(K2E1Y)與阿爾塔什組(E1a)的顯示特征的差異說明協同數據影像地層劃分情況。研究區主要巖性單元巖性描述如表2。

表2 研究區巖性單元具體巖性表Table 2 Lithology of lithology unit in study area
圖7為烏魯克恰特組、英吉莎群與阿爾塔什組在Worldview-2、Landsat-8影像上的圖像特征,由圖可見,在Worldview-2、Landsat-8影像上烏魯克恰特組(K1w)、英吉莎群(K2E1Y)與阿爾塔什組(E1a)的界線都比較清晰。但對于英吉莎群(K2E1Y)內部庫克拜組(K2k)、烏依塔克組(K2w)、依格孜牙組(K2y)以及吐依洛克組(E1t)4個組劃分效果有所差異。

圖7 烏魯克恰特組、英吉莎群與阿爾塔什組在Worldview-2、Landsat-8影像上的圖像特征Fig.7 Image feature of Worldview-2,Landsat-8 image
圖7-a中為Landsat-8數據734波段組合影像,圖像上烏魯克恰特組(K1w)石英砂巖呈深紫色,阿爾塔什組(E1a)白色塊狀硬石膏巖呈淺藍色,英吉莎群(K2E1Y)整體呈綠色或淺紅色。三個巖石單元走向近南北向,接觸關系均為整合接觸。英吉莎群(K2E1Y)內部庫克拜組(K2k)灰綠色頁巖與灰黃色頁巖互層的巖性在Landsat-8影像上表現為淺綠色中夾雜著黃色。其上伏的烏依塔克組(K2w)的長石砂巖與泥巖互層與依格孜牙組(K2y)的灰色泥晶灰巖在影像均顯示為深綠色,難以劃分。Landsat-8影像上,英吉莎群(K2E1Y)上部的吐依洛克組(E1t)的紫紅色厚層泥巖呈紅色、暗紅色與周圍巖性單元區分比較明顯。
圖7-b中為Worldview-2數據432波段組合影像,圖像上各巖性單元的顏色差距并不明顯。烏魯克恰特組(K1w)石英砂巖呈暗紅色,英吉莎群(K2E1Y)與阿爾塔什組(E1a)顏色比較接近,色調都較淺,但英吉莎群(K2E1Y)層理比較發育,阿爾塔什組(E1a)則是塊狀的硬石膏巖,以此為劃分標準,也比較容易將兩個巖性單元分開。英吉莎群(K2E1Y)下部的庫克拜組(K2k)和烏依塔克組(K2w)在Worldview-2影像上色調比較接近,結構上層理并不發育,所以難以從影像上進行劃分。依格孜牙組(K2y)的灰色泥晶灰巖水平層理比較發育,在Worldview-2影像上表現為一組較密集的平行紋理。以此為特征比較容易與其下伏烏依塔克組(K2w)區分。而依格孜牙組(K2y)與吐依洛克組(E1t)顏色上均呈深灰色,內部紋理也比較接近,所以比較難劃分。
圖8為烏魯克恰特組(K1w)、英吉莎群(K2E1Y)與阿爾塔什組(E1a)在協同影像上的表現。可以看出三個組之間的接觸關系都非常清晰。

圖8 烏魯克恰特組、英吉莎群與阿爾塔什組在協同影像上的圖像特征Fig.8 Image feature of synergistic image
烏魯克恰特組(K1w)整體呈紫色,與上伏英吉莎群(K2E1Y)為整合接觸關系[15]。從下部到上部,紫色逐漸加深,這與該區域1∶25萬地質調查報告上表述的從棕紅色粉砂質泥巖—淺褐色厚層狀含泥礫鈣質細粒長石石英砂巖夾細礫巖、棕紅色極薄層狀長石石英粉砂巖的巖性相符。另外紫色加深,與該組不同層位石英含量有關,石英礦物的特征吸收譜帶在 2.195 μm[16],對應協同數據的第10波段,也就是影像上假彩色合成的紅通道,所以巖性上石英含量增高,圖像上表現為紫色的加深。結構上從下部向上部巖層變薄,水平層理在一些層位比較發育。
英吉莎群(K2E1Y)在影像上表現為不同顏色的幾個條帶,與其上伏阿爾塔什組(E1a)為整合接觸關系。英吉莎群內部分為庫克拜組(K2k)、烏依塔克組(K2w)、依格孜牙組(K2y)以及吐依洛克組(E1t)四個組,所以在影像上表現為4個顏色、結構都不相同的條帶。庫克拜組(K2k)在英吉莎群中色調比較淺,與其下伏的烏魯克恰特組容易區分。灰綠色、灰黃色頁巖互層巖性在影像上表現為白黃色相間的條紋狀紋理,同時水平層理比較發育。烏依塔克組(K2w)巖性為褐紅色中層狀細粒巖屑長石砂巖與泥巖互層,在影像上表現為比較均一的淺紫色,層理也比較發育。色調上的差異也是它與下伏的庫克拜組比較容易區分。依格孜牙組(K2y)主體為灰色泥晶灰巖,在影像上表現為整體深灰色。另外該組灰巖中水平層理十分發育,在影像上表現為間距相等的一組平行紋理。依據特征可將其與下伏烏依塔克組(K2w)進行劃分。吐依洛克組(E1t)在影像上表現為比較特殊的暗紅色,其巖性為棕紅、紫紅色厚層泥巖。依據顏色上的特征,使其比較容易與下伏的依格孜牙組區分。
阿爾塔什組(E1a),在影像上為比較均勻的淺藍色。該組巖性為白色塊狀硬石膏巖,所以在影像上層理并不發育,并且顏色比較均一。結合地質調查資料和影像特征可以判定其與下伏的英吉莎群(K2E1Y)為整合接觸關系。
綜上,Landsat-8數據更能體現巖性單元間的光譜差異,在英吉莎群(K2E1Y)內依格孜牙組(K2y)和吐依洛克組(E1t)的劃分上優勢明顯。Worldview-2數據則在反應巖性單元間結構信息差異上有優勢,如沉積紋理、韻律。在英吉莎群(K2E1Y)內烏依塔克組(K2w)和依格孜牙組(K2y)的劃分上具有優勢。協同數據集合Worldview-2數據與Landsat-8數據的空間結構信息優勢和光譜信息優勢,在實驗區的巖性劃分效果上優于單獨使用Worldview-2數據或Landsat-8數據。
從協同理論和Worldview-2、Landsat-8數據的特點出發,將協同理論應用于多源遙感數據遙感地質領域,以提高遙感巖性識別精度。
(1)分析了Worldview-2、Landsat-8數據對于遙感巖性識別的光譜探測能力和空間探測能力,并從兩種數據光譜和空間探測能力的優劣出發,將協同理論引入遙感地質領域,闡述了協同運用多源遙感數據進行巖礦信息提取的理論基礎。提出了相應的結構協同方案和光譜協同方案,通過對方案的優化,確定了最終的協同方案。
(2)分別使用Worldview-2、Landsat-8和協同數據影像對同一組巖性單元進行劃分,并對比劃分結果,表明Worldview-2、Landsat-8數據對于不同巖性單元的劃分各有優劣;協同數據的巖性劃分能力要優于單獨使用Worldview-2和Landsat-8數據。
本次研究針對巖層出露較好的沉積巖區,驗證了多源遙感數據協同進行巖性劃分方法的可行性,研究結果表明數據協同有效地提高了遙感數據的巖性劃分能力,然而對于巖漿巖區和變質巖區的巖性劃分能力有待進一步研究。
[1] 張玉君,楊建民,陳薇.ETM+(TM)蝕變遙感異常提取方法研究與應用—地質依據和波譜前提[J].國土資源遙感,2002(4):30-37.
[2] Lawrence C.Rowan.Lithologic mapping in the Mountain Pass,California area using Advanced Spaceborne Thermal Emission and Reflection Radiometer(ASTER)data[J].Remote Sensing of Environment,2003(84):350-366.
[3] Yamaguchi Y.Spectral indices for lithologic discrimination and mapping by using the ASTER SWIR bands[J].INT.J.REMOTESENSING,2003(24):4311-4323.
[4] 李培軍.用ASETR圖像和地統計學紋理進行巖性分類[J].礦物巖石,2004(3):116-120.
[5] XianFeng Zhang.Lithologic and mineral imformation extraction for gold exploration using Aster data in the south chocolate mountains(California)[J].Journal of Photogrammetry & Remote Sensing,2007(62):271-282.
[6] 余海闊,李培軍.運用LANDSAT ETM+和ASTER數據進行巖性分類[J].巖石學報,2010(1):345-350.
[7] Rowan L C,Hook S,Abrams MJ,et al.Mapping hydrothermally altered rocks at CUPRITE,NEVADA,using the Advanced Spaceborne Thermal Emission and Reflection Radiometer(ASTER),a new satelliteimaging system[J].Economic Geology,2003(98):1019-1027.
[8] ROKOS D,ARGIALAS D,MAVRANTZA R,et al.Structural analysis for gold mineralization using remote sensing and geochemical techniques in a GIS environment:island of Lesvos,Hellas[J].Natural Resources Research,2000,9(4):277-293.
[9] 郭治安,沈小峰.協同論[M].太原:山西經濟出版社,1991:19-21.
[10] 馬艷華.高空間分辨率和高光譜分辨率遙感圖像的融合[J].紅外,2003(10):11-16.
[11] GR Hunt.Spectral signatures of particulate minerals in the visible and nearinfrared[J].Geophysics,1977(3):501-513.
[12] Crosta A,Moore J.Enhancement of Landsat Thematic Mapper imagery for residual soil mapping in SW Minas Gerais State,Brazil-A prospecting case history in greenstone belt terrain[Z].1990:1173-1187.
[13] Murphy S W,Wright R,Oppenheimer C.MODIS and ASTER synergy for characterizing thermal volcanic activity[J].Remote Sensing of Environment,2013,131(4):195-205.
[14] 余健.ASTER與worldview-2影像協同應用提取巖性信息方法研究[D].武漢:中國地質大學(武漢),2012.
[15] 盧書煒,杜鳳軍,方懷賓,等.新疆 1∶25萬英吉沙縣幅(J43C002003)區域地質調查報告[R].鄭州:河南省地質調查院,2004.
[16] Hunt,R G.Near-infrared(1.3 ~ 2.4 μm)spectra of alteration minerals;potential for use in remote sensing[J].Geophysics,1979,44(12):1974-1986.