白 洋，宋唐雷，賈玉娜，李孟倩，康會(huì)濤

(1.華北理工大學(xué)礦業(yè)工程學(xué)院，河北 唐山 063210；2.河北省地質(zhì)礦產(chǎn)勘查開發(fā)局第七地質(zhì)大隊(duì)，河北 三河 065201)

尾礦庫(kù)是礦山開采過程中形成的最大危險(xiǎn)源之一，嚴(yán)重威脅著尾礦庫(kù)周圍的生態(tài)環(huán)境及周邊居民的生命財(cái)產(chǎn)安全[1-3]。因此，實(shí)時(shí)、快速、高效地掌握尾礦庫(kù)的分布情況，對(duì)提高礦產(chǎn)資源循環(huán)利用能力具有十分重要的現(xiàn)實(shí)意義。傳統(tǒng)的尾礦庫(kù)監(jiān)測(cè)方法多以實(shí)地調(diào)查、勘測(cè)為主，耗費(fèi)大量的人力、物力，且監(jiān)測(cè)結(jié)果存在一定的主觀性和局限性，無(wú)法滿足大面積尾礦監(jiān)測(cè)。隨著遙感技術(shù)在空間分辨率、光譜分辨率和時(shí)間分辨率方面的顯著提升[4-6]，利用遙感技術(shù)探測(cè)范圍廣、數(shù)據(jù)獲取快捷、受自然條件約束小等優(yōu)勢(shì)，可大大提高尾礦庫(kù)監(jiān)測(cè)的成效[7-8]。閻永忠等[9]利用遙感技術(shù)和GIS技術(shù)系統(tǒng)對(duì)尾礦庫(kù)進(jìn)行監(jiān)測(cè)分析；申彥科等[10]依據(jù)TM影像總結(jié)出尾礦庫(kù)增容會(huì)對(duì)周圍的生態(tài)環(huán)境產(chǎn)生的影響；方雪娟等[11]對(duì)陳貴鎮(zhèn)小型尾礦庫(kù)周圍的4類地物進(jìn)行了危害分級(jí)分析；南竣祥等[12]對(duì)宣威市煤礦開采情況展開遙感解譯，提升了尾礦庫(kù)遙感監(jiān)測(cè)技術(shù)的作用；張?jiān)朴⒌萚13]利用高分一號(hào)影像對(duì)尾礦庫(kù)進(jìn)行目視解譯最佳波段組合，為運(yùn)用遙感影像進(jìn)行尾礦庫(kù)研究提供了依據(jù)；凌子燕[14]利用GF-1遙感衛(wèi)星，采取支持向量機(jī)(support vector machine，SVM)建立遙感反演模型，提高了尾礦庫(kù)遙感監(jiān)測(cè)精度；曹蘭杰等[15]對(duì)典型地區(qū)鐵尾礦樣本進(jìn)行了光譜測(cè)量，并利用SAM(spectral angle mapper)方法和SFF(spectral feature fitting)方法進(jìn)行匹配，優(yōu)化了河北省兩種礦床類型鐵尾礦的識(shí)別窗口；郝利娜等[16]從尾砂的光譜、紋理特征等因素考慮，根據(jù)WorldView-2遙感影像對(duì)鄂東南尾礦庫(kù)規(guī)模進(jìn)行了尾礦識(shí)別因子研究；馬國(guó)超等[17]集合高分遙感衛(wèi)星數(shù)據(jù)和無(wú)人機(jī)遙感數(shù)據(jù)，利用其快速排查，精準(zhǔn)監(jiān)控的優(yōu)勢(shì)，極大提高了礦山安全監(jiān)測(cè)效率；汪金花等[18-19]對(duì)比分析了差異化粒徑、干濕狀態(tài)下的鐵尾砂樣本光譜特征，擬合確定了鐵尾礦高光譜識(shí)別的有效窗口，經(jīng)光譜特征參量及光譜匹配分析獲得了鐵尾礦特征識(shí)別波段，得到了鐵尾礦多源信息提取的可行性。多源遙感數(shù)據(jù)具有一定的優(yōu)勢(shì)，但同樣存在局限性，如多光譜數(shù)據(jù)空間信息豐富但無(wú)法區(qū)分“同物異譜”或“同譜異物”，而高光譜數(shù)據(jù)光譜信息豐富但數(shù)據(jù)量大，信息冗余，存在Hughes現(xiàn)象。因此，如何利用高光譜數(shù)據(jù)和多光譜數(shù)據(jù)的優(yōu)勢(shì)進(jìn)行尾礦庫(kù)信息提取尤為重要。

鑒于此，本文以唐山市司家營(yíng)尾礦庫(kù)為例，利用實(shí)測(cè)尾砂波譜信息分別與多光譜Landsat8 OLI影像數(shù)據(jù)和高光譜珠海一號(hào)影像數(shù)據(jù)的端元波譜進(jìn)行匹配，并采用光譜角分類的方法進(jìn)行尾礦信息填圖，獲取尾礦庫(kù)的位置信息。

1 研究區(qū)概況與數(shù)據(jù)源

本文研究區(qū)域是河北省唐山市灤州市司家營(yíng)尾礦庫(kù)，位于灤州市城南10 km處，屬于高硅鞍山型鐵尾礦。礦體南北長(zhǎng)10 km，東西2 km，礦石資源總儲(chǔ)量為23.48億t，屬于我國(guó)整裝特大型鐵礦床之一。司家營(yíng)尾礦庫(kù)距選礦廠約12 km，尾礦庫(kù)匯水面積為2.5 km2，平均坡度4.1°。尾礦庫(kù)由東、西兩個(gè)溝組成，分別建有初期壩和副壩。尾礦堆積壩最終堆積標(biāo)高為128 m，總壩高93 m，總庫(kù)容約8 714萬(wàn)m3[20]。

2 數(shù)據(jù)源介紹與預(yù)處理

2.1 數(shù)據(jù)源介紹

Landsat8衛(wèi)星于2013年2月發(fā)射成功，是目前世界范圍內(nèi)應(yīng)用領(lǐng)域最廣的民用對(duì)地觀測(cè)衛(wèi)星，本次數(shù)據(jù)選用采集于2019年10月29日的一幅OLI影像，共9個(gè)波段，包括8個(gè)30 m的多光譜波段和1個(gè)15 m的全色波段。珠海一號(hào)高光譜衛(wèi)星于2018年4月發(fā)射成功，其采集的數(shù)據(jù)空間分辨率達(dá)到10 m，光譜分辨率達(dá)到2.5 nm，共256個(gè)波段，均在400～1 000 nm之間。本文使用的60個(gè)尾砂樣本均取自河北省灤州市司家營(yíng)尾礦庫(kù)，由ASD Field Spec Pro光譜儀對(duì)各樣本進(jìn)行測(cè)量，每個(gè)樣本均測(cè)量5次，每5 min測(cè)一次，白板進(jìn)行校正，最終得到60個(gè)樣本的平均波譜曲線。

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

根據(jù)數(shù)據(jù)源的具體情況，對(duì)Landsat8 OLI影像數(shù)據(jù)進(jìn)行了輻射定標(biāo)、大氣校正等預(yù)處理，提高了影像的定位精度、空間分辨率及光譜分辨率，有利于進(jìn)行尾礦庫(kù)的精準(zhǔn)識(shí)別。以研究區(qū)的Landsat8 OLI影像文件和珠海一號(hào)影像文件為基準(zhǔn)，對(duì)去除包絡(luò)線后的波譜曲線進(jìn)行重新采樣，確保實(shí)測(cè)波譜與影像波譜范圍相同，方便后續(xù)開展實(shí)測(cè)波譜與端元波譜的匹配。

3 光譜匹配實(shí)驗(yàn)

影像數(shù)據(jù)經(jīng)最小噪聲分離(minimum noise fractions，MNF)達(dá)到降維和去噪目的，通過MNF特征值曲線選擇信息量較大的幾個(gè)分量，進(jìn)行純凈像元指數(shù)(PPI)計(jì)算，通過PPI求得的端元信息，在N維可視化中找到緊密抱團(tuán)的端元，并繪制其特征波譜曲線，經(jīng)與實(shí)測(cè)的尾礦自定義波譜庫(kù)進(jìn)行匹配，并進(jìn)行尾礦庫(kù)填圖。

3.1 端元波譜提取

端元波譜是地物識(shí)別分類中十分重要的參考光譜，端元波譜的質(zhì)量直接影響地物識(shí)別結(jié)果的精度。預(yù)處理后的影像進(jìn)行MNF變換后可以得到各波段的特征信息，如圖1和圖2所示。從圖1和圖2中可以看出沿X軸方向各特征數(shù)包含的特征信息逐漸遞減。圖1為L(zhǎng)andsat8 OLI影像波段經(jīng)過MNF變換后的特征圖，其中前3個(gè)分量的特征值比較大，第4分量以后特征曲線的特征值幾乎都在5以下，因此，選取前3個(gè)特征波段進(jìn)行PPI計(jì)算。圖2為珠海一號(hào)影像的特征曲線，其中前4個(gè)波段涵蓋了大部分的圖像信息，從第5分量開始曲線基本保持平穩(wěn)，后面波段的影像信息幾乎全部為噪聲，所以后續(xù)PPI計(jì)算時(shí)只選取了前4個(gè)波段，以大幅度減少處理工作量和處理時(shí)間。

圖1 Landsat8 OLI MNF變換特征曲線Fig.1 MNF transformation characteristic curve of Landsat8 OLI

圖2 珠海一號(hào)MNF變換特征曲線Fig.2 MNF transformation characteristic curve of Zhuhai-1

純凈像元指數(shù)(pixel purity index，PPI)是在數(shù)據(jù)降維處理后進(jìn)行純凈像元提取的過程，PPI越大，說明與之對(duì)應(yīng)的像元越純凈，對(duì)應(yīng)的PPI結(jié)果圖就越亮。在影像進(jìn)行MNF變換之后，經(jīng)過多次實(shí)驗(yàn)對(duì)Landsat8 OLI影像進(jìn)行PPI計(jì)算采用了10 000次迭代，閾值系數(shù)設(shè)置為2.5，珠海一號(hào)影像迭代次數(shù)為10 000次，閾值系數(shù)設(shè)置為2.5，如圖3和圖4所示，其中，高亮區(qū)域?yàn)镻PI提取后的極值像元的分布位置，這些純凈像元多聚集在尾礦庫(kù)中央的位置，還有一些零散的分布在尾礦庫(kù)的外圍。

圖3 Landsat8 OLI影像PPI結(jié)果圖Fig.3 PPI results of Landsat8 OLI images

圖4 珠海一號(hào)影像PPI結(jié)果圖Fig.4 PPI results of Zhuhai-1 images

將MNF和PPI計(jì)算后的結(jié)果輸入到N維可視化中，選取不同的波段進(jìn)行旋轉(zhuǎn)，找到多個(gè)獨(dú)立的點(diǎn)集群，如圖5(a)和圖6(a)所示，圖中有明顯聚集的4個(gè)團(tuán)簇，這4個(gè)團(tuán)簇代表了4種“純凈像元”，以此定義為端元，圖5(b)和圖6(b)為得到的4種端元的波譜曲線，通過現(xiàn)場(chǎng)采集的尾砂樣本制作的自定義波譜庫(kù)可以識(shí)別出這4種像元代表的物質(zhì)和分布的區(qū)域。

圖5 Landsat8 OLI影像N維可視化圖Fig.5 N-dimensional visualization of Landsat8 OLI images

圖6 珠海一號(hào)影像N維可視化圖Fig.6 N-dimensional visualization of Zhuhai-1 images

3.2 光譜匹配算法

在端元提取出來(lái)之后，還不能明確各端元代表哪種地物類型，需要通過波譜分析對(duì)端元進(jìn)行識(shí)別，同時(shí)要修改波譜曲線分辨率的縮放系數(shù)，使端元波譜與波譜庫(kù)具有相同的分辨率，此次研究分別使用光譜角度匹配(spectral angle match，SAM)和波譜特征擬合(spectral feature fitting，SFF)2種方法進(jìn)行識(shí)別。

3.2.1 光譜角度匹配法

光譜角度匹配(spectral angle match，SAM)是將待驗(yàn)證光譜和參考光譜均作為高維向量，通過計(jì)算二者的夾角大小來(lái)確定匹配程度[21]，夾角越小說明相似程度越高，是同類地物的可能性越大。

計(jì)算夾角通常通過兩個(gè)向量的余弦來(lái)求得，公式見式(1)。

(1)

式中：X*為像元光譜的向量；Xi為參考光譜的向量。

3.2.2 波譜特征擬合法

波譜特征擬合(spectral feature fitting，SFF)是一種利用波譜曲線的吸收特征對(duì)像元波譜和參考波譜進(jìn)行匹配的方法，使用最小二乘法逐次迭代以評(píng)價(jià)兩者相似程度[22]。SFF會(huì)對(duì)輸入的像元波譜進(jìn)行包絡(luò)線去除、刪去背景信息、增強(qiáng)特征波段的吸收特性，比較吸收谷的位置和深度，將參考波譜與像元波譜相匹配。

4 匹配結(jié)果與分析

4.1 端元匹配結(jié)果

N維可視化后的端元波譜需要通過波譜分析與已知物質(zhì)的波譜進(jìn)行匹配，方可得知哪一個(gè)端元最能代表需填圖的地物類型，并用最終得分的高低來(lái)衡量匹配的可信度，分值越高匹配度越好。本研究使用兩種方法參與權(quán)重分析，分別是光譜角度匹配(SAM)和波譜特征擬合(SFF)，表1為L(zhǎng)andsat8 OLI數(shù)據(jù)得到的4種端元與實(shí)測(cè)光譜的識(shí)別結(jié)果。從表1中可以看出n-D Class #1與尾礦實(shí)測(cè)波譜的匹配程度最高達(dá)到1.618，因此，確定端元n-D Class #1為尾礦波譜。表2為珠海一號(hào)影像數(shù)據(jù)得到的4種端元與實(shí)測(cè)光譜的識(shí)別結(jié)果。從表2可以看出n-D Class #1與尾礦實(shí)測(cè)波譜的匹配程度最高達(dá)到2.069，因此，確定端元n-D Class #1為尾礦波譜。從綜合波譜匹配的結(jié)果來(lái)看珠海一號(hào)影像數(shù)據(jù)的匹配精度要優(yōu)于Landsat8 OLI影像。

表1 Landsat8 OLI影像端元波譜與實(shí)測(cè)波譜的匹配得分表Table 1 Matching score of end element spectrum and measured spectrum of Landsat8 OLI images

表2 珠海一號(hào)影像端元波譜與實(shí)測(cè)波譜的匹配得分表Table 2 Matching scores of end element spectrum and measured spectrum of Zhuhai-1 images

4.2 尾礦識(shí)別填圖

針對(duì)Landsat8 OLI影像數(shù)據(jù)進(jìn)行尾礦地物識(shí)別，確定影像分類閾值分別為0.10、0.13、0.15，圖7中深色區(qū)域?yàn)樽R(shí)別出的尾礦庫(kù)，其中閾值為0.10時(shí)，在非尾礦庫(kù)區(qū)域識(shí)別出的深色區(qū)域最少，即錯(cuò)分最少；閾值為0.15時(shí)，尾礦位置的識(shí)別面積最大，此時(shí)識(shí)別尾礦位置的像元數(shù)量最多。針對(duì)珠海一號(hào)影像數(shù)據(jù)識(shí)別的結(jié)果如圖8所示，最優(yōu)閾值分別為0.28、0.29、0.30，都取得了較為精準(zhǔn)的識(shí)別結(jié)果，基本覆蓋整個(gè)尾礦庫(kù)。對(duì)比Landsat8 OLI影像和珠海一號(hào)影像數(shù)據(jù)尾礦識(shí)別結(jié)果，提取的尾礦空間分布基本一致，珠海一號(hào)影像數(shù)據(jù)填圖結(jié)果優(yōu)于Landsat8 OLI的填圖結(jié)果。

圖7 Landsat8 OLI數(shù)據(jù)尾礦識(shí)別填圖Fig.7 Tailings identification mapping of Landsat8 OLI data

圖8 珠海一號(hào)數(shù)據(jù)尾礦識(shí)別填圖Fig.8 Tailings identification mapping of Zhuhai-1 data

5 結(jié) 論

本文對(duì)實(shí)測(cè)高光譜數(shù)據(jù)重采樣后，利用多源遙感影像進(jìn)行尾礦庫(kù)識(shí)別研究，得到結(jié)論如下所述。

1) 依據(jù)影像波段范圍對(duì)實(shí)測(cè)波譜進(jìn)行重采樣，在降低數(shù)據(jù)維數(shù)的同時(shí)，還能確保在相同波譜范圍開展波譜匹配實(shí)驗(yàn)。

2) 對(duì)比高光譜數(shù)據(jù)尾礦識(shí)別結(jié)果，Landsat8 OLI影像能夠有效識(shí)別尾礦信息且空間位置基本一致，從精度來(lái)說，高光譜影像識(shí)別尾礦庫(kù)的解譯效果更接近真實(shí)地物。