內(nèi)蒙古本巴圖鈾資源勘查遠(yuǎn)景區(qū)航空高光譜數(shù)據(jù)不同大氣校正方法評(píng)價(jià)及應(yīng)用

張川，葉發(fā)旺，童勤龍，郭幫杰，李新春，淦清清

（核工業(yè)北京地質(zhì)研究院 遙感信息與圖像分析技術(shù)國(guó)家級(jí)重點(diǎn)實(shí)驗(yàn)室，北京 100029）

高光譜遙感成像技術(shù)的光譜分辨率一般達(dá)到5～15 nm，能夠區(qū)分具有診斷性光譜特征的地表地物。在地質(zhì)勘探領(lǐng)域，利用高光譜圖像光譜特征與標(biāo)準(zhǔn)光譜數(shù)據(jù)庫(kù)中光譜的相似性進(jìn)行巖礦識(shí)別和精細(xì)分類(lèi)，在地質(zhì)礦產(chǎn)填圖方面具有廣泛的應(yīng)用前景［1-2］。與衛(wèi)星高光譜相比，航空高光譜兼具高光譜和高空間分辨率的雙重優(yōu)勢(shì)，單個(gè)像素的地物混合效應(yīng)大為降低，對(duì)提取多種類(lèi)型的巖石礦物更為有利，尤其是與礦產(chǎn)密切相關(guān)的蝕變礦物。當(dāng)前，遙感地質(zhì)調(diào)查一般先采用衛(wèi)星影像進(jìn)行初步分析和潛力評(píng)價(jià)，然后對(duì)篩選出的成礦有利區(qū)帶進(jìn)行航空高光譜調(diào)查，以獲取更為精細(xì)的巖石礦物信息。

在光學(xué)遙感常用的0.4～2.5 μm 波長(zhǎng)范圍，成像光譜儀測(cè)得的地面目標(biāo)的總輻射亮度并不是地表真實(shí)反射率的反映，其中包含了由大氣吸收，尤其是散射作用造成的輻射量誤差。大氣校正是消除這些由大氣影響所造成的輻射誤差，反演地物真實(shí)反射率的過(guò)程。因此，大氣校正是還原或增強(qiáng)航空高光譜影像像元地物診斷性光譜特征的一項(xiàng)關(guān)鍵技術(shù)，大氣校正的質(zhì)量將直接影響后續(xù)基于航空高光譜反射率數(shù)據(jù)的巖礦填圖精度。

自20 世紀(jì)70 年代起，眾多學(xué)者在遙感大氣校正方面開(kāi)展了研究工作，產(chǎn)生了多種大氣校正理論方法，并開(kāi)發(fā)了相應(yīng)的軟件［3］。其中，部分方法僅適用于低分辨率的衛(wèi)星影像，對(duì)反射率光譜重建的精度要求亦不嚴(yán)格。例如，暗目標(biāo)法，用于消除程輻射［4］；去霧算法，通常視為一種相對(duì)大氣校正［5］。由于遙感礦物填圖需要以較為嚴(yán)格的反射率反演為前提，因此，以上算法對(duì)地質(zhì)領(lǐng)域的航空高光譜影像礦物填圖并不適用。為此，部分學(xué)者針對(duì)航空高光譜的大氣校正開(kāi)發(fā)了多種方法。本研究以?xún)?nèi)蒙古巴丹吉林盆地東部的本巴圖鈾資源勘查遠(yuǎn)景區(qū)為試驗(yàn)區(qū)，針對(duì)礦物填圖應(yīng)用需求，開(kāi)展不同大氣校正模型對(duì)航空高光譜影像反射光譜重建的評(píng)價(jià)研究，通過(guò)典型礦物反射光譜的綜合對(duì)比，評(píng)價(jià)不同模型的校正效果，并用以指導(dǎo)試驗(yàn)區(qū)航空高光譜影像大氣校正和礦物填圖。

1 大氣校正模型

高光譜遙感圖像大氣校正模型方法主要包括圖像統(tǒng)計(jì)學(xué)模型法和大氣輻射傳輸模型法［6-8］。圖像統(tǒng)計(jì)學(xué)模型法包括平場(chǎng)域法（Flat Field Transfer，F(xiàn)FT）、內(nèi)部平均法（Internal Average Relative Reflectance，IARR）、對(duì)數(shù)殘差法（Log Residuals，LR）和經(jīng)驗(yàn)線(xiàn)性法（Empirical Line Calibration，ELC）。前3 種方法不需要進(jìn)行實(shí)際地面光譜及大氣環(huán)境參數(shù)的測(cè)量，而是直接從圖像特征本身出發(fā)消除大氣影響，進(jìn)行反射率反演，基本屬于數(shù)據(jù)歸一化的范疇，雖然簡(jiǎn)便快捷，但只得到了反射率的相對(duì)值。航空高光譜遙感巖礦信息提取基于真實(shí)的巖礦反射光譜特征，因此，一般采用能夠獲得絕對(duì)反射率的經(jīng)驗(yàn)線(xiàn)性法和大氣輻射傳輸模型法進(jìn)行處理。對(duì)于大氣輻射傳輸模型法，MODTRAN 和6S 是目前最常用的輻射傳輸模型，實(shí)測(cè)大氣參數(shù)可用于輔助實(shí)現(xiàn)精確的輻射傳輸計(jì)算。目前，在缺少實(shí)際地面測(cè)量光譜時(shí)最為常用的 FLAASH（Fast-Line-of-sight Atmospheric Analysis of Spectral Hypercubes）模型法，是在MODTRAN 的基礎(chǔ)上發(fā)展的。6S 模型法和經(jīng)過(guò)輻射傳輸模型簡(jiǎn)化的快速大氣校正法（Quick Atmospheric Correction，QUAC），亦常被用于處理高光譜遙感圖像。本文的方法試驗(yàn)包括FLAASH 模型法、6S 模型法、快速大氣校正模型法和基于實(shí)際地面測(cè)量光譜的經(jīng)驗(yàn)線(xiàn)性法。

1.1 FLAASH 模型法

FLAASH 模型是由美國(guó)光譜科學(xué)研究所（Spectral Sciences Inc.）和空軍研究實(shí)驗(yàn)室（Air Force Research Laboratory）共同研制開(kāi)發(fā)的，F(xiàn)LAASH 模型法采用MODTRAN 4+輻射傳輸模型的代碼，基于查找表（Look-up Table）和插值方法計(jì)算，能夠精確補(bǔ)償大氣影響，是目前普遍認(rèn)為精度較高的大氣輻射傳輸模型方法。該方法是基于像素級(jí)的校正，集成在ENVI 遙感處理軟件中，能對(duì)400～2 500 nm 波長(zhǎng)范圍內(nèi)的遙感影像進(jìn)行大氣校正［9］。

FLAASH 模型法不是在預(yù)先計(jì)算好的模型數(shù)據(jù)庫(kù)中加入輻射傳輸參數(shù)來(lái)進(jìn)行大氣校正，而是直接結(jié)合MODTRAN 4+的大氣輻射傳輸源碼，通過(guò)特殊譜段反演氣溶膠厚度和水汽含量，進(jìn)而為每一幅影像生成一個(gè)唯一的MODTRAN 解決方案，同時(shí)，標(biāo)準(zhǔn)的MODTRAN大氣模型和氣溶膠類(lèi)型也可被直接使用。

1.2 6S 模型法

6S 模型是法國(guó)里爾科技大學(xué)大氣光學(xué)實(shí)驗(yàn)室在太陽(yáng)光譜波段模擬信號(hào)程序5S（Simulation of the Satellite Signal in the Solar Spectrum）的基礎(chǔ)上發(fā)展而來(lái)的，該模型考慮了地表的非朗伯特性，并且考慮了一些新的氣體影響，在計(jì)算透過(guò)率時(shí)加入了3 種新氣體（CH4、N2O、CO），采用了最新近似（state of the art）和逐次散射SOS（Successive Orders of Scattering）算法來(lái)計(jì)算散射和吸收，提高了瑞利散射與氣溶膠散射的計(jì)算精度［10-11］。該模型受目標(biāo)物類(lèi)型以及研究區(qū)特點(diǎn)的影響較小，具有較高的精度，且參數(shù)設(shè)置方面考慮了多種實(shí)際情況，適用于對(duì)400～4 000 nm 波長(zhǎng)范圍內(nèi)的遙感影像進(jìn)行大氣校正。缺點(diǎn)是不能處理球形大氣和limb（臨邊）觀測(cè)。

1.3 快速大氣校正模型法

快速大氣校正是從輻射傳輸模型出發(fā)，引入先驗(yàn)知識(shí)簡(jiǎn)化輻射傳輸模型并求解相關(guān)參數(shù)。該方法基于遙感圖像自身進(jìn)行校正，只借助中心波段即可從影像中獲取大氣補(bǔ)償參數(shù)從而快速實(shí)現(xiàn)大氣校正，不需要考慮其他光譜參數(shù)（太陽(yáng)天頂角、衛(wèi)星天頂角等）以及氣象信息（大氣氣溶膠厚度、風(fēng)速等），自動(dòng)從圖像上收集不同物質(zhì)的波譜信息，獲取經(jīng)驗(yàn)值完成高光譜和多光譜圖像的快速大氣校正。一般來(lái)說(shuō)，得到結(jié)果的精度近似FLAASH 或者其他基于輻射傳輸模型的＋/-15%［12］。目前支持的多光譜和高光譜波譜范圍是400～2 500 nm。

1.4 經(jīng)驗(yàn)線(xiàn)性法

經(jīng)驗(yàn)線(xiàn)性法是一個(gè)比較簡(jiǎn)便的定標(biāo)算法，國(guó)內(nèi)外已多次成功地利用該模型進(jìn)行遙感定標(biāo)實(shí)驗(yàn)。首先假設(shè)地面目標(biāo)的反射率與遙感器探測(cè)的信號(hào)之間具有線(xiàn)性關(guān)系，通過(guò)獲取遙感影像上特定地物的灰度值及其成像時(shí)相應(yīng)的地面目標(biāo)反射光譜的測(cè)量值，建立兩者之間的線(xiàn)性回歸方程式（1），在此基礎(chǔ)上對(duì)整幅遙感影像進(jìn)行輻射校正［13-14］。

式中：DNk（i）—第k個(gè)波段、第i個(gè)像元的灰度值；Refk（i）—對(duì)應(yīng)第k個(gè)波段、第i個(gè)像元反射率值；Ak—乘積因子，與大氣透過(guò)率和儀器增益有關(guān)；Bk—加減項(xiàng)，與暗電流和大氣程輻射有關(guān)。ELC 法的原理如圖1 所示。該方法數(shù)學(xué)和物理意義明確，計(jì)算簡(jiǎn)單，但必須以大量野外光譜測(cè)量為前提，因此成本較高，對(duì)野外工作依賴(lài)性強(qiáng)，且對(duì)地面定標(biāo)點(diǎn)的要求比較嚴(yán)格。

圖1 經(jīng)驗(yàn)線(xiàn)性模型原理示意圖Fig.1 Principal diagram of empirical linear model

鑒于航空高光譜數(shù)據(jù)航帶多、處理流程繁瑣的特點(diǎn)，針對(duì)上述4 種大氣校正模型方法開(kāi)發(fā)了適用于多景高光譜影像的大氣校正批處理軟件。軟件基于IDL（Interactive Data Language）語(yǔ)言，采用ENVI+IDL 模式進(jìn)行開(kāi)發(fā)，便于調(diào)用ENVI 已有的函數(shù)，提高了開(kāi)發(fā)的效率。上述4 種大氣校正模型方法均集成于軟件中，以供選擇使用，軟件界面如圖2 所示，為本研究試驗(yàn)區(qū)的航空高光譜數(shù)據(jù)不同大氣校正方法評(píng)價(jià)提供了高效便捷的軟件工具。

圖2 航空高光譜影像大氣校正批處理軟件界面Fig.2 Batch processing software interface for airborne hyperspectral image atmospheric correction

2 試驗(yàn)區(qū)與數(shù)據(jù)源

2.1 試驗(yàn)區(qū)概況

本巴圖試驗(yàn)區(qū)是具有鈾資源勘查潛力的遠(yuǎn)景區(qū)，位于巴丹吉林盆地東部，處于塔里木、哈薩克斯坦、西伯利亞、華北等4 大板塊的結(jié)合部位，區(qū)域構(gòu)造背景復(fù)雜多變。區(qū)內(nèi)經(jīng)歷了晚加里東、海西、印支等多期次的構(gòu)造巖漿活化，巖漿巖發(fā)育，活動(dòng)期次多，其中以華力西中晚期中酸性巖為主，印支期的花崗巖和加里東晚期中酸性巖次之，燕山期花崗巖呈零星分布，中酸性火山巖分布廣泛。區(qū)內(nèi)花崗巖分布面積約占基巖出露面積的三分之二，蝕變發(fā)育主要有絹云母化、綠泥石化、高嶺石化、碳酸鹽化等。

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

本巴圖地區(qū)的航空高光譜數(shù)據(jù)包含了本巴圖坳陷和圖克木隆起區(qū)，地理范圍大致為105°32'49"～106°0'57"E；40°41'21"～41°4'30"N，面積約1 702 km2；前人在區(qū)內(nèi)已發(fā)現(xiàn)多個(gè)鈾礦異常點(diǎn)，鈾成礦條件良好，具有較好的鈾找礦前景。

航空高光譜數(shù)據(jù)采集是利用引進(jìn)的加拿大ITRES 公司的SASI 成像光譜儀，數(shù)據(jù)采集時(shí)間為2021 年9—10 月。所采集的航空高光譜數(shù)據(jù)參數(shù)如表1 所示。

表1 本巴圖地區(qū)航空高光譜數(shù)據(jù)技術(shù)參數(shù)Table 1 Technical indexes of airborne hyperspectral data in Bembatu

在航空高光譜數(shù)據(jù)采集過(guò)程中，利用FieldSpec ASD 地物光譜儀同步獲取了明、暗地物反射光譜曲線(xiàn)，用于經(jīng)驗(yàn)線(xiàn)性法大氣校正。

3 試驗(yàn)結(jié)果

大氣校正之前，對(duì)試驗(yàn)區(qū)航空高光譜數(shù)據(jù)進(jìn)行了光譜和輻射校正、幾何校正等預(yù)處理。大氣校正采用上文中所述的4 種模型（FLAASH、6S、QUAC、ELC）對(duì)本巴圖地區(qū) 的SASI 航空高光譜影像進(jìn)行處理和對(duì)比評(píng)價(jià)。

3.1 大氣校正結(jié)果對(duì)比

基于試驗(yàn)區(qū)高光譜遙感礦物填圖應(yīng)用目的，重點(diǎn)評(píng)價(jià)大氣校正后獲得的各類(lèi)礦物的反射光譜是否準(zhǔn)確。由于本次SASI數(shù)據(jù)具有1.9 m的空間分辨率，混合像元問(wèn)題雖然無(wú)法避免，但與衛(wèi)星高光譜數(shù)據(jù)通常30 m 的空間分辨率相比，其像元混合效應(yīng)已大幅降低。因此，本次采用的評(píng)判依據(jù)參考USGS 標(biāo)準(zhǔn)光譜庫(kù)中的礦物標(biāo)準(zhǔn)光譜，重點(diǎn)對(duì)比不同方法校正后的光譜整體形態(tài)是否有差異，診斷性光譜特征位置是否準(zhǔn)確等。

圖3-圖7 分別為大氣校正后的試驗(yàn)區(qū)5 種典型蝕變礦物在航空高光譜影像中的反射光譜曲線(xiàn)與USGS 標(biāo)準(zhǔn)礦物光譜庫(kù)中的光譜曲線(xiàn)的對(duì)比情況，由于這些礦物的診斷性光譜特征均集中在2 000 nm 之后，故圖中展示的是2 000～2 450 nm 的波長(zhǎng)范圍。為了更好地進(jìn)行對(duì)照，在評(píng)價(jià)之前將USGS 標(biāo)準(zhǔn)礦物光譜重采樣為與SASI 航空高光譜數(shù)據(jù)一致的光譜分辨率。

由圖3a 可見(jiàn)，標(biāo)準(zhǔn)的高嶺石礦物的診斷性光譜特征位于2 165 nm 和2 210 nm，且兩處形成雙峰特征，在2 165 nm 處表現(xiàn)為一特殊的吸收肩狀特征。如圖3b 所示，4 種大氣校正模型校正后的結(jié)果均能夠反映高嶺石的上述特征，整體光譜形態(tài)趨勢(shì)亦基本一致，在2 400 nm 后的差異為公知的大氣中的水汽噪音所致。僅在反射率的高低存在一定差異，其中，F(xiàn)LAASH模型和6S 模型校正后的結(jié)果在反射率的高低上基本一致，QUAC 模型校正后的結(jié)果的反射率整體稍低，ELC 模型校正后的結(jié)果的反射率最低。

圖3 高嶺石光譜曲線(xiàn)對(duì)比圖Fig.3 Comparison of kaolinite spectral curves

如圖4a 所示，標(biāo)準(zhǔn)的碳酸鹽礦物方解石的診斷性光譜特征位于2 330～2 345 nm 附近，表現(xiàn)為非常明顯的單峰特征。如圖4b 所示，4 種大氣校正模型校正后的結(jié)果亦在2 330～2 345 nm具有明顯的單峰特征，且整體光譜形態(tài)趨勢(shì)亦基本一致，僅6S 模型在2 300 nm 處出現(xiàn)一處噪音。反射率的高低差異與高嶺石表現(xiàn)較一致，反射率的高低表現(xiàn)為FLAASH 模型≈6S 模型＞QUAC 模型＞ELC 模型。

圖4 方解石光譜曲線(xiàn)對(duì)比圖Fig.4 Comparison of calcite spectral curves

如圖5a 所示，標(biāo)準(zhǔn)的綠泥石礦物的診斷性光譜特征與碳酸鹽具有一定的重合，差異表現(xiàn)為綠泥石礦物在2 255 nm 附近具有次級(jí)吸收特征。如圖5b 所示，4 種大氣校正模型校正后的結(jié)果在綠泥石的2 255 nm 次級(jí)吸收特征均能夠獲得一定的反映。校正后整體光譜形態(tài)趨勢(shì)與標(biāo)準(zhǔn)光譜基本一致，僅FLAASH 模型與其他3 種模型結(jié)果存在少許形態(tài)差異。反射率的高低差異與前面兩種礦物一致，反射率的高低表現(xiàn)為FLAASH 模型≈6S 模型＞QUAC 模型＞ELC 模型。

圖5 綠泥石光譜曲線(xiàn)對(duì)比圖Fig.5 Comparison of chlorite spectral curves

如圖6a 所示，標(biāo)準(zhǔn)的中鋁絹云母礦物的診斷性光譜特征位于2 210 nm 和2 345 nm 附近，前者為Al-OH 特征，為高鋁、中鋁、低鋁絹云母區(qū)分的主要特征位置［15］。如圖6b 所示，4 種大氣校正模型校正后的結(jié)果的Al-OH 特征波長(zhǎng)位置均位于2 210 nm，次級(jí)特征位置與標(biāo)準(zhǔn)光譜基本一致，整體光譜形態(tài)6S 模型和ELC 模型更具有相似度，F(xiàn)LAASH 模型和QUAC 模型稍遜，但趨勢(shì)基本一致。反射率的高低差異與其他礦物一致，反射率的高低表現(xiàn)為FLAASH 模型≈6S 模型＞QUAC 模型＞ELC 模型。

圖6 中鋁絹云母光譜曲線(xiàn)對(duì)比圖Fig.6 Comparison of medium-Al sericite spectral curves

如圖7a 所示，標(biāo)準(zhǔn)的低鋁絹云母礦物的診斷性光譜特征位于2 225 nm 和2 345 nm 附近，前者為主要的Al-OH 特征。如圖7b 所示，4 種大氣校正模型校正后的結(jié)果的Al-OH 特征波長(zhǎng)位置均位于2 225 nm，次級(jí)特征位置均向長(zhǎng)波方向有所偏移，但4 種方法的偏移度基本一致。整體光譜形態(tài)6S 模型和ELC 模型更具有相似度，F(xiàn)LAASH 模型和QUAC 模型稍遜，但趨勢(shì)基本一致。反射率的高低差異與其他礦物一致，反射率的高低表現(xiàn)為FLAASH 模型≈6S模型＞QUAC 模型＞ELC 模型。

圖7 低鋁絹云母光譜曲線(xiàn)對(duì)比圖Fig.7 Comparison of low-Al sericite spectral curves

綜上所述，從試驗(yàn)區(qū)主要的5 種典型蝕變礦物光譜對(duì)比結(jié)果來(lái)看，在診斷性光譜特征方面，4 種大氣校正模型的校正結(jié)果均基本準(zhǔn)確。在光譜形態(tài)方面，ELC 相對(duì)更佳，F(xiàn)LAASH、6S、QUAC在部分礦物上表現(xiàn)較好，但在部分礦物上表現(xiàn)稍遜。在反射率高低差異方面，所有礦物在不同方法中的趨勢(shì)是一致的，反映了各模型較好的穩(wěn)定性。

3.2 填圖結(jié)果

選擇對(duì)礦物光譜形態(tài)保持效果較好的經(jīng)驗(yàn)線(xiàn)性法，對(duì)本巴圖試驗(yàn)區(qū)所有航帶影像進(jìn)行大氣校正和反射光譜重建。本巴圖試驗(yàn)區(qū)總計(jì)有4 個(gè)時(shí)相的航空高光譜數(shù)據(jù)，不同時(shí)相的航空高光譜影像分別采用該時(shí)相的地面定標(biāo)場(chǎng)測(cè)量數(shù)據(jù)進(jìn)行經(jīng)驗(yàn)線(xiàn)性法大氣校正?；诤娇崭吖庾V數(shù)據(jù)大氣校正批處理軟件，分4 次分別在軟件界面中設(shè)置好輸入和輸出文件路徑，填寫(xiě)有關(guān)參數(shù)，運(yùn)行軟件即可調(diào)用大氣校正批處理模式，之后的處理均在無(wú)人值守下自動(dòng)進(jìn)行。

應(yīng)用經(jīng)大氣校正后的航空高光譜反射率影像數(shù)據(jù)，基于光譜沙漏（Spectral Hourglass）方法提取礦物端元，結(jié)合光譜相似性匹配方法，完成了本巴圖試驗(yàn)區(qū)航空高光譜礦物填圖（圖8），為分析本巴圖鈾資源勘查遠(yuǎn)景區(qū)主要目標(biāo)層的次生蝕變、成礦構(gòu)造等提供了重要信息。

圖8 本巴圖試驗(yàn)區(qū)航空高光譜遙感礦物填圖Fig.8 Airborne hyperspectral remote sensing mineral mapping in Benbatu test area

4 結(jié)論

以巴丹吉林盆地東部本巴圖鈾資源勘查遠(yuǎn)景區(qū)為試驗(yàn)區(qū)，面向地質(zhì)礦產(chǎn)勘查中的航空高光譜遙感礦物填圖應(yīng)用目的，基于FLAASH、6S、快速大氣校正和經(jīng)驗(yàn)線(xiàn)性法4 種經(jīng)典大氣校正模型開(kāi)發(fā)了航空高光譜數(shù)據(jù)大氣校正批處理軟件，對(duì)試驗(yàn)區(qū)SASI 航空高光譜影像進(jìn)行了大氣校正和礦物反射光譜重建評(píng)價(jià)試驗(yàn)。結(jié)果表明：4 種大氣校正模型方法均能較準(zhǔn)確地反演5 種典型蝕變礦物的診斷性光譜特征，其中經(jīng)驗(yàn)線(xiàn)性法對(duì)礦物的光譜形態(tài)校正效果相對(duì)更佳。最后，在大氣校正的基礎(chǔ)上完成了本巴圖鈾資源勘查遠(yuǎn)景區(qū)礦物填圖，為分析和評(píng)價(jià)該地區(qū)鈾成礦潛力提供了技術(shù)支撐。