基于实测光谱和SRF的稀疏植被区MODIS积雪信息提取

刘 艳，李 杨，张 璞

(1．中国气象局乌鲁木齐沙漠气象研究所，乌鲁木齐 830002;2．乌鲁木齐气象卫星地面站，乌鲁木齐 830011)

0 引言

地物光谱数据可分为在实验室非常严格条件下获取的材料光谱(如矿物)、在野外相对均一条件下测量得到的端元光谱(如积雪)及由遥感平台上的传感器测量得到的像元光谱(像元对应地物光谱信号的综合)［1－2］。由于传感器空间分辨率的限制和被探测目标的多样性，在获取的遥感图像中往往会产生混合像元，影响图像分析的精度［3］。国内外许多学者提出了几何光学模型、线性及非线性光谱模型、随机几何模型等混合像元分解模型，其中，对线性光谱混合模型(linear spectral mixture model，LSMM)研究最多、应用最广［4－9］。应用 LSMM 时，端元(组成混合像元的纯像元)类型、数量及端元光谱值综合决定了混合像元分解的精度［10］。在已有研究中，端元数目的确定和端元选择基本上是利用星载或机载遥感数据、基于像元纯度指数(pixel purity index，PPI)等算法从遥感图像中提取的，由于图像空间分辨率低，端元提取精度不高［11］。使用Field Spec Pro FR(full range)全波段光谱仪(analytical spectral devices，ASD)采集的反射光谱曲线因波段间隔非常小(在350～1 000 nm谱段，采样间隔1．4 nm，光谱分辨率3 nm;在1 000～2 500 nm谱段，采样间隔2 nm，光谱分辨率10 nm)，可被近似地看作是一条连续变化曲线。要完成端元波谱向像元波谱转换、直接获取混合像元内高精度的端元光谱，多光谱成像传感器光谱响应函数(spectral response function，SRF)是一个可利用的重要手段。SRF是波长的函数，是在每一波长上传感器记录辐射能量与入射辐射能量的比值［2，12］。

许多学者以MODIS遥感数据为基础，对我国新疆及天山南、北地区进行积雪、冰川及融雪径流研究［13－15］，其积雪面积提取结果以像元为单位，很少考虑到亚像元级中的混合像元问题。陈晓娜等［16］以MOD 02 HKM数据为基础，通过LSMM对新疆天山中段MODIS图像进行混合像元分解，从中提取了积雪面积;延昊等［17］使用混合像元分解法，根据甚高分辨率辐射计(advanced very high resolution radiometer，AVHRR)数据反演了积雪盖度;Painter等［18］利用机载可见光/红外成像光谱仪(airborne visible infrared imaging spectrometer，AVIRIS)图像结合离散纵标辐射传输(discrete ordinates radiative transfer，DISORT)模型，并考虑到积雪亚像元对雪粒径反演的影响，发展了多端元雪覆盖面积和雪粒大小(multiple endmember snow covered area and grain size，MEMSCAG)模型同步反演雪粒径，取得了较好的效果。但针对稀疏植被区积雪覆盖的亚像元级遥感精确估算方法研究目前还处于空白。本文选用MODIS作为参考传感器，在稀疏植被区野外观测样方内使用ASD准同步实测积雪端元光谱数据，根据MODIS1—7波段的SRF，从能量角度完成了端元光谱到像元光谱的转换，提取出高精度端元光谱数据;将其与MODIS影像端元光谱进行比较，优选影像端元光谱值并将其应用到LSMM混合像元分解中，完成了稀疏植被区积雪遥感信息的提取。

1 研究方法

1．1 遥感图像及预处理

本文选用的遥感图像为2008年12月20日获取的、覆盖新疆准噶尔盆地古尔班通古特沙漠的中等分辨率成像光谱仪(MODIS 1B)数据，空间分辨率为500 m，包含7个波段(MODIS1—7)。对 MODIS 1B数据进行了 FLAASH(fast line－of－sight atmospheric analysis of spectral hypercube)大气校正、辐射校正和几何纠正等预处理。图像中有一定量的积雪覆盖，适用于积雪遥感信息提取研究。

1．2 试验设计

在古尔班通古特沙漠建立了积雪参数野外采集场，布设360 m×360 m的采集样方，并分割成16个90 m×90 m的样方用于MODIS准同步观测(图1(a))。采集场内地势平坦，沙漠地表植被主要为梭梭，高度约为20～30 cm，盖度约为30%(图1(b))。2008年12月2—3日新疆北疆地区出现较大降雪天气过程，积雪参数野外采集期(2008年12月20日)雪深6 cm，积雪覆盖率在90%以上(图1(c))。

图1 MODIS准同步试验场和地表类型图Fig．1 Maps of quasi－synchronism test ground for MODIS and the surface types

1．3 光谱测定

2008年12月20日，天气晴朗无风，使用ASD光谱仪(工作谱段为350～2 500 nm)在样方内进行了反射光谱准同步观测，观测时间为MODIS过顶前后5 min。采集得到梭梭端元光谱曲线6条、积雪端元光谱曲线15条(图2)。测量时ASD光谱仪的传感器探头垂直于观测目标，距离目标30 cm，所测量的地物光谱反射率经专用平面白板标准化。对每个测量点位进行10次地物反射光谱测量，最后取其均值作为观测目标的反射率。

图2 积雪参数采集样方内观测端元(左)及其光谱曲线(右)Fig．2 End members in acquisition quadrats for snow parameter(left)and their spectral curves(right)

1．4 端元光谱处理

传感器的每个波段都有一定宽度的光谱响应范围，但光谱响应范围内的入射辐亮度不是完全被传感器接收并记录的。在理论上，每个波段实际接收到的能量是该波段范围内各波长处接收能量的总和［3］。将野外测量端元光谱转换为遥感图像像元光谱，相当于对实测光谱按波长做加权平均，每个波长处的权重为，其转换模型［2，12］为

式中:Ri为端元光谱R(λ)经光谱响应函数fi(λ)校正后得到的第i个波段的像元光谱反射率值;R(λ)为端元各波长处的反射率值;fi(λ)为MODIS传感器光谱响应函数;λmax和λmin分别为端元光谱取值范围上界和下界。根据野外采集端元光谱(图3)和MODIS各波段的SRF(图4)，利用式(1)计算得到MODIS1—7波段的像元光谱(图5)。

图3 ASD采集的积雪反射光谱Fig．3 Reflectance spectra of snow collected with ASD

图4 MODIS1—7波段光谱响应函数Fig．4 SRF of MODIS 1—7 bands

图5 野外采集积雪端元光谱转换得到的MODIS1—7波段像元光谱Fig．5 Pixel spectra of MODIS1—7 bands converted from the snow end members collected in the field

1．5 LSMM 分析

LSMM中像元在某一波段的反射率是以构成该像元的端元反射率和端元所占像元面积比例为权重的线性组合，其计算公式［10］为

式(2)—(4)中:i为光谱波段数，i=1，2，…，m;k为端元数，k=1，2，…，n;Ri为第 i波段的像元反射率(已知)，包含一种或多种端元;fk为1个像元内第k端元所占分量(待求);cik为1个像元内端元k在第i波段上的反射率;ERi为第i波段的误差;n为端元数;m为可用波段数(n≤m+1)。

各影像像元中端元所占像元面积比例可以通过最小二乘法来得到。LSMM的基本思想是求取使误差ERi最小的解。通过引入不同的约束条件，得到相应的LSMM分解策略:①若不考虑约束条件式(3)和(4)，仅使误差最小，是纯粹的最小二乘求解过程，称此策略为无约束条件的LSMM;②若满足约束条件式(3)和(4)，则使误差最小的策略被称之为带全约束条件的LSMM。

1．6 最优MODIS影像端元选择

混合像元内影像端元的光谱值可以根据实地测量光谱确定或直接从光谱数据库获取，也可以从遥感图像分析中获得。本文采用最小噪声分离(minimum noise fraction，MNF)变换和像元纯度指数(PPI)，结合目视解译选择端元［19－20］。通过对 MNF变换特征值的分析可知，端元信息主要集中在MNF变换后的前2个分量(所占比例分别为69%和24%)。结合实验区的实地调查，分别解译出云、积雪和裸土3种端元(图6)。

图6 选取的MODIS影像端元光谱Fig．6 End member spectral values selected from MODIS image

将用上述方法得到的最优MODIS影像端元光谱与野外采集端元光谱经转换得到的MODIS1—7波段积雪像元光谱进行比较(图7)，发现所提取影像端元在MODIS1波段的光谱值接近1，远大于转换光谱值(约为0．65);在MODIS2—7波段，影像端元光谱值与转换值趋于接近。

图7 实测端元与提取的影像端元光谱比较Fig．7 Comparison between the end member spectral values measured in the field and extracted from MODIS image

因此，最终确定选取MODIS2—7波段影像端元光谱值(图8)输入LSMM，进行像元分解。

图8 确定用于LSMM的端元光谱Fig．8 Determined end member spectrum used for LSMM

2 结果分析

2．1 LSMM光谱分解结果

图9为 MODIS6(R)，4(G)，3(B)假彩色合成图像(蓝色区域为积雪覆盖区域)。

图9 研究区MODIS 6(R)，4(G)，3(B)假彩色合成图像Fig．9 MODIS 6(R)，4(G)，3(B)false color composite image of the study area

将确定的积雪、裸土和云端元输入LSMM进行像元分解，得到积雪、裸土和云分量值(图10)。

图10 积雪、裸土和云分量值图Fig．10 Component diagrams of snow，bare soil and cloud

对比图9和图10可以看出，积雪、裸土和云端元分量值的分布与它们的实际分布范围吻合，积雪分量值较大的区域(图10(a))对应实际积雪覆盖较为丰富的区域(图9);同时，积雪分量值较大的区域(图10(a))对应的裸土分量值也相应较小(图10(b))。

在观测样方内，全约束条件LSMM分解法获取的积雪分量(图11(a))和云分量(图11(b))分布范围符合实际观测结果，全约束条件LSMM比无约束LSMM在稀疏植被区积雪亚象元分解(图11(c)和(d))中更加有效。

图11 观测样方内积雪和云分量值Fig．11 Component values of snow and cloud in the acquisition quad rats

2．2 端元优选比较

将应用MODIS1—7波段端元LSMM分解所得积雪端元分量值与归一化差分积雪指数(normalized difference snow index，NDSI)进行相关分析(图12(a))，发现其相关系数(R=0．28)远低于应用MODIS2—7波段端元无约束LSMM分解所得积雪端元分量值与NDSI相关分析(图12(b))的相关系数(R=0．47)和全约束条件LSMM分解所得积雪端元分量值与NDSI相关分析(图12(c))的相关系数(R=0．92)。

图12 MODIS端元LSMM分解积雪端元分量值与NDSI的相关性Fig．12 Correlation between snow component values using LSMM for MODIS and NDSI

图12表明，剔除MODIS1波段后LSMM分解的结果较理想，其中全约束条件LSMM所提取的积雪分量与NDSI的相关性最高，故全约束条件LSMM所提取的积雪分量可以作为估算积雪的典型指数。

3 结论

1)针对新疆北疆地区冬季多云特点和稀疏植被区积雪遥感信息多为混合像元的现状，结合ASD实测光谱和MODIS传感器光谱响应函数(SRF)，从能量角度对野外采集端元光谱进行拟合，将拟合离散光谱值与MODIS影像端元光谱值进行比较，发现经FLASSH大气校正后MODIS1波段的积雪反射值接近1，存在过饱和现象，远大于拟合离散光谱值。

2)在LSMM分解过程中，分别使用最小噪声分离(MNF)变换和像元纯度指数(PPI)法获取了对应于MODIS1—7和MODIS2—7波段的积雪、裸土和云3类端元。结果表明，剔除MODIS1波段后线性光谱混合模型(LSMM)分解的结果较理想，提取的积雪分量值与NDSI的相关性显著提高，无约束LSMM所得积雪端元分量值与NDSI相关系数为0．47，全约束条件 LSMM所得积雪端元分量值与NDSI相关系数高达0．92，全约束条件LSMM所提取的积雪分量可以作为估算积雪的典型指数。

［1］ 王锦地，李小文，张立新，等．中国典型地物标准波谱数据库建设［J］．遥感学报，2003，7(sl):13－16．Wang JD，Li X W，Zhang L X，et al．Database construction for standard spectral of typical objects in China［J］．Journal of Remote Sensing，2003，7(sl):13－16．

［2］ 张 婷，丁建丽，王 飞．基于实测端元光谱的多光谱图像光谱模拟研究［J］．光谱学与光谱分析，2010，30(11):2889－2892．Zhang T，Ding J L，Wang F．Simulation of image multi－spectral using field measured end member spectrum［J］．Spectroscopy and Spectral Analysis，2010，30(11):2889－2892．

［3］ 潘明忠，元洪兴，肖功海，等．基于高精度端元的混合像元线性分解模型研究［J］．红外与毫米波学报，2010，29(5):357－361．Pan M Z，Yuan H X，Xiao G M，et al．Study of mixed pixel linear unmixing model based on high accurate end member［J］．Journal of Infrared and Millimeter Waves，2010，29(5):357－361．

［4］ Chabrillat S，Pinet P C，Ceuleneer G，et al．Ronda peridotite massif:Methodology for its geological mapping and lithological discrimination from airborne hyper spectral data［J］．International Journal of Remote Sensing，2000，21(12):2363－2388．

［5］ 李 慧，陈健飞，余 明．线性光谱混合模型的ASTER影像植被应用分析［J］．地球信息科学，2005，7(1):103－106．Li H，Chen J F，Yu M．Auto－detection of land cover changes based satellite image fusion technique［J］．Geo－information Science，2005，7(1):103－106．

［6］ Small C．Estimation of urban vegetation abundance by spectral mixture analysis［J］．International Journal of Remote Sensing，2001，22(7):1305－1334．

［7］ 张熙川，赵英时．应用线性光谱混合模型快速评价土地退化的方法研究［J］．中国科学院研究生院学报，1999，16(2):170－172．Zhang X H，Zhao Y S．Application of line spectral mixture model to rapid assessment of land degradation in semiarid area［J］．Journal of Graduate School，Academia Sinica，1999，16(2):169－176．

［8］ 万 军，蔡运龙．应用线性光谱分离技术研究喀斯特地区土地覆被变化——以贵州省关岭县为例［J］．地理研究，2003，22(4):440－443．Wan J，Cai Y L．Applying linear spectral unmixing approach to the research of land cover change in karst area:A case in Guanling county of Guizhou Province［J］．Geographical Research，2003，22(4):440－443．

［9］ 邹 蒲，王云鹏，王志石，等．基于ETM+图像的混合像元线性分解方法在澳门植被信息提取中的应用及效果评价［J］．华南师范大学学报:自然科学版，2007(7):131－136．Zou P，Wang Y P，Wang Z S，etal．Accessing the linear spectral un－mixing approach for extracting vegetation information using Landsat ETM+data in Macao［J］．Journal of South China Normal University:Natural Science Edition，2007(7):131－136．

［10］ 饶 萍．EOS－MODIS像元组分分解中端元的选择与改进［C］//《测绘通报》测绘科学前沿技术论坛摘要集，2008:1－10．Rao P．Choice and improvement of end members in EOS－MODIS pixel component unmixing［C］//Bulletin of Surveying and Mapping．Abstracts of Forum on Advanced Surveying and Mapping，2008:1－10．

［11］ 高晓惠，相里斌，魏儒义，等．基于光谱分类的端元提取算法研究［J］．光谱学与光谱分析，2011，31(7):1995－1998．Gao X H，Xiang L B，Wei R Y，et al．Research on end member extraction algorithm based on spectral classification［J］．Spectroscopy and Spectral Analysis，2011，31(7):1995－1998．

［12］ 王福民，黄敬峰，唐延林，等．采用不同光谱波段宽度的归一化植被指数估算水稻叶面积指数［J］．应用生态学报，2007，18(11):2444－2450．Wang FM，Huang H F，Tang Y L，et al．Estimation of rice LAI using NDVI at different spectral bandwidths［J］．Chinese Journal of Applied Ecology，2007，18(11):2444－2450．

［13］ 冯琦胜，张学通，梁天刚．基于MOD10A1和AMSR－E的北疆牧区积雪动态监测研究［J］．草业学报，2009，18(1):125－133．Feng Q S，Zhang X T，Liang T G．Dynamic monitoring of snow cover based on MOD10A1 and AMSR－E in the north of Xinjiang Province，China［J］．Acta Prataculturae Sinica，2009，18(1):125－133．

［14］ 裴 欢，房世峰，覃志豪，等．基于遥感的新疆北疆积雪盖度及雪深监测［J］．自然灾害学报，2008，17(5):52－57．Pei H，Fang SF，Jia Z H，et al．Remote sensing－based monitoring of coverage and depth of snow in northern Xinjiang［J］．Journal of Natural Disasters，2008，17(5):52－57．

［15］ 李宝林，张一驰，周成虎．天山开都河流域雪盖消融曲线研究［J］．资源科学，2004，26(6):23－29．Li B L，Zhang Y C，Zhou C H．Snow cover depletion curve in Kaidu River basin，Tianshan Mountains［J］．Resources Science，2004，26(6):23－29．

［16］ 陈晓娜，包安明，张红利，等．基于混合像元分解的MODIS积雪面积信息提取及其精度评价——以天山中段为例［J］．应用气象学报，2004，15(6):665－671．Chen X N，Bao A M，Zhang H L，et al．A study on methods and accuracy assessment for extracting snow covered areas from MODIS images based on pixel unmixing:A case on the middle of the Tianshan Mountain ［J］．Journal of Applied Meteorological Science，2004，15(6):665－671．

［17］ 延 昊，张国平．混合像元分解法提取积雪盖度［J］．应用气象学报，2004，15(6):655－671．Yan H，Zhang G P．Unmixing method applied to NOAA－AVHRR data for snow cover estimation［J］．Journal of Applied Meteorology，2004，15(6):655－671．

［18］ Painter T H，Dozier J，Roberts D A，et al．Retrieval of subpixel snow－covered area and grain size from imaging spectrometer data［J］．Remote Sensing of Environment，2003，85(1):64－77．

［19］ 崔耀平，刘 彤，赵志平，等．干旱荒漠区植被覆盖变化的遥感监测分析［J］．地理信息学报，2011(3):305－311．Cui Y P，Liu T，Zhao Z P，et al．Using multi spectral remote sensing data to extract and analyze the vegetation change of the western Gurbantunggut desert［J］．地理信息学报，2011(3):305－311．

［20］ 崔耀平，王让会，刘 彤，等．基于光谱混合分析的干旱荒漠区植被遥感信息提取研究——以古尔班通古特沙漠西缘为例［J］．中国沙漠，2010，30(2):335－341．Cui Y P，Wang R H，Liu T，et al．Extraction of vegetation information in arid desert area based on spectral mixture analysis［J］．Journal of Desert Research，2010，30(2):335－341．