航空高光譜地表氧化鐵礦物含量評估技術

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

(1. 核工業北京地質研究院 遙感信息與圖像分析技術國家級重點實驗室，北京 100029；2.中國地質大學(北京) 國土資源與高新技術研究中心，北京 100083)

0 引言

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

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

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

1 方法

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

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

1.1 數據預處理

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

由于航空數據采集的特殊性，每次成像光譜儀的安裝和飛行過程的差異，會造成實驗室定標參數的誤差，成像光譜儀本身的smile效應等因素也可能會造成反射率數據的錯誤。可以根據氣體吸收特征等方法迭代進行定標校正，提高反射率的精度。

1.2 連續統去除

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

1.3 吸收峰位置和深度計算

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

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

(1)

公式中擬合系數[a，b，c]可以將連續統去除后的3個波段數據通過直接帶入法計算。

(2)

(3)

(4)

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

(5)

(6)

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

(7)

1.4 回歸分析

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

2 實驗與分析

2.1 研究區概況

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

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

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

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

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

2.2 航空高光譜數據獲取和預處理

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

表1 航空高光譜基本技術參數

2.3 地面光譜數據獲取與分析測試

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

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

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

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

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

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

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

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

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

圖4 樣品光譜薄片鑒定結果

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

2.4 吸收特征計算與回歸分析

選擇850～1 000 nm波譜對航空高光譜數據進行連續統去除，對航空和地面數據采用二次擬合方法，計算航空光譜數據的吸收深度位置和相對應的吸收深度，利用線性回歸方法，對采樣點的光譜吸收深度和Fe2O3含量進行反演。圖5顯示了不同采樣區的反演結果。

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

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

3 實驗結果

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

Y=59.42+122.94X1-237.49X2

(8)

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

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

表3 氧化鐵含量回歸分析預測結果及誤差

圖6 黑溝鐵礦區赤鐵礦豐度圖

4 結束語

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

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