余先川，李建廣，2，徐金東，張立保，胡 丹

(1.北京師范大學(xué)信息科學(xué)與技術(shù)學(xué)院，北京 100875;2.中國傳媒大學(xué)信息工程學(xué)院，北京 100024)

0 引言

高光譜遙感圖像(簡稱高光譜圖像)是光譜分辨率在10 nm數(shù)量級范圍內(nèi)的圖像[1]。由于其光譜分辨率高，能夠提供更為豐富的地球表面信息，因此倍受國內(nèi)外學(xué)者的關(guān)注和應(yīng)用[2]。然而，高光譜圖像的空間分辨率普遍較低，混合像元現(xiàn)象嚴(yán)重，這是傳統(tǒng)的像元級圖像分類和面積量測精度難以達(dá)到使用要求的主要原因。為了提高高光譜圖像的分類精度，就必須解決混合像元的分解問題，使遙感應(yīng)用由像元級達(dá)到亞像元級。光譜解混的目的就是得到每個(gè)混合像元內(nèi)純地物(端元)的光譜及每種地物所占的比例(豐度)。

光譜解混模型從本質(zhì)上分為線性混合模型和非線性混合模型2種，區(qū)分二者的標(biāo)準(zhǔn)就是光子有沒有在地物間發(fā)生多次散射。通常宏觀尺度上的混合可認(rèn)為是線性的，而微觀尺度上的混合被認(rèn)為是非線性的[3]。目前，基于線形混合模型的光譜解混理論方法已得到了廣泛的研究與應(yīng)用，但由于受實(shí)際地物間復(fù)雜關(guān)系以及大氣散射的影響，光譜混合都是非線性的，這就使得應(yīng)用傳統(tǒng)的基于線性光譜混合模型的解混結(jié)果難以滿足精度要求[4－10]。

為了在一定程度上解決線性混合模型解混精度不高的問題，本文根據(jù)地物散射特點(diǎn)，提出了一種基于二次散射的非線性混合模型——二次散射模型(secondary scattering model，SSM)，并應(yīng)用該模型進(jìn)行了解混實(shí)驗(yàn)，對實(shí)驗(yàn)結(jié)果進(jìn)行了分析。

1 非線性混合模型

1.1 廣義的非線性混合模型

為了克服線性光譜混合模型的上述不足，許多學(xué)者利用非線性模型對高光譜圖像的光譜混合方式進(jìn)行了描述。

鑒于線性混合是非線性混合在光學(xué)多次散射被忽略情況下的特例[11]，本文定義了廣義的非線性光譜混合模型，其表達(dá)式為

式中:y=[y1，…，yL]T為 L波段的高光譜圖像;f是一未知的非線性函數(shù)向量;M=[m1，…，mR]為端元光譜矩陣;a=[a1，…，aR]T為端元的豐度向量(R為端元的個(gè)數(shù));n為0均值且方差為σ2的高斯噪聲向量。

由于f函數(shù)是未知的，因此基于非線性光譜混合模型的光譜解混是個(gè)病態(tài)問題，需要給出合適且確定的非線性混合模型才能進(jìn)行解混。為此，本文提出了基于二次散射的非線性混合模型即二次散射模型。

1.2 二次散射模型(SSM)

光子在進(jìn)入傳感器前一般會在地物間發(fā)生多次碰撞和散射，為了使計(jì)算簡化，本文忽略了3次以上的散射，只考慮光子在進(jìn)入傳感器前與地物發(fā)生2次撞擊的情況，提出二次散射模型，表達(dá)式為

1.3 模型推導(dǎo)

假設(shè)某一混合像元內(nèi)只存在A和B 2種地物，且地物的空間位置及比例都未知;設(shè)2種地物的光譜向量分別為m1和m2，豐度向量分別為a1和a2，即 M=[m1，m2]，a=[a1，a2]T。

光譜混合方式分線性和二次散射2大類，如圖1所示。

圖1 地物散射情況Fig.1 Ground object scattering case

從圖1可以看出，當(dāng)光子入射到某種地物上后，將以概率c直接反射到傳感器，以概率d(d=1－c)散射到其他地物上。以概率d散射的這部分光子又分別以概率a1和a2碰撞到地物A或B后再反射到傳感器中。該混合像元得到的光譜為

式中:前6項(xiàng)分別對應(yīng)圖1(a)—(f)6種情況;第7 項(xiàng)為誤差項(xiàng)。則式(3)可進(jìn)一步化簡為

1.4 二次散射系數(shù)的確定

光譜混合從形式上分為聚合式、整合式和致密式3種，如圖2所示。多次散射大多發(fā)生在致密式分布的地物中。

二次散射模型是半確定性模型，需要首先確定二次散射系數(shù)d。二次散射系數(shù)能直觀反映地物分布的復(fù)雜程度:地物分布越簡單，二次散射系數(shù)越低(極限情況d=0，為線性混合);相反，地物分布越復(fù)雜，非線性程度越強(qiáng)，則二次散射系數(shù)越高。

圖2 光譜混合的3種形式Fig.2 3 forms of spectrum mixture

假設(shè)某混合區(qū)由2種純地物組成，如果某個(gè)像元只存在其中一種地物，可認(rèn)為二次散射系數(shù)d=0(線性情況)。假設(shè)像元內(nèi)有2種地物且均等分布，則地物的交界地帶最多，可認(rèn)為此種情況下二次散射系數(shù)最大(圖3(a));如果2種地物不均等分布，比如地物A占10%，地物B占90%，則認(rèn)為此種情況地物間存在交界，但交界地帶相對較少，二次散射系數(shù)較低(圖3(b))。

圖3 地物分布復(fù)雜度Fig.3 Complexity of ground object distribution

地物豐度方差能表征地物的復(fù)雜程度，當(dāng)像元中只存在一種純地物，則豐度方差最大，對應(yīng)的地物復(fù)雜度最低，即二次散射系數(shù)最小(d=0);當(dāng)像元中地物分布均等時(shí)，則地物豐度方差為0，但地物分布最復(fù)雜，二次散射系數(shù)最大?？梢姡紊⑸湎禂?shù)與地物豐度方差成反比。

如果忽略地物組分的顆粒度、粗糙度及透明度等微觀因素，則二次散射系數(shù)可定義為

式中:d(i，j)∈[0，dmax]，dmax為最大的二次散射系數(shù)(本文設(shè)為0.5);var(i，j)為地物豐度方差，(i，j)為圖像坐標(biāo);varmax為各混合像元中地物豐度方差的最大值。

2 方法實(shí)驗(yàn)及結(jié)果分析

2.1 應(yīng)用模擬數(shù)據(jù)的實(shí)驗(yàn)及結(jié)果分析

2.1.1 模擬數(shù)據(jù)的生成

從ENVI USGS Spectral Libraries中找出天然堿、瀉利鹽和黝簾石3種地物光譜數(shù)據(jù)[12]。該光譜數(shù)據(jù)的波長區(qū)間為395～2 560 nm，光譜分辨率為10 nm，波段數(shù)為210個(gè)。分別根據(jù)光譜線性混合模型和二次散射模型生成大小為100像元×100像元且包含210個(gè)波段的高光譜圖像，并增加不同信噪比的高斯白噪聲。3種純地物的光譜曲線如圖4所示。

圖4 3種純地物的光譜曲線Fig.4 Spectral curves of 3 pure mineral

根據(jù)模擬數(shù)據(jù)生成的3種純地物豐度圖像如圖5所示。

圖5 模擬生成的3種純地物豐度圖Fig.5 Simulated abundance image of 3 minerals

將每種純地物豐度圖像變成一個(gè)1×10 000的行向量，然后將3個(gè)行向量組成一個(gè)3×10 000的地物豐度矩陣。3個(gè)地物的光譜數(shù)據(jù)矩陣為210行×3列，每一列對應(yīng)一種純地物的光譜向量，每種純地物光譜向量有210個(gè)采樣點(diǎn)。

首先，對模擬數(shù)據(jù)進(jìn)行光譜線型混合:將光譜矩陣與地物豐度矩陣相乘，得到210×10 000的混合矩陣。該矩陣表示高光譜圖像每1行為1個(gè)波段，共210個(gè)波段。圖6為第100和150波段的圖像，不同顏色代表不同反射率值。

圖6 線性混合模擬生成的高光譜圖像Fig.6 Simulated hyperspectral image based on linear mixture

根據(jù)式(1)對模擬數(shù)據(jù)進(jìn)行光譜非線性混合。圖7是二次散射系數(shù)d=0.3時(shí)模擬得到的第100和150波段的圖像。不同顏色代表不同反射率值。

圖7 非線性混合模擬生成的高光譜圖像Fig.7 Simulated hyperspectral image based on nonlinear mixture

在實(shí)驗(yàn)中，對模擬得到的數(shù)據(jù)加入不同信噪比的噪聲，信噪比計(jì)算公式為

式中var(s)和var(n)分別代表信號和噪聲方差。

圖8為非線性光譜混合后再加入50 dB噪聲后的第100波段圖像。

圖8 非線性光譜混合后加入50 dB噪聲的第100波段圖像Fig.8 100thband simulated hyperspectral image based on nonlinear mixture(RSNR=50 dB)

2.1.2 實(shí)驗(yàn)方法及結(jié)果

將模擬得到的數(shù)據(jù)分別用全約束最小二乘法(fully constraind least squares，F(xiàn)CLS)[13]和非線性變換后的最小二乘法進(jìn)行解混，得到估計(jì)的地物豐度數(shù)據(jù)。利用式(7)計(jì)算估計(jì)值的均方誤差[14]，即

式中:yn表示各純地物的實(shí)際豐度值向量;y^n表示各純地物的估計(jì)豐度值向量;N為采樣點(diǎn)個(gè)數(shù)，在本實(shí)驗(yàn)中N=60 000。

表1為在相同信噪比、不同二次散射系數(shù)情況下，基于線性和非線性模型解混得到的均方誤差。

表1 相同信噪比(RSNR=50)、不同散射系數(shù)下的解混均方誤差Tab.1 MSE of solved abundance at same RSNR(RSNR=50)but different scatter coefficient

由表1可以看出，隨著二次散射系數(shù)的增大，即地物分布復(fù)雜情況的加重，基于線性混合模型解混得到的均方誤差值增大?？梢?，地物分布越復(fù)雜，線性混合模型的適用性越差;而基于非線性混合模型的解混效果好得多。同時(shí)可以發(fā)現(xiàn)，隨著二次散射系數(shù)的提高，非線性解混效果也變得更好。

表2為在相同二次散射系數(shù)(d=0.3)、不同信噪比情況下，基于線性和非線性模型解混得到的均方誤差。

表2 相同二次散射系數(shù)(d=0.3)、不同信噪比下的解混均方誤差Fig.2 MSE of solved abundance at same d(d=0.3)but different RSNR

由表2可以看出，不管是基于線性混合模型還是基于非線性混合模型，信噪比RSNR越低，解混效果越差。與線性解混結(jié)果相比，非線性的解混效果更優(yōu)，提高的百分比隨著信噪比的增大而增大，表明非線性模型抗噪性更好。

圖9為在二次散射系數(shù)d=0.3、信噪比RSNR=40情況下，基于線性混合模型解混得到的地物A豐度圖像及其誤差圖像。

圖9 基于線性混合模型解混得到的地物A豐度圖像(左)及其誤差圖像(右)(d=0.3，RSNR=40)Fig.9 Solved abundance image(left)of mineral A and it’s error image(right)based on linear mixture(d=0.3，RSNR=40)

圖10 為在二次散射系數(shù) d=0.3、信噪比為RSNR=40情況下，基于非線性混合模型解混得到的地物A豐度圖像及其誤差圖像。

圖10 基于非線性混合模型解混得到的地物A豐度圖像(左)及其誤差圖像(右)(d=0.3，RSNR=40)Fig.10 Solved abundance image(left)of mineral A and it’s error image(right)based on nonlinear mixture(d=0.3，RSNR=40)

基于圖9，10可以發(fā)現(xiàn)，基于非線性混合模型的解混誤差明顯小于基于線性混合模型的解混誤差。

2.2 應(yīng)用AVIRIS數(shù)據(jù)的實(shí)驗(yàn)及結(jié)果分析

2.2.1 數(shù)據(jù)的選取

選取美國內(nèi)華達(dá)州南部Cuprite地區(qū)的AVIRIS高光譜數(shù)據(jù)[15]。圖像大小為250像元×191像元，包含224個(gè)波段，波長區(qū)間為0.389～2.467μm。除去低信噪比和水氣吸收波段(波段1～2，104～113，148～167，221～224)，最后可用數(shù)據(jù)為188個(gè)波段數(shù)據(jù)。圖11為第100波段高光譜圖像。

圖11 Cuprite地區(qū)AVIRIS第100波段圖像Fig.11 100thband image of AVIRIS of Cuprite district

Cuprite區(qū)主要分布6種礦物[16]，分別為Desert-Varnish(沙漠巖漆)、Alunite(明礬石)、Montmorillonite(蒙脫石)、Chalcedony(玉髓)、Kaolinite(高嶺石)#1和Kaolinite(高嶺石)#2。

2.2.2 實(shí)驗(yàn)方法及結(jié)果

利用頂點(diǎn)成分分析算法[17]提取6種礦物端元，基于線性混合模型、利用FCLS方法[13]得到的6種礦物豐度圖像如圖12所示。

圖12－1 Cuprite地區(qū)6種礦物線性解混結(jié)果Fig.12－1 Solved abundance results of 6 minerals based on linear mixture in Cuprite district

圖12－2 Cuprite地區(qū)6種礦物線性解混結(jié)果Fig.12－2 Solved abundance results of 6 minerals based on linear mixture in Cuprite district

利用式(5)確定基于該解混結(jié)果的地物復(fù)雜度，設(shè)定最大二次散射系數(shù)dmax=0.5，得到的地物復(fù)雜度分布如圖13所示。

基于非線性混合模型(式(1))，使用上面得到的二次散射系數(shù)對每個(gè)混合像元進(jìn)行解混，得到的解混結(jié)果如圖14所示。

對比圖12，14可以看出，基于二次散射模型的解混效果更好，圖像上道路影像特征清晰(圖14);而基于線性模型解混的道路影像特征則模糊(圖12)。主要原因是道路與周圍地物是主要邊界地帶，而邊界地帶容易發(fā)生二次散射。

圖13 Cuprite地區(qū)地物分布復(fù)雜度圖Fig.13 Complexity of ground object distribution in Cuprite district

圖14 Cuprite地區(qū)6種礦物非線性解混結(jié)果Fig.14 Solved abundance results of 6 minerals based on nonlinear mixture in Cuprite district

圖15 顯示了在不同散射系數(shù)下(d=0～0.5)， 蒙脫石的分布情況。

圖15 不同二次散射系數(shù)下的蒙脫石分布情況Fig.15 Solved abundance of Montmorillonite in different secondary scatter coefficient

從圖15可以發(fā)現(xiàn)，基于二次散射模型解混的道路影像特征更加清晰連續(xù)，而基于線性混合模型(d=0)的道路影像特征則比較模糊，并出現(xiàn)道路影像局部中斷現(xiàn)象。

3 結(jié)論

本文對高光譜遙感圖像的光譜解混問題進(jìn)行了研究，得出以下結(jié)論:

1)提出了廣義的非線性混合模型。

2)針對線性混合模型解混精度不高的問題，提出了一種基于二次散射的非線性混合模型，即二次散射模型(SSM)。

3)基于地物豐度方差，定義了地物混合復(fù)雜度。

4)通過計(jì)算解混后的地物豐度與實(shí)際地物豐度的均方誤差發(fā)現(xiàn)，與基于線性混合模型解混結(jié)果相比，基于二次散射模型的解混精度能提高近一個(gè)數(shù)量級。

本論文豐富了高光譜遙感圖像的光譜解混理論，為光譜解混的進(jìn)一步應(yīng)用提供了理論基礎(chǔ)。今后將針對不同區(qū)域地形地貌特點(diǎn)，對二次散射模型的二次散射系數(shù)做進(jìn)一步研究，以期得到更精確的解混結(jié)果。

[1]Schweizer S M，Moura J M F.Hyperspectral imagery:Clutter adaptation in anomaly detection[J].IEEE Tran on Information Theory，2000，46(5):1855－1871.

[2]童慶禧，張 兵，鄭蘭芬.高光譜遙感——原理、技術(shù)與應(yīng)用[M].北京:高等教育出版社，2006:166－169.Tong Q X，Zhang B，Zheng L F.Hyperspectral remote sensing[M].Beijing:Higher Education Press，2006:166－169.

[3]Keshava N，Mustard J.Spectral unmixing[J].IEEE Signal Processing Mag，2002，19(1):44－57.

[4]Craig M.Minimum volume transforms for remotely sensed data[J].IEEE Trans Geosci and Remote Sensing，1994，32(3):542－552.

[5]Nascimento J M P，Bioucas－ Dias J M.Does independent component analysis play a role in unmixing hyperspectral data?[J].IEEE Trans Geosci and Remote Sensing，2005，43(1):175－187.

[6]Jia S，Qian Y.Constrained nonnegative matrix factorization for hyperspectral unmixing[J].IEEE Trans Geosci and Remote Sensing，2009，47(1):161－173.

[7]Zortea M，Plaza A.Spatial preprocessing for endmember extraction[J].IEEE Trans Geosci and Remote Sensing，2009，47(8):2679－2693.

[8]Ghosh G，Kumar S，Saha S K.Hyperspectral satellite data in mapping salt－affected soils using linear spectral unmixing analysis[J].Journal of the Indian Society of Remote Sensing，2012，40(1):129－136.

[9]Liu K H，Wong E，Du E Y，et al.Kernel－ based linear spectral mixture analysis [J].IEEE Trans Geosci and Remote Sensing，2012，9(1):129－133.

[10]Deng Y B，F(xiàn)an F L，Chen R R.Extraction and analysis of impervious surfaces based on a spectral unmixing method using Pearl River Delta of China landsat TM/ETM+imagery from 1998 to 2008[J].Sensors，2012，12(2):1846－1862.

[11]Iordache M D，Bioucas－Dias J M.Sparse unmixing of hyperspectral data[J].IEEE Trans Geosci and Remote Sensing，2011，49(6):2014－2039.

[12]RSI(Research Systems Inc.).ENVI user’s guide version 4.0[S].Boulder，CO 80301 USA，2003.

[13]Heinz D C，Chang C I.Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery[J].IEEE Tran GRS，2001，39(3):529－ 545.

[14]Plaza J，Plaza A，Martinez P，et al.Nonlinear mixture models for analyzing laboratory simulated－ forest hyperspectral data[J].Proc SPIE Image and Signal Processing for Remote Sensing IX，2004，52(38):480－487.

[15]AVIRIS data of Cuprite distuict，America.http://aviris.jpl.nasa.gov/html/aviris.freedata.html.

[16]Vane G，Green R，Chrien T G，et al.The airborne visible/infrared imaging spectrometer(AVIRIS)[J].Remote Sensing of the Environment，1993，44(2/3):127－ 143.

[17]Nascimento J，Dias J B.Vertex component analysis:A fast algorithm to unmix hyperspectral data[J].IEEE Trans Geosci and Remote Sens，2005，43(4):175－187.