秦凱，陳建平，趙英俊，朱玲，楊國防

(1. 核工業(yè)北京地質(zhì)研究院 遙感信息與圖像分析技術(shù)國家級重點(diǎn)實(shí)驗(yàn)室，北京 100029；2.中國地質(zhì)大學(xué)(北京) 國土資源與高新技術(shù)研究中心，北京 100083)

0 引言

赤鐵礦、針鐵礦和褐鐵礦等含氧化鐵礦物分布廣泛，在可見-近紅外(VNIR)波段產(chǎn)生診斷性的吸收光譜特征[1-2]。這些吸收特性來源于電子過程，主要包括兩種類型：第一種是兩個(gè)不同離子間的價(jià)間電荷轉(zhuǎn)移(intervalence charge transfer，IVCT)，例如三價(jià)鐵和氧離子的電荷轉(zhuǎn)移。在光譜上，主要表現(xiàn)為480～550 nm的強(qiáng)吸收。第二種是離子的晶體場或能級的電荷轉(zhuǎn)移。對于蝕變礦物巖石光譜，三價(jià)鐵離子晶體場轉(zhuǎn)移(crystal field transitions，CFT)會在630～715 nm和850～1 000 nm之間產(chǎn)生吸收特征。亨特和阿什莉觀察到，鐵氧化物的含量與在850～1 000 nm之間的晶體場吸收光譜的深度特征有關(guān)[1]，這也被Cudahy和Ramanaidou的實(shí)驗(yàn)工作進(jìn)行了驗(yàn)證，得出結(jié)論：Fe2O3的含量與850～1 000 nm的吸收深度相關(guān)[2-3]。

可是，采用航空高光譜技術(shù)獲取的巖礦光譜還與傳感器采樣間隔、噪聲和地形等因素有關(guān)。這些因素都可能影響吸收深度，降低氧化鐵定量評估的精度。在之前的研究中，Rodger等設(shè)計(jì)了一種簡單二次多項(xiàng)式擬合方法，可以對不同光譜分辨率的高光譜數(shù)據(jù)進(jìn)行光譜吸收峰位置和深度信息的提取[4]。Murphy等指出地形會影響氧化鐵礦物的吸收深度[5-6]。連續(xù)統(tǒng)去除(continuum removal)方法可以有效去除地形對光譜吸收特征的影響[7-9]。除此之外，不同巖性的理化特性可能會對氧化鐵在光譜響應(yīng)上具有不同的效果，所以對于氧化鐵的光譜反演，需要不同巖性、不同礦物的混合光譜的先驗(yàn)知識。張琴等以美國地質(zhì)調(diào)查局(USGS)、噴氣推進(jìn)實(shí)驗(yàn)室(JPL)和約翰霍普金斯大學(xué)(JHU)提供的標(biāo)準(zhǔn)數(shù)字光譜庫為基礎(chǔ)數(shù)據(jù)，基于去連續(xù)統(tǒng)后小波包分解高頻和低頻部分運(yùn)用多元逐步回歸法、偏最小二乘法建模建立了Fe2O3含量的13個(gè)高光譜反演模型[10]。閆柏琨等以亞利桑那州國家大學(xué)(Arizona State University，ASU)光譜庫為數(shù)據(jù)源，研究了硅酸鹽礦物發(fā)射率光譜克里斯琴森特征與SiO2含量的定量關(guān)系[11]。

前人的研究大多基于實(shí)驗(yàn)室光譜或者地面光譜，針對航空高光譜數(shù)據(jù)獲取、處理和分析的特點(diǎn)，基于趙英俊等在祁連山黑溝鐵礦等典型礦區(qū)的航空高光譜礦物填圖和礦產(chǎn)預(yù)測研究基礎(chǔ)[12]，結(jié)合連續(xù)統(tǒng)去除、二次擬合、多元回歸分析等方法，構(gòu)建了基于航空高光譜數(shù)據(jù)的地表氧化鐵含量評估方法，通過鐵離子晶體場轉(zhuǎn)移光譜吸收特征和氧化鐵含量的回歸分析，獲得了研究區(qū)氧化鐵礦物含量分布圖。

1 方法

通過借鑒前人方法和結(jié)合航空高光譜特點(diǎn)，提出了一套地表氧化鐵含量的反演方法，該方法的流程：首先基于大氣特征譜線的定標(biāo)參數(shù)進(jìn)行數(shù)據(jù)預(yù)處理；其次利用連續(xù)統(tǒng)去除和二次擬合的方法獲取鐵離子CFT光譜特征；最后根據(jù)地面實(shí)測氧化鐵含量進(jìn)行回歸分析。流程圖如圖1所示。

圖1 航空高光譜數(shù)據(jù)的地表氧化鐵含量評估流程圖

1.1 數(shù)據(jù)預(yù)處理

航空高光譜預(yù)處理就是對成像光譜儀獲取的原始輻射數(shù)據(jù)進(jìn)行定標(biāo)和大氣輻射校正。預(yù)處理的目的主要是消除大氣和光照等因素對地物反射的影響，獲得地表反射率數(shù)據(jù)。航空高光譜數(shù)據(jù)的反射率計(jì)算主要基于大氣輻射傳輸模型來計(jì)算，其中大氣成分的參數(shù)主要通過地面實(shí)測和吸收線評估獲取，水汽成分的評估可以采用大氣預(yù)處理微分吸收法逐像元的評估水汽含量，二氧化碳可以采用地面實(shí)測數(shù)據(jù)，氣溶膠采用鄉(xiāng)村模型。各參數(shù)確定之后，就可以計(jì)算出航空高光譜圖像數(shù)據(jù)的各像元的反射率光譜。

由于航空數(shù)據(jù)采集的特殊性，每次成像光譜儀的安裝和飛行過程的差異，會造成實(shí)驗(yàn)室定標(biāo)參數(shù)的誤差，成像光譜儀本身的smile效應(yīng)等因素也可能會造成反射率數(shù)據(jù)的錯(cuò)誤?？梢愿鶕?jù)氣體吸收特征等方法迭代進(jìn)行定標(biāo)校正，提高反射率的精度。

1.2 連續(xù)統(tǒng)去除

在計(jì)算光譜特征并反演物理量之前，需進(jìn)行光譜的連續(xù)統(tǒng)去除。利用直線連接光譜起始與終止波長處反射率，搜尋整個(gè)譜段上反射率大于連接直線且差值最大的波段，用該波段的反射率代替連接直線在該波段的值，生成連接起始波長、終止波長與該波段的折線，上述過程反復(fù)迭代循環(huán)，直至在所有波段上，反射率均小于或等于生成的折線，最終形成的折線即為連續(xù)統(tǒng)。在可見至短波紅外波譜范圍內(nèi)，礦物的吸收特征可以看作為礦物吸收和背景吸收2個(gè)部分。其中背景吸收(連續(xù)統(tǒng))是疊加在礦物吸收特征之上，會造成吸收峰位置的偏移。因此，采用去連續(xù)的方法來減少背景信息和混合礦物對光譜的影響。連續(xù)統(tǒng)去除的方法可以減少地形因素的影響，巖石的反射光譜和連續(xù)統(tǒng)光譜都包含相同的地形影響因子，在比值算法計(jì)算下，地形影響被消除。

1.3 吸收峰位置和深度計(jì)算

本文采用二次擬合的方法計(jì)算鐵離子吸收峰位置和吸收深度，該方法不受高光譜設(shè)備的光譜采樣間隔、信噪比等因素影響，波長和連續(xù)統(tǒng)去除后的反射率可以用公式(1)來表示：

ρcr(λ)=aλ2+bλ+c

(1)

公式中擬合系數(shù)[a，b，c]可以將連續(xù)統(tǒng)去除后的3個(gè)波段數(shù)據(jù)通過直接帶入法計(jì)算。

(2)

(3)

(4)

對公式(1)計(jì)算一階導(dǎo)數(shù)，計(jì)算導(dǎo)數(shù)為零值的波長，將直接帶入法計(jì)算出的擬合系數(shù)帶入導(dǎo)數(shù)公式。

(5)

(6)

λab為吸收峰位置的波長，進(jìn)一步求解吸收深度：

(7)

1.4 回歸分析

根據(jù)地面采樣分析測試的Fe2O3含量與對應(yīng)光譜的吸收特征，基于線性回歸的方法建立特定吸收波段吸收深度與氧化鐵含量的回歸方程。應(yīng)用回歸方程對航空成像光譜數(shù)據(jù)進(jìn)行氧化鐵含量的反演。

2 實(shí)驗(yàn)與分析

2.1 研究區(qū)概況

黑溝鐵礦位于甘肅省嘉峪關(guān)市鏡鐵山東南14 km處，有簡易公路可以到達(dá)，地理位置：97°57′31″E，39°18′24″N。

黑溝鐵礦區(qū)出露地層為長城系樺樹溝組和薊縣系鏡鐵山群，由淺變質(zhì)的雜色千枚巖、石英巖及鐵礦層組成，總體上為一套復(fù)式含鐵建造，呈北西—南東向展布，為西部樺樹溝礦區(qū)含鐵巖系東延部分。礦區(qū)內(nèi)巖漿巖僅見石英閃長玢巖及輝綠巖。石英閃長玢巖見于礦區(qū)南部，侵入灰綠色千枚巖及鐵礦層中，與巖層、礦層走向基本平行，局部斜切，呈脈狀產(chǎn)出。輝綠巖為巖墻產(chǎn)出，見于礦區(qū)西部，侵入灰綠色千枚巖中(圖2)。

注：此圖據(jù)甘肅省地質(zhì)局鏡鐵山地質(zhì)資料修編。圖2 黑溝礦床地質(zhì)圖

黑溝礦區(qū)本身為一向斜構(gòu)造，位于樺樹溝復(fù)向斜東南端，軸向北西西，總體上向南東傾沒。礦區(qū)范圍內(nèi)，黑溝向斜兩端則向中間傾伏，向斜兩翼次一級小型褶皺發(fā)育，次一級褶皺軸向與主向斜軸向基本一致，屬于同期構(gòu)造。

礦區(qū)內(nèi)斷裂構(gòu)造有北東東和北北東兩組，斷層面多傾向北北西或北西西。破碎帶不甚發(fā)育且延長不遠(yuǎn)，斷層整體上南東盤相對升起，在垂直陷落的同時(shí)，斷層北西盤普遍向南西方向位移，對于礦體破壞不大。礦區(qū)內(nèi)巖層及礦層中，節(jié)理發(fā)育，其中以鐵礦層及灰色千枚巖中的節(jié)理尤為發(fā)育。

2.2 航空高光譜數(shù)據(jù)獲取和預(yù)處理

采用空中國王飛機(jī)搭載CASI/SASI航空高光譜傳感器，在2014年7月完成了祁連山黑溝鐵礦研究區(qū)的航空高光譜數(shù)據(jù)。飛行前完成了成像光譜儀的實(shí)驗(yàn)室定標(biāo)，所獲取的航空高光譜數(shù)據(jù)的參數(shù)見表1。完成了工作區(qū)CASI/SASI數(shù)據(jù)的輻射校正、正射校正、幾何校正、色差處理、影像拼接和大氣校正等數(shù)據(jù)預(yù)處理工作。主要包括6個(gè)步驟：輻射校正；傳感器姿態(tài)數(shù)據(jù)處理；GPS定位數(shù)據(jù)處理；姿態(tài)數(shù)據(jù)與定位數(shù)據(jù)時(shí)間同步與集成；精細(xì)的幾何校正；地形校正；大氣校正和光譜重建。其中，采用的是5 m分辨率衛(wèi)星DEM數(shù)據(jù)完成了正射校正。大氣成分的參數(shù)主要通過地面實(shí)測和吸收線評估獲取，其中水汽成分的評估采用大氣預(yù)處理微分吸收法逐像元的評估水汽含量，二氧化碳采用地面實(shí)測數(shù)據(jù)357.5 ppm(CO2混合比)，氣溶膠采用鄉(xiāng)村模型。大氣參數(shù)確定之后計(jì)算出研究區(qū)CASI/SASI的各像元反射率光譜。

表1 航空高光譜基本技術(shù)參數(shù)

2.3 地面光譜數(shù)據(jù)獲取與分析測試

采用可見光-短波紅外ASD FieldSpec Pro FR光譜儀采集了研究區(qū)地面光譜數(shù)據(jù)。典型巖石的光譜曲線見圖3。

在黑溝鐵礦開展了多條地質(zhì)剖面的調(diào)查和光譜測量，獲取了地面典型地質(zhì)單元的光譜數(shù)據(jù)，用于礦物高光譜礦物信息提取和分析。以黑溝礦區(qū)的“3790開采平臺”高光譜地面測量剖面為例，主要測量了千枚巖絹云母粉砂質(zhì)板巖(QLYZ-7A QLYZ-12A)、赤鐵礦化砂巖(QLYZ-7B QLYZ-8B)、綠泥粉砂質(zhì)板巖(QLYZ-8A)、鐵白云巖(QLYZ-9B、QLYZ-10B、QLYZ-11A)和鐵礦鐵礦石(QLYZ-11B)等巖石的光譜。有以下特征：

(1)千枚巖絹云母粉砂質(zhì)板巖的光譜曲線顯示具有3個(gè)特征吸收位置，分別位于2 210 nm、2 345 nm和2 440 nm，其主要吸收峰位于2 210 nm，與絹云母礦物的鋁羥基基團(tuán)振動(dòng)有關(guān)。

圖3 典型蝕變巖石光譜吸收特征

(2)赤鐵礦化砂巖由于其赤鐵礦的三價(jià)鐵離子的電子作用，光譜顯示特征吸收峰位于880 nm，同時(shí)在其左肩存在770 nm的強(qiáng)反射峰。

(3)綠泥粉砂質(zhì)板巖的光譜顯示了鐵羥基2 255 nm和2 360 nm的吸收特征組合。

(4)黑溝礦區(qū)碳酸鹽主要為鐵白云石，碳酸根特征吸收位置分別位于2 330 nm和二價(jià)鐵的吸收谷從1 000～1 300 nm。

(5)鐵礦石主要礦物為鏡鐵礦及少量磁鐵礦，無明顯吸收特征，該樣品表面存在少量褐鐵礦化現(xiàn)象。

從以上光譜分析可以得出，礦區(qū)內(nèi)的巖性光譜差異較大，赤鐵礦、絹云母、綠泥石、白云石等礦物光譜特征明顯，以上光譜特征位置、深度和組合關(guān)系可以作為航空高光譜信息提取的依據(jù)和參考。對測量點(diǎn)進(jìn)行了地面采樣，并開展了地面樣品的薄片鑒定(圖4)和分析測試，部分測試結(jié)果見表2。

圖4 樣品光譜薄片鑒定結(jié)果

表2 典型礦床剖面樣品主要化學(xué)成分含量%

2.4 吸收特征計(jì)算與回歸分析

選擇850～1 000 nm波譜對航空高光譜數(shù)據(jù)進(jìn)行連續(xù)統(tǒng)去除，對航空和地面數(shù)據(jù)采用二次擬合方法，計(jì)算航空光譜數(shù)據(jù)的吸收深度位置和相對應(yīng)的吸收深度，利用線性回歸方法，對采樣點(diǎn)的光譜吸收深度和Fe2O3含量進(jìn)行反演。圖5顯示了不同采樣區(qū)的反演結(jié)果。

鐵礦層的主要巖石含礦砂巖、含礦硅質(zhì)巖，圖5(a)為12件樣品的光譜與氧化鐵含量的相關(guān)關(guān)系，鐵礦石光譜鐵離子晶體場吸收深度特征與Fe2O3含量具有明顯的線性關(guān)系，相關(guān)系數(shù)R2=0.871 4；在鐵礦層周邊的褐鐵礦化、赤鐵礦化白云巖、赤鐵礦化千枚巖、赤鐵礦化砂巖采樣18件，圖5(b)顯示的光譜鐵離子晶體場吸收深度特征與Fe2O3含量具有明顯的線性關(guān)系，相關(guān)系數(shù)R2=0.833 4。礦床區(qū)域千枚巖氧化鐵含量較低，線性相關(guān)度低，相關(guān)系數(shù)R2=0.181 8(圖5(c))。

圖5 不同采樣區(qū)吸收深度 (吸收峰位置900 nm)與氧化鐵含量相關(guān)關(guān)系

3 實(shí)驗(yàn)結(jié)果

對黑溝鐵礦的航空高光譜數(shù)據(jù)進(jìn)行了數(shù)據(jù)預(yù)處理、連續(xù)統(tǒng)去除、二次擬合等處理，根據(jù)地面和航空高光譜數(shù)據(jù)二次擬合計(jì)算的吸收峰位置和深度結(jié)果，獲取了850～1 000 nm的吸收峰位置和深度的鐵離子晶體場轉(zhuǎn)移光譜特征。通過對不同采樣區(qū)域的回歸結(jié)果分析，引入其他波段的光譜變量，開展了多元回歸迭代分析，最終，建立了研究區(qū)的航空高光譜地表氧化鐵含量的評估公式：

Y=59.42+122.94X1-237.49X2

(8)

式中：Y為氧化鐵含量預(yù)測值;X1為吸收峰位置900 nm的光譜吸收深度;X2為2 100～2 280 nm的平均光譜。回歸系數(shù)R2=0.965 2，標(biāo)準(zhǔn)誤差為4.44wt%。表3為部分采樣點(diǎn)回歸分析的預(yù)測值和殘差。

根據(jù)回歸分析的預(yù)測結(jié)果，計(jì)算了研究區(qū)地表氧化鐵含量豐度(圖6)，可以更加清晰地展示出黑溝鐵礦該處礦體的分布規(guī)律，氧化鐵含量的區(qū)域特征與礦體分布形態(tài)基本一致，豐度較高的礦物信息表現(xiàn)為長軸北西向的封閉環(huán)形。

表3 氧化鐵含量回歸分析預(yù)測結(jié)果及誤差

圖6 黑溝鐵礦區(qū)赤鐵礦豐度圖

4 結(jié)束語

本文面向航空高光譜遙感影像，充分利用氧化鐵礦物的吸收光譜特征，提出了一種地表氧化鐵含量的評估方法。實(shí)驗(yàn)結(jié)果表明該方法能有效地反演出地表氧化鐵含量。本文方法在祁連山黑溝鐵礦進(jìn)行了應(yīng)用，對礦區(qū)內(nèi)的鐵含量評估較為準(zhǔn)確，該方法仍然需要地面的高光譜數(shù)據(jù)和地質(zhì)背景知識，如何實(shí)現(xiàn)在未知區(qū)域上快速分析和評估地表氧化鐵含量將是進(jìn)一步研究的方向。

致謝：特此感謝核工業(yè)北京地質(zhì)研究院劉德長研究員、中國國土資源航空物探遙感中心王潤生教授的對工作的指導(dǎo)和支持。