趙恒謙，高尉，楊天藝，劉帥，趙學勝

1.中國礦業(yè)大學(北京)地球科學與測繪工程學院，北京 100083；2.煤炭資源與安全開采國家重點實驗室，北京 100083

蝕變信息提取是遙感在地質(zhì)領域最成功也是最重要的應用之一[1]。根據(jù)蝕變巖和正常巖石在遙感影像中特征波段的光譜響應差異，可以提取對找礦有重要指示意義的蝕變礦物信息，如鐵染蝕變、羥基蝕變、碳酸巖蝕變等[2-6]。但是，已有研究主要在裸露或半裸露巖層區(qū)域開展[7]。在植被覆蓋區(qū)域遙感影像中的蝕變礦物信息非常微弱，植被信息的干擾在很大程度上影響遙感蝕變信息的提取精度，制約了遙感找礦技術的進一步應用推廣[8]。

為了解決植被信息對遙感蝕變信息提取的影響，許多學者以ETM+/TM和ASTER為數(shù)據(jù)源，在植被覆蓋地區(qū)開展了多種遙感蝕變信息提取方法，如最佳波段+主成分分析、無損線性拉伸+抑制干擾因素+波段比值法、主成分分析+密度分割、MNF變換+礦物標識等方法[9-11]。然而，這些研究方法都是針對中低植被覆蓋區(qū)，對于高植被覆蓋區(qū)的蝕變信息提取效果不明顯。在高植被覆蓋區(qū)，遙感圖像分析方法提取的巖礦信息多是小目標、弱信息[12]。因此，又有學者在分析礦物、植被、土壤地面光譜的基礎上，采用Hyperion高光譜衛(wèi)星數(shù)據(jù)對高植被覆蓋區(qū)進行了植被抑制和遙感蝕變信息提取研究，證實基于差分的植被抑制方法可以有效地抑制植被信息并提取巖礦弱信息[13]。然而，Hyperion高光譜衛(wèi)星數(shù)據(jù)幅寬較窄，無法獲取大面積區(qū)域的蝕變礦物分布情況。

目前主流的遙感衛(wèi)星傳感器是多光譜傳感器，多光譜數(shù)據(jù)空間分辨率和時間分辨率都較高，幅寬面積大，數(shù)據(jù)質(zhì)量好且數(shù)據(jù)源眾多[14]，其中Sentinel-2數(shù)據(jù)是近年來比較新且有代表性的衛(wèi)星多光譜數(shù)據(jù)。Sentinel-2由兩顆相同的衛(wèi)星哨兵2號A(Sentinel-2A)與B(Sentinel-2B)組成，具有多光譜、重訪周期短(雙星組網(wǎng)5 d)、空間分辨率高、數(shù)據(jù)幅寬范圍大及免費獲取等優(yōu)點，已經(jīng)在裸露巖層區(qū)蝕變礦物信息提取研究中得到應用[15]。然而，在高植被覆蓋區(qū)，如何去除植被干擾獲取蝕變礦物信息仍是一個難題。因此，本文采用Sentinel-2A作為衛(wèi)星多光譜數(shù)據(jù)源，以黑龍江呼瑪?shù)貐^(qū)作為實驗區(qū)，通過植被與礦物反射率做差分計算的方法來抑制植被信息，利用主成分分析方法提取大面積高植被覆蓋區(qū)的蝕變礦物信息。

1 實驗區(qū)和數(shù)據(jù)處理

本研究選取黑龍江呼瑪?shù)貐^(qū)作為研究區(qū)。該區(qū)域位于E126°—127°/N50°—52°的大興安嶺東麓，地勢呈西北高(丘陵區(qū))、東南低(平原區(qū))，平均海拔350 m。研究區(qū)植被以林木為主，植被覆蓋率達95%以上；土地肥沃，有棕壤、黑土、草甸沼澤土和草甸土等豐富的土壤資源；同時，礦產(chǎn)資源也十分豐富，盛產(chǎn)金、鐵、煤、石墨、石英、云母以及高嶺土等。

Sentinel-2A衛(wèi)星于2015年6月23日發(fā)射升空，設計壽命為7年，軌道高度786 km[16]。Sentinel-2A衛(wèi)星上攜帶的高分辨率多光譜成像儀(MSI)，覆蓋可見光、近紅外和短波紅外13個光譜段(表1)，衛(wèi)星重訪周期為10 d，成像幅寬達290 km。本研究使用的影像獲取時間為2017年10月26日，產(chǎn)品等級為已幾何精校正的L1C級，下載自歐洲太空局官網(wǎng)(https：//scihub.copernicus.eu/dhus/#/home)。采用整幅Sentinel-2A影像作為高植被覆蓋區(qū)總體研究范圍(約11 100 km2，圖1)，其中Hyperion數(shù)據(jù)所覆蓋的重合區(qū)域為重點研究區(qū)[13](約337 km2，圖1黃框區(qū))。

圖1 本研究高植被覆蓋區(qū)范圍(Sentinel-2A數(shù)據(jù))Fig.1 High vegetation coverage area in this study (Sentinel-2A data)

為了保證蝕變礦物信息提取的精度，使用ENVI 5.5軟件對原始數(shù)據(jù)進行預處理，包括輻射定標、大氣校正、波段融合以及實驗范圍裁剪。

表1 Sentinel-2多光譜傳感器波段信息

2 蝕變信息提取方法

對于已知礦物類別的地區(qū)，可采用光譜角匹配和支持向量機方法直接提取蝕變信息。但對基礎地質(zhì)工作少、地面作業(yè)困難且有植被因素影響的地區(qū)，更多是使用比值與主成分分析方法來進行蝕變信息的識別與提取。

將植被抑制后的影像結合主成分分析作散點圖，提取蝕變礦物信息。主成分變換是目前應用最多的特征壓縮和特征選擇的經(jīng)典方法之一[17]。對于高植被覆蓋區(qū)蝕變礦物信息提取，前提工作是去除植被干擾。劉彥麗[13]通過對比分析植被與土壤/礦物光譜曲線特征差異，提出了對波段做差分運算的差分植被抑制方法。另外，ENVI軟件中也提供了一個簡單快速的植被抑制工具Vegetation Suppression。

2.1 植被抑制方法

2.1.1 多光譜差分植被抑制法

差分植被抑制法是指結合特定傳感器波段設置，選擇植被反射率相同而土壤/礦物反射率差異較大的波段組合，將波段組合的反射率做差分計算，從而對影像中的植被信息進行抑制，為通過主成分分析散點圖提取蝕變信息奠定基礎。在高光譜植被差分抑制中，由于傳感器具有較高光譜分辨率，很多波段組合都符合差分植被抑制法的要求，如1 336～742 nm、2 072～548 nm、2 314～701 nm、1 699～721 nm、2 203～681 nm、2 183～671 nm等。由于波段組合較多，只選取部分最優(yōu)差分波段組合進行展示，如圖2(a)所示。

圖2 植被抑制差分波段選擇示意圖Fig.2 The schematic diagram of the selection of differential bands for vegetation suppression

值得注意的是，經(jīng)差分植被抑制法處理之后，還需要對差分后的波段進行主成分分析并繪制散點圖，才能提取蝕變礦物信息，因此至少需要3對差分波段組合才能滿足多光譜差分植被抑制法的波段設置要求。對于最優(yōu)化差分波段組合來講，植被光譜在一對組合內(nèi)的兩個波段中應具有相近的反射率。將實驗區(qū)Sentinel-2多光譜傳感器的植被和土壤/礦物光譜曲線進行對比，如圖2(b)所示?？梢钥吹?，Sentinel-2傳感器具備3對波段組合滿足最優(yōu)化差分波段設置條件，從而可以通過MDVS處理得到3對差分波段(664～496 nm，2 202～560 nm，945～783 nm)。相比之下，Landsat 7/8、ASTER等其他多光譜衛(wèi)星傳感器不滿足最優(yōu)化差分波段設置要求(表2)。在實際應用中，也可以結合不同傳感器波段設置，選取反射率接近的波段組合進行差分處理。本研究中，采用波段設置條件最優(yōu)的Sentinel-2傳感器數(shù)據(jù)進行分析和處理。

表2 常用多光譜傳感器植被抑制差分波段選取結果

2.1.2 ENVI植被抑制工具

ENVI軟件中提供了植被抑制工具Vegetation Suppression。該工具利用遙感影像的紅光波段和近紅外波段消除或削弱多光譜和高光譜影像中的植被光譜信息，對圖像進行植被抑制[18]。利用植被抑制工具Vegetation Suppression對影像進行植被抑制后，波段數(shù)不會改變。研究使用的植被抑制工具是ENVI 5.5 版本中的Vegetation Suppression。

2.2 蝕變信息提取方法

對MDVS法植被抑制后的結果再進行主成分變換，分離植被與蝕變信息。選擇貢獻值符號相反且絕對值最大的兩個主成分分量，做2D散點圖并圈出異常范圍。散點圖中所圈定的異常范圍(蝕變信息)分布在遙感影像中，經(jīng)過波譜分析確定蝕變礦物類型。

對于ENVI植被抑制工具處理得到的結果，羥基蝕變礦物在Sentinel-2A影像的波長1 600 nm附近(Band 11)有強反射峰，而在波長2 200 nm附近(Band 12)有強吸收谷[19]。因此，選擇Sentinel-2A數(shù)據(jù)的Band 2、4、11、12這四個波段組合進行主成分分析，提取羥基蝕變信息。Sentinel-2A數(shù)據(jù)中鐵離子的光譜特征表現(xiàn)為：在Band 3的反射率最低，在Band 4的反射率相對Band 3高，在Band 8a的反射率下降且呈吸收谷特征，在Band 11的反射率快速上升形成反射峰，隨后反射率又逐漸下降[20]。因此，提取鐵染蝕變異常信息選取Sentinel-2A 數(shù)據(jù)的Band 3、4、8a、11這四個波段進行主成分分析。主成分變換后，同樣使用貢獻值符號相反且絕對值最大的兩個主分量做散點圖，最后圈定出蝕變礦物信息。

2.3 蝕變信息種類匹配及精度驗證

蝕變信息提取后，對異常光譜進行匹配。光譜匹配的常用方法有光譜角匹配法(SAM)、光譜特征擬合法(SFF)以及二維編碼法等。本研究使用SFF法，該方法的特點是對選定的匹配光譜特征區(qū)間內(nèi)的光譜整體波形、吸收位置及吸收深度作對比，從而判斷光譜的匹配程度[21]。經(jīng)SFF光譜匹配后，識別出各異常代表的蝕變類型。

重點研究區(qū)地面驗證點的礦物類型，是通過在研究區(qū)采集巖石樣品、實驗室切片、鏡下分析鑒定而得到的(表3)。將影像上提取的蝕變信息光譜與USGS光譜庫中光譜進行匹配，結合表3已知驗證點位的礦物種類進行提取精度的驗證。重點研究區(qū)的地面驗證點位分布如圖3所示。

表3 重點研究區(qū)地面驗證點礦物種類[13]

圖3 重點研究區(qū)地面驗證點位分布圖Fig.3 Distribution map of field verification points in key experimental areas

3 數(shù)據(jù)處理結果與分析

3.1 植被抑制結果和分析

多光譜差分植被抑制共選擇3對波段，將包含植被光譜特征的綠峰和紅谷附近對應的2對波段(2 202～560 nm)和(664～496 nm)做散點圖；選擇Band 12(2 202 nm)和Band4(664 nm)做ENVI植被抑制工具的植被抑制后散點圖；將兩幅散點圖與原始影像散點圖做直觀對比，分析MDVS法與ENVI植被抑制工具法的植被抑制效果(圖4)。

紅圈為蝕變異常信息圖4 植被抑制前后散點圖對比Fig.4 Comparison of scatter plots before and after vegetation suppression

結果表明，原始影像的散點圖呈扁橢圓形狀，植被信息和礦物蝕變信息混雜一起，總體表現(xiàn)為植被背景信息。而植被反射率相同且土壤/礦物反射率差異大的兩個波段做差分計算后，所求出的新波段間散點圖相比于原始影像散點圖，其異常信息在數(shù)據(jù)空間分布上突出明顯，范圍更大(圖中紅色虛線所圈范圍)。ENVI植被抑制工具進行植被抑制的散點圖相比原始數(shù)據(jù)，在數(shù)據(jù)空間分布上變化較小，但也有一小部分的異常范圍，抑制植被有一定的效果。與ENVI植被抑制工具相比，MDVS方法得到的散點圖中的異常信息更突出且更多，植被抑制效果更好。

3.2 蝕變信息提取結果

3.2.1 MDVS蝕變信息提取結果

利用MDVS方法將Sentinel-2數(shù)據(jù)選出的3對波段做差分，并計算各差分波段的主成分特征值，結果見表4。

表4 MDVS處理后差分波段的主成分特征值Tab.4 Principal component eigenvalues of the differential bands after MDVS processing

從表4可以看出，主成分1(PC1)與主成分2(PC2)特征值絕對值較大，且差分波段1和差分波段3特征值正負號相反，能夠較好地反映不同差分波段的信息，因此適合選擇這兩個主成分進行散點圖分析。

在理想狀況下，散點圖形狀呈三角形，拐點分布在其3個頂點。在實際情況中，散點圖形狀并非呈三角形，其拐點為散點圖周圍的逸散部分，而蝕變異常信息往往分布在這些逸散部分。因此，在PC1和PC2的二維散點圖中，先圈定外層的信息設為一級異常，再圈定內(nèi)層的信息設為二級異常，分別求一級異常和二級異常的平均光譜，然后對這兩個光譜進行光譜匹配，根據(jù)匹配分數(shù)確定是否為同一級異常，最后逐漸往更內(nèi)層圈定。這樣可將不同位置的逸散部位圈定出來，并使用不同顏色標注，即可得到不同類型的蝕變異常信息，如圖5所示。

紅、綠、藍3種顏色代表光譜不同的異常信息圖5 MDVS結果主成分散點圖Fig.5 The scatter plot of principal components of MDVS results

將從散點圖圈定的不同蝕變異常信息標記到遙感影像中，可以得到蝕變信息分布圖(圖6)。蝕變信息提取結果圖中的紅、綠、藍顏色，與相應散點圖中的異常信息顏色分別對應。為了進一步確定不同顏色區(qū)域代表的蝕變異常種類，需要對各個區(qū)域樣本光譜與已有光譜庫進行匹配，對其礦物類型進行識別。

紅、綠、藍3種顏色代表光譜不同的異常信息圖6 重點研究區(qū)MDVS蝕變信息提取結果Fig.6 Results of alteration information extraction using MDVS method in the key experimental area

通過與地面驗證信息進行對比，發(fā)現(xiàn)地面驗證點位1、2、3恰好分別位于蝕變信息分布圖中紅、綠、藍色3種不同異常信息范圍內(nèi)。選取驗證點位的Sentinel-2A光譜與USGS礦物光譜庫中的所有光譜重采樣到相同分辨率的光譜，再進行SFF法光譜匹配，確定了3種顏色的異常信息所對應的蝕變礦物類型，結果如圖7所示。從圖7可以看到，驗證點衛(wèi)星數(shù)據(jù)光譜匹配結果和USGS光譜庫重采樣光譜特征一致，匹配效果非常好。將光譜匹配結果和表2中地面考察結果進行對比發(fā)現(xiàn)，點1、點2提取的礦物種類識別準確，點3提取結果相對較差，這可能與衛(wèi)星空間分辨率和地面考察在尺度上有很大差異有關。

圖7 驗證點光譜與USGS光譜庫重采樣 光譜匹配結果對比Fig.7 Spectra at field verification points and the matched spectra in the resampled USGS spectral library

黏土類礦物蝕變含有羥基離子，鐵染蝕變含有鐵離子，這些離子在遙感圖像上表現(xiàn)出獨特的光譜特征，從而使蝕變礦物區(qū)別于非蝕變礦物。在本文的重點研究區(qū)圍巖蝕變類型中，絹云母化和黏土化(高嶺土化)都含有羥基團，可歸類為羥基蝕變；而角閃石是鐵的還原化物，歸類為鐵染蝕變。圖6中提取結果的綠色和藍色異常信息歸為羥基蝕變，紅色異常信息為鐵染蝕變。因此，在多光譜影像上最終的蝕變信息可歸為羥基蝕變和鐵染蝕變兩大類，如圖8所示。

圖8 重點研究區(qū)MDVS法蝕變信息提取結果Fig.8 Results of alteration information extraction by MDVS method in key experimental areas

3.2.2 ENVI植被抑制蝕變信息提取結果

根據(jù)羥基蝕變和鐵染蝕變礦物在Sentinel-2A影像上不同波段具有不同的光譜吸收特征，選擇具有對應波譜特征的波段做主成分變換，其結果見表5和表6。

表5 羥基蝕變提取波段的主成分特征貢獻值Tab.5 Principal component eigenvectors of hydroxyl alteration extraction bands

表6 鐵染蝕變提取波段的主成分特征貢獻值

從表5可以看出，PC2滿足羥基蝕變信息主成分分量選擇標準(Band 11波段和Band 12波段特征值正負號相反，且貢獻率最大)。對PC2進行密度分割提取羥基蝕變信息，如圖9(a)圖像中的藍色部分。

同理，根據(jù)表6可以確定鐵染蝕變信息主成分分量，只有PC2滿足鐵染蝕變信息主成分分量選擇標準(Band 8a波段與Band 11波段特征值正負號相反，且貢獻率最大)。對PC2進行密度分割提取鐵染蝕變信息，如圖9(b)圖像中的紅色區(qū)域。

圖9 重點研究區(qū)ENVI植被抑制工具蝕變信息提取結果Fig.9 Results of alteration information extraction by ENVI vegetation suppression tool method in key experimental areas

為了更直觀對比兩種不同方法提取的重點研究區(qū)蝕變礦物信息結果，現(xiàn)將提取的羥基蝕變和鐵染蝕變信息疊加進行考察，如圖10所示。

圖10 重點研究區(qū)蝕變信息疊加結果Fig.10 Overlapping results of alteration information in the key experimental area

利用ENVI植被抑制工具Vegetation Suppression，根據(jù)影像的紅波段和近紅外波段信息對Sentinel-2A影像進行植被抑制，然后選取特征波段做主成分變換，提取蝕變信息，如圖10(b)所示?？梢钥吹?，蝕變信息的分布位置與差分植被抑制后提取的蝕變信息分布較一致，點1、2和3也包含差分植被抑制法得到的蝕變信息，但信息范圍較小。用Sentinel-2A提取的蝕變信息與已知的實驗區(qū)Hyperion影像提取結果(圖11)作對比，可以看出，蝕變信息大部分被提取出來，并且MDVS法進行植被抑制后提取的蝕變范圍和位置更準確，信息量也豐富。

圖11 重點研究區(qū)Hyperion影像礦物 信息提取結果[13]Fig.11 Extraction of mineral information from Hyperion image in key experimental areas[13]

為了進一步驗證Sentinel-2A數(shù)據(jù)和MDVS植被抑制方法的適用性，對包括實驗區(qū)在內(nèi)的更大范圍研究區(qū)影像進行了蝕變礦物信息提取，結果如圖12所示。從圖12可以看到，利用反射率波段做差分的MDVS植被抑制方法，能夠在比Hyperion數(shù)據(jù)幅寬更大的多光譜影像中進行更大空間范圍的蝕變信息提取，且結果與高光譜數(shù)據(jù)相一致，這也證明了該方法在多光譜影像上的適用性和Sentinel-2數(shù)據(jù)提取蝕變信息在幅寬上的優(yōu)勢。

黃色虛線框內(nèi)范圍為重點研究區(qū)圖12 高植被覆蓋區(qū)MDVS蝕變信息提取疊加結果Fig.12 Overlapping result of the alteration information extraction in the high vegetation coverage area using MDVS

4 結 論

針對大范圍高植被覆蓋區(qū)遙感找礦面臨的技術難題，在差分法植被抑制原理的基礎上，結合多光譜傳感器的波段設置進行差分波段選擇，取得了以下成果：

(1) 基于Sentinel-2傳感器波段設置特點，提出了多光譜差分植被抑制法(MDVS)，為實現(xiàn)高植被覆蓋區(qū)蝕變礦物信息提取由高光譜數(shù)據(jù)到多光譜數(shù)據(jù)處理推廣，提供了重要的技術前提。

(2) 典型實驗區(qū)Sentinel-2A影像數(shù)據(jù)分析結果表明，MDVS法可以取得比已有方法更好的植被抑制效果，且基于該方法處理后的數(shù)據(jù)能夠有效提取蝕變信息。

(3) 獲取了大范圍高植被覆蓋區(qū)的蝕變信息分布圖，實現(xiàn)了高植被覆蓋區(qū)蝕變礦物信息提取推廣到多光譜傳感器這一重要突破，對于遙感找礦技術的發(fā)展具有重要意義。

此次實驗中的高植被覆蓋區(qū)影像的蝕變信息混合了大量植被信息，光譜波形和純凈礦物有很大差異。如何在提取蝕變異常大類的基礎上進一步匹配識別具體礦物種類，仍有待進一步研究。另外，村莊、道路等信息在散點圖中也容易被識別為異常信息，如何將其與蝕變信息進行區(qū)分需要展開更多實驗分析?？傊?，在MDVS方法的基礎上，如何進一步改善高植被覆蓋區(qū)蝕變信息提取的精度，將是下一步研究工作的重點。