一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法

摘要

一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法，所述方法包括对阴影区进行识别和提取、非阴影相似像元搜寻以及阴影区光谱信息恢复步骤。本发明所提供的一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法，用于进行遥感图像阴影区的光谱信息恢复，该方法充分结合了阴影像元本身的弱光谱信息及其邻近相似像元的光谱信息，摆脱了传统地形纠正算法对DEM数据的依赖，当DEM数据的分辨率和精度不能满足地形纠正算法得要求时，该算法可以提供一种行之有效的阴影处理方法。

说明书

技术领域

本发明涉及一种对Landsat TM/ETM+图像进行处理的方法，特别是一种 Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法。

背景技术

一直以来，在遥感成像领域，Landsat TM/ETM+图像有着大量的应用，美 国陆地卫星历经三代，第一代是Landsat-1卫星、Landsat-2卫星和Landsat-3 卫星，第二代是Landsat-4卫星和Landsat-5卫星，第三代是Landsat-6卫 星和Landsat-7卫星。TM是指美国陆地卫星4～5号专题制图仪Thematic  Mapper所获取的多波段扫描影像，TM图像有7个波段，其波谱范围是：TM-1 为0.45～0.52微米，TM-2为0.52～0.60微米，TM-3为0.63～0.69微米， 以上为可见光波段；TM-4为0.76～0.90微米，为近红外波段；TM-5为1.55～ 1.75微米，TM-7为2.08～2.35微米，为中红外波段；TM-6为10.40～12.50 微米，为热红外波段。因TM影像具较高空间分辨率、波谱分辨率、极为丰富 的信息量和较高定位精度，从20世纪80年代中后期开始成为世界各国广泛 应用的重要的地球资源与环境遥感数据源。

LANDSAT-7卫星的改进型主题测绘仪Enhanced Thematic Mapper Plus， 也就是ETM+，是在TM基础上改进的。ETM+相对TM增加了1个金色谱段和2 个增益区域，增加了太阳定标器，并提高了红外谱段的分辨率。

在Landsat TM/ETM+图像中，阴影现象普遍存在，阴影由高出地面的物 体遮挡太阳直射光而形成，可以分为投射阴影，即是高物体遮挡光源而投射 到地面的阴影，和自身阴影，即是物体自身未受到光源直射部分的阴影。可 以参见Salvador，E.，Cavallaro，A.，& Ebrahimi，T.(2001).Shadow  identification and classification using invariant color models.In  IEEE International Conference on Acoustics，Speech and Signal Processing，3，1545-1548.以及Yao，J.，& Zhang，Z.(2006). Hierarchical shadow detection for color aerial images.Computer Vision  and Image Understanding，102，60-69.文献描述。在估测建筑物的形状和 高度等三维结构重建方面，阴影有着积极的意义，然而在大部分图像处理和 工程应用中，由于遥感图像中阴影区DN值相对偏低，光谱信息量相对较少， 因此给地物类型的解译和地表覆盖分类等工作带来了很大的困难。尤其在山 地区域，复杂崎岖的地形造成了山体阳坡和阴坡之间接收直射光照量的巨大 差异，阴影现象尤为严重。

地形纠正是解决阴影问题的一种有效手段。在过去数十年间，研究者们 提出了很多地形纠正算法：针对朗伯表面，Teillet提出了余弦法，C纠正法 和经验统计模型；参见Teillet，P.M.，Guindon，B.，& Goodenough，D.G. (1982).On the slope-aspect correction of multispectral scanner data. Canadian Journal of Remote Sensing，8，84-106.。针对非郎伯表面，应 用最广泛的是Minnaert纠正法，参见Minnaert，M.(1941).The reciprocity  principle in lunar photometry.Astrophysics Journal，93，403-410.。 随后Richter考虑了不同地物类型的影响，将此方法发展为新Minnaert纠正 法，参见Richter，R.，Kellenberger，T.，& Kaufmann，H.(2009). Comparison of topographic correction methods.Remote Sensing，1， 184-196.。除了上述基于几何光学模型的算法之外，Dozier等基于大气和 地表间的辐射传输模型提出的IPW算法也是一种很有效的地形纠正方法，参 见Dozier，J.，& Frew，J.(1981).Atmospheric corrections to satellite  radiometric data over rugged terrain.Remote Sensing of Environment， 11，191-205.；Dozier，J.(1989).Spectral signature of alpine snow  cover from the Landsat Thematic Mapper.Remote Sensing of Environment， 28，9-22.；Dozier，J.& Frew，J.(1990).Rapid calculation of terrain  parameters for radiation modeling from digital elevation data.IEEE  Transactions on Geoscience and Remote Sensing，28，963-969。同时， 近些年也有很多更为高级的纠正算法先后被提出。然而，这些算法大部分都 需要数字高程模型数据，即DEM数据的辅助，但获取符合这些算法要求的， 具有合适的空间分辨率和较高的精度的DEM数据比较困难。例如，在对 Landsat TM/ETM+图像进行地形纠正时，需要空间分辨率为30米的DEM数据， 而目前尚无全球范围内的30m分辨率DEM数据；而且在大多数情况下，即使 我们拥有足够的DEM数据能够覆盖我们的研究区域，其分辨率也往往与待处 理的遥感图像并不匹配。因此，在进行DEM数据重采样以及DEM数据与遥感 图像数据几何配准等预处理过程中，则会无可避免地引入额外的噪声和误差， 从而损坏地形纠正的结果，导致纠正后图像呈现破碎和条纹等现象。

发明内容

本发明提供了一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复 方法，该方法可减少或避免前面所提到的问题。

为解决上述问题，本发明提出了一种Landsat TM/ETM+图像中山体阴影 区的光谱信息恢复方法，其包括如下步骤，

A、对阴影区进行识别和提取：该步骤利用Landsat TM/ETM+图像中阴影 区和非阴影区在表观亮度上的明显差异来对二者进行识别区分，从而对山体 阴影区进行提取，

首先，对需要进行处理的图像进行缨帽变换，并将变换后所得的图像的 各个像元的亮度波段数值作为每一个像元的亮度值，然后根据所述的每个像 元的亮度值对所述变换后所得的图像进行分割，将所述变换后所得的图像划 分为多个对象区域，所述对象区域内部的像元具有相似的空间性质和光谱性 质；

计算出每个对象区域内部像元的亮度平均值，将所述亮度平均值作为该 对象区域的亮度值，然后利用自适应阈值算法得到阈值，根据所述阈值对阴 影区进行提取；

B、非阴影相似像元搜寻：该步骤根据图像分割以及阴影区提取的结果， 进行相似像元的寻找，对阴影区的每一个阴影像元，利用用于光谱吸收特征 分析的亮度标准化技术“连续统去除”，根据像元光谱曲线的形状特点信息， 在其邻近空间范围内搜索出一定数量的与所述阴影像元具有相似的光谱曲线 形状特点的非阴影相似像元，以便后继用于对阴影像元光谱信息的恢复。

在图像分割后，阴影区以对象区域的形式被提取出来，对于每一个阴影 对象区域，首先在其周围建立至少两个像元宽度的缓冲区，然后在此缓冲区 范围内，对应所述对象区域中每个目标阴影像元选出与所述目标阴影像元光 谱相似度最高的至少N个非阴影像元，作为所述目标阴影像元的非阴影相似 像元，此处，相似度由像元间光谱曲线的标准偏差RMSD来衡量：

$\mathrm{RMS}{D}_i=\sqrt{\frac{\underset{b=1}{\overset{n}{\Sigma }}{(\mathrm{CR}(x_i,y_i,b)-\mathrm{CR}(x,y,b))}^2}{n}}$

其中CR(xi，yi，b)是位于(xi，yi)处的像元在第b波段的连续统去 除信息CR值，CR(x，y，b)则是目标阴影像元在第b波段的连续统去除信 息CR值，n为总的波段数。RMSD越小，则像元间的相似度越高，反之则越低。

被选择作为所述目标阴影像元的所述非阴影相似像元满足下列条件

${\mathrm{RMSD}}_i\le [\underset{b=1}{\overset{n}{\Sigma }}\sigma \left(b\right)\times 2/m]/n$

其中σ(b)表示整个图像在第b波段的连续统去除信息CR值的标准差，m 表示图像中的地物类别数。

如果在两个像元宽度的缓冲区范围内，满足条件的所述非阴影相似像元 数量个数N无法满足大于等于20的要求，则以两个像元为步长扩充所述缓冲 区的宽度，直到此要求被满足为止。

C、阴影区光谱信息恢复：该步骤通过利用所述非阴影相似像元提供的平 均亮度信息替换所述阴影像元的亮度信息，从而实现对阴所述影像元的光谱 信息进行反演，所述阴影像元的亮度信息来自所述非阴影相似像元连续统去 除后得到的连续统信息CI，经此处理后，所述阴影区像元的光谱特征得到恢 复，其表观亮度也将被提高到与其周围非阴影像元相近的水平，从而使得所 述阴影区的形状特征得以表现，

对于每个非阴影相似像元，根据该像元与对应的所述阴影像元的空间距 离以及光谱相似度，将二者分别标准化后相结合，即可得出表示所述非阴影 相似像元对用于进行光谱恢复的平均亮度值的贡献程度的权重值W，根据所 述非阴影相似像元的所述权重值W，即可计算出所述目标阴影像元的对应的 非阴影相似像元亮度的加权平均值CIwavg(b)，所述加权平均值CIwavg(b) 的计算公式如下

${\mathrm{CI}}_{\mathrm{wavg}\left(b\right)}=\underset{j=1}{\overset{N}{\Sigma }}{W}_j\times {\mathrm{CI}}_{j\left(b\right)}$

其中权重值Wj表示第j个非阴影相似像元对平均亮度值的贡献程度， CIj(b)表示第j个非阴影相似像元在第b波段的连续统信息，根据所述加权 平均值CIwavg(b)即可得到所述目标阴影像元在第b波段的经过恢复后的DN 值，计算公式如下，

OI_re(b)＝CI_wavg(b)×CR_(b)

其中OI_re(b)表示目标阴影像元在第b波段的经过恢复后的DN值，CR_(b)是该 像元在第b波段的CR值。

优选地，步骤A中通过eCognition Developer 8.0软件对所述变换后所 得的图像进行分割，在使用eCognition Developer 8.0软件对图像进行图像 分割时，形状因子Ssh设值为0.1，紧致度因子Scm设值为0.5，根据图像 中地物的异质或匀质程度设定尺度因子Ssc的值。

优选地，步骤A中的自适应阈值算法满足如下条件，

σ_B²＝ω₀(μ₀-μ_T)²+ω₁(μ₁-μ_T)²

其中，μ₀和μ₁分别表示由根据特定阈值进行分割后得到阴影区和非阴 影区像元亮度的平均值，μ_T表示原始图像整体的亮度平均值，ω₀和ω₁则分 别表示阴影和非阴影像元个数占总像元数的比例，σ_B²表示阴影区和非阴影区 的差异性，当σ_B²最大时，所述特定阈值就是最佳分割阈值。

优选地，步骤B中缓冲区的宽度上限设为10，如果在最大的缓冲区内相 似像元数量仍然无法满足大于等于20的要求，则将相似像元实际数量作为相 似像元数。

本发明所提供的一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢 复方法，用于进行遥感图像阴影区的光谱信息恢复，该方法充分结合了阴影 像元本身的弱光谱信息及其邻近相似像元的光谱信息，摆脱了传统地形纠正 算法对DEM数据的依赖，当DEM数据的分辨率和精度不能满足地形纠正算法 得要求时，该算法可以提供一种行之有效的阴影处理方法。

附图说明

以下附图仅旨在于对本发明做示意性说明和解释，并不限定本发明的范 围。其中，

图1为根据本发明的一个具体实施例的一种Landsat TM/ETM+图像中山 体阴影区的光谱信息恢复方法的流程示意图；

图2为图1所示的Landsat TM/ETM+图像中山体阴影区的光谱信息恢复 方法的非阴影相似像元搜寻示意图；

图3为一幅Landsat TM影像图片；

图4a为对图3所示的Landsat TM影像图片使用根据本发明的一个具体 实施例的一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法进行 阴影恢复后的图片；

图4b为对图3所示的Landsat TM影像图片使用ATCOR模型的地形纠正 进行阴影恢复后的图片；

图4c为图4a的局部放大图；

图4d为图4b的局部放大图。

具体实施方式

为了对本发明的技术特征、目的和效果有更加清楚的理解，现说明本发 明的具体实施方式。

图1为根据本发明的一个具体实施例的一种Landsat TM/ETM+图像中山 体阴影区的光谱信息恢复方法的流程示意图；参照图1所示，下面详细说明 根据本发明的一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法 的原理。所述方法包括如下步骤：

步骤A对阴影区进行识别和提取；

步骤B非阴影相似像元搜寻；

步骤C建立在阴影像元与其非阴影相似像元间的光谱相似性基础上的阴 影区光谱信息恢复。

考虑到大气纠正可能会导致阴影像元弱光谱信息的丢失，本发明优选用 于对大气纠正前的原始图像，即DN值图像进行处理。

步骤A：对阴影区进行识别和提取

对Landsat TM/ETM+图像中阴影区的准确识别和提取是对其进行光谱信 息恢复的先决条件，研究者们也提出了许多行之有效的阴影识别方法，例如：

Salvador，E.，Cavallaro，A.，& Ebrahimi，T.(2001).Shadow  identification and classification using invariant color models.In  IEEE International Conference on Acoustics，Speech and Signal  Processing，3，1545-1548.)；

Giles，P.T.(2001).Remote sensing and cast shadows in mountainous  terrain.Photogrammetric Engineering and Remote Sensing，67，833-839；

Tsai，V.J.D.(2006).A comparative study on shadow compensation  of color aerial images in invariant color models.IEEE Transactions  on Geoscience and Remote Sensing，44，1661-1671.；

Chen，Y.，Wen，D.，Jing，L.，& Shi，P.(2007).Shadow information  recovery in urban areas from very high resolution satellite imagery. International Journal of Remote Sensing，28，3249-3254；

然而大部分方法比较繁琐复杂，或者仍然需要DEM数据的辅助。

本发明利用Landsat TM/ETM+图像中阴影区和非阴影区在表观亮度上的 明显差异来对二者进行识别区分，从而对山体阴影区进行提取。

首先，对需要进行处理的图像进行缨帽变换，也即是TC变换，可参见 Crist，E.P.，& Cicone，R.C.(1984).Application of the Tasseled Cap  Concept to Simulated Thematic Mapper Data.Photogrammetric  Engineering and Remote Sensing，50，343-352.，并将变换后所得的图像 的各个像元的亮度波段数值作为每个像元的亮度值，当所述图像上水体面积 较大时，考虑到水体在图像上的表观亮度较低，容易与阴影区混淆，因此可 以先借助所述图像所在的地表区域已有的水体图像进行水体掩膜。

因为Landsat TM/ETM+图像中山地地区的山体阴影往往具有一定的形状 和大小，而并非呈离散的点状分布，因此在对Landsat TM/ETM+图像中阴影 区进行识别和提取时，以对象区域尺度进行识别和提取比以像元尺度进行识 别和提取更加合理。

对所述图像进行缨帽变换后，根据变换计算得到的像元的亮度值，对所 述变换后所得的图像进行分割，将所述变换后所得的图像划分为多个对象区 域，其中每一个对象区域内部的像元具有相似的空间性质和光谱性质，

在对所述图像进行分割时，可以采用eCognition Developer 8.0软件进 行，在使用eCognition Developer 8.0软件对图像进行图像分割时，形状因 子Ssh设值为0.1，这样可以给予光谱特征更高的权重，同时可以将紧致度 因子Scm设值为0.5，这样可以对分割出的所述对象区域的边界给予相同程 度的光滑度和粗糙度。

尺度因子Ssc控制着分割出的所述对象区域的尺度大小水平，由于最佳 的尺度水平会随着具体图像和研究地区的不同而具有很大的变化范围，因此 通常可以根据具体图像中地物的异质或匀质程度来经验性地设定尺度因子 Ssc的值，即图像异质性越高，该参数值越小。还可以采用基于分割质量评 判准则提出的尺度因子最优化算法来进行该参数的确定，例如Chabrier，S.， Emile，B.，Rosenberger，C.，Laurent，H.(2006).Unsupervised  performance evaluation of image segmentation.EURASIP Journal on  Applied Signal Processing，1-12所记录的算法。也就是选择一个最合适 的尺度因子Ssc，使得分割后图像整体的对象内匀质性最高，而整体的对象 间异质性最高。所述图像整体的对象内匀质性和对象间异质性，可以分别用 全局加权方差wVar和全局空间自相关指数MI计算得出，所述全局加权方差 wVar和全局空间自相关指数MI的原理可参见Espindola，G.M.，Camara，G.， Reis，I.A.，Bins，L.S.，& Monteiro，A.M.(2006).Parameter selection  for region-growing image segmentation algorithms using spatial  autocorrelation.International Journal of Remote Sensing，27， 3035-3040.，所述全局加权方差wVar公式如下：

$\mathrm{wVar}=\frac{\underset{i=1}{\overset{n}{\Sigma }}{a}_i*v_i}{\underset{i=1}{\overset{n}{\Sigma }}{a}_i}$

其中n表示分割图像中总的对象区域个数，v_i表示第i个对象区域内部 像元值的方差，a_i表示第i个对象区域内的像元个数；所述全局空间自相关 指数MI公式如下：

$\mathrm{MI}=\frac{n\underset{i=1}{\overset{n}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}{w}_{\mathrm{ij}}(y_i-\bar{y})(y_j-\bar{y})}{\underset{i=1}{\overset{n}{\Sigma }}{(y_i-\bar{y})}^2\left(\underset{i\ne j}{\Sigma }\Sigma w_{\mathrm{ij}}\right)}$

其中n表示分割图像中总的对象区域个数，y_i表示第i个对象区域内部 像元的平均值，表示整个图像像元的总体平均值，w_ij用于衡量第i个对象 区域和第j个对象区域间的空间临近性，当二者为相邻对象区域时，w_ij＝1， 否则为w_ij＝0。

将所述全局加权方差和所述全局空间自相关指数结合起来，便可得到一 个衡量图像分割整体效果的综合性指数，即全局分数GS，所述全局分数GS 公式如下：

GS＝V_norm+MI_norm

其中Vnorm和MInorm分别表示全局加权方差和全局空间自相关指数的极值标 准化值。

当GS取得最小值的时候，此时的尺度因子Ssc即为最优值。

在对所述图像分割后，可以计算出每个所述对象区域内部像元的亮度平 均值，将所述亮度平均值作为该对象区域的亮度值，然后可以利用自适应阈 值算法得到阈值，根据所述阈值对阴影区进行提取。

所述自适应阈值算法可以采用Otsu自适应阈值算法，本领域技术人员参 照文献Otsu，N.(1979).A threshold selection method from grey-level histograms.IEEE Transactions on Systems，Man，and Cybernetics，9， 62-66即可通过该算法获得一个阈值，使得一个整体在使用该阈值进行分割 的情况下，得到的两个子数据集具有最高的差异性。

具体来说，在本发明中，参照上述Otsu自适应阈值算法，所述自适应阈 值满足如下条件：

σ_B²＝ω₀(μ₀-μ_T)²+ω₁(μ₁-μ_T)²

其中，μ₀和μ₁分别表示由根据阈值进行分割后得到阴影区和非阴影区 像元亮度的平均值，μ_T表示原始图像整体的亮度平均值，ω₀和ω₁则分别表 示阴影和非阴影像元个数占总像元数的比例，σ_B²表示阴影区和非阴影区的 差异性，当σ_B²最大时，所述阈值就是最佳分割阈值。

步骤B：非阴影相似像元搜寻

虽然直射光的缺乏导致对Landsat TM/ETM+图像中阴影像元的表观亮度 一般较低，但由于天空散射光的存在，所述阴影像元仍然保留了较弱的光谱 信息，而这种弱信息所提供的光谱形状与同类别地物的非阴影像元的光谱形 状具有较高的相似性。可以将对Landsat TM/ETM+图像中像元的光谱曲线分 离成两部分，一部分表征其DN值或反射率的绝对高低水平，另一部分则仅体 现其曲线形状特点，那么后者的相似性便可用于衡量对Landsat TM/ETM+图 像中阴影像元与邻近的非阴影像元间地物类别的相似程度。

可以使用连续统去除方法进行相似像元比对，连续统去除方法原是一种 用于光谱吸收特征分析的亮度标准化技术，相关文献为

Clark，R.N.，& Roush，T.L.(1984).Reflectance spectroscopy： Quantitative analysis techniques for remote sensing applications. Journal of Geophysical Research，89，6329-6340.；

Youngentob，K.N.，Roberts，D.A.，Held，A.A.，Dennison，P.E.， Jia，X.，& Lindenmayer，D.B.(2011).Mapping two Eucalyptus subgenera using multiple endmember spectral mixture analysis and  continuum-removed imaging spectrometry data.Remote Sensing of  Environment，115，1115-11.

该方法将一条光谱曲线分离为两部分：连续统信息CI表征绝对亮度水 平；连续统去除信息CR表征光谱的形状特征。具体的，连续统是一条连接原 始光谱各个局部峰值点的平滑曲线，亦即上包络线；在每一个波段或波长上， 求取原始光谱值OI与连续统值的商，则形成一条取值范围在0到1之间的新 的曲线，即连续统去除曲线。对应计算公式为

${\mathrm{CR}}_{\lambda }=\frac{{\mathrm{OI}}_{\lambda }}{{\mathrm{CI}}_{\lambda }}$

由于相同地物类别的阴影像元和非阴影像元具有相同或相似的光谱形状 特征，因此它们的CR曲线之间则应具有相对较高的相似性，这便是非阴影相 似像元搜寻的理论基础。

具体在本发明中，根据步骤A中图像分割的结果，进行相似像元的寻找， 对每一个阴影像元，在其邻近空间范围内搜索出一定数量的非阴影相似像元， 以便后继用于对阴影像元光谱信息的恢复。

在所述变换后所得的图像上建立平面二维x-y坐标系，这样就能够用二 维坐标来确定所述变换后所得的图像上每个像元的位置。

图2为图1所示的Landsat TM/ETM+图像中山体阴影区的光谱信息恢复 方法的非阴影相似像元搜寻示意图；参照图2所示，在图像分割后，阴影区 以对象区域的形式被提取出来，对于每一个阴影对象区域，首先在其周围建 立至少两个像元宽度的缓冲区，然后在此缓冲区范围内，对应所述对象区域 中每个目标阴影像元选出与所述目标阴影像元光谱相似度最高的至少N个非 阴影像元，作为所述目标阴影像元的非阴影相似像元，此处，相似度由像元 间光谱曲线的标准偏差RMSD来衡量：

$\mathrm{RMS}{D}_i=\sqrt{\frac{\underset{b=1}{\overset{n}{\Sigma }}{(\mathrm{CR}(x_i,y_i,b)-\mathrm{CR}(x,y,b))}^2}{n}}$

其中CR(x_i，，y_i，b)是位于(x_i，，y_i)处的像元在第b波段的连续统去除信 息CR值，CR(x，y，b)则是目标阴影像元在第b波段的连续统去除信息CR值， n为总的波段数。RMSD越小，则像元间的相似度越高，反之则越低。

被选择作为所述目标阴影像元的所述非阴影相似像元满足下列条件

${\mathrm{RMSD}}_i\le [\underset{b=1}{\overset{n}{\Sigma }}\sigma \left(b\right)\times 2/m]/n$

其中σ(b)表示整个图像在第b波段的连续统去除信息CR值的标准差，m 表示图像中的地物类别数，其值可以根据目视判别和经验，也可以是根据已 有的土地覆盖类型图进行预先人为设定。

考虑到混合像元和随机误差的存在，非阴影相似像元的数量N必须要达 到一个最低要求，以保证据此进行的阴影像元光谱信息恢复具有统计稳定性， 相关文献可参见Chen，J.，Zhu，X.，Vogelmann，J.E.，Gao，F.，& Jin， S.(2011).A simple and effective method for filling gaps in Landsat  ETM+SLC-off images.Remote Sensing of Environment，115，1053-1064。

如果在两个像元宽度的缓冲区范围内，满足条件的所述非阴影相似像元 数量个数N无法满足大于等于20的要求，则以两个像元为步长扩充缓冲区的 宽度，直到此要求被满足为止。

为了节约运算资源，也可以缓冲区的宽度上限设为10，如果在最大的缓 冲区内相似像元数量仍然无法满足大于等于20的要求，则将相似像元实际数 量作为相似像元数。

步骤C：阴影区光谱信息恢复

考虑到阴影像元与其非阴影相似像元的光谱曲线具有相似的形状特征和 不同的亮度水平，因此可以用所述非阴影相似像元提供的平均亮度信息替换 所述阴影像元的亮度信息，从而实现对所述阴影像元的光谱信息进行反演， 所述阴影像元的亮度信息来自所述非阴影相似像元连续统去除后得到的连续 统信息CI，经此处理后，所述阴影像元的光谱特征得到了恢复，其表观亮度 也被提高到与其周围非阴影像元相近的水平，从而使得所述阴影区的形状特 征得以表现，

对于每个非阴影相似像元，根据该像元与对应的所述阴影像元的空间距 离以及光谱相似度，将二者分别标准化后相结合，即可得出表示所述非阴影 相似像元对用于进行光谱恢复的平均亮度值的贡献程度的权重值W，例如， 对于第j个非阴影相似像元，其对平均亮度值的贡献程度的权重值W_j可由下 式得出$W_j=\frac{1/({D_j}^{*}\times {{\mathrm{RMSD}}_j}^{*})}{\underset{j=1}{\overset{N}{\Sigma }}1/({D_j}^{*}\times {{\mathrm{RMSD}}_j}^{*})}$

上述公式中，和分别是该第j个非阴影相似像元与对应的目标 所述阴影像元的空间距离D_j以及它们之间的光谱相似度RMSD_j的标准化之后 所得的数值；

第j个非阴影相似像元与目标所述阴影像元的空间距离D_j的计算公式为

$D_j=\sqrt{{(x_j-x)}^2+{(y_j-y)}^2}$

的计算公式为：$D_j^{*}=\frac{D_j-D_{\min }}{D_{\max }-D_{\min }}+1$

的计算公式为：${\mathrm{RMSD}}_j^{*}=\frac{{\mathrm{RMSD}}_j-{\mathrm{RMSD}}_{\min }}{{\mathrm{RMSD}}_{\max }-{\mathrm{RMSD}}_{\min }}+1$

上述公式中下脚标max和min分别表示非阴影相似像元中该量的最大值 和最小值

根据每个所述非阴影相似像元的所述权重值W，即可计算出所述阴影像 元的对应的非阴影相似像元亮度的加权平均值CI_wavg(b)，所述加权平均值 CI_wavg(b)的计算公式如下

${\mathrm{CI}}_{\mathrm{wavg}\left(b\right)}=\underset{j=1}{\overset{N}{\Sigma }}{W}_j\times {\mathrm{CI}}_{j\left(b\right)}$

其中权重值W_j表示第j个非阴影相似像元对平均亮度值的贡献程度， CIX表示第j个非阴影相似像元在第b波段的连续统信息，也即是第j个 非阴影相似像元的绝对亮度水平

根据所述加权平均值CI_wavg(b)即可得到所述阴影像元在第b波段的经过恢 复后的DN值，计算公式如下，

OI_re(b)＝CI_wavg(b)×CR_(b)

其中OI_re(b)表示目标阴影像元在第b波段的经过恢复后的DN值，CR_(b)是 该像元在第b波段的CR值。

如上所述，得到目标阴影像元在第b波段的经过恢复后的DN值之后，即 可使得阴影像元的亮度值被提高到与其周围非阴影像元相近的水平，也就能 够恢复阴影像元部位的地物光谱特征。

为了更好的说明本发明的技术效果，针对一幅原始Landsat TM影像，分 别利用本发明所提供的一种Landsat TM/ETM+图像中山体阴影区的光谱信息 恢复方法以及传统的ATCOR模型进行阴影恢复，然后对阴影恢复后的结果进 行比较。ATCOR模型的详细描述参见Richter，R.(1998).Correction of  satellite imagery over mountainous terrain.Applied Optics，37， 4004-4015.

图3为一幅Landsat TM影像图片，该图片为2009年6月30日获取的 Landsat TM影像，影像区域位于东经109°1′24″，北纬33°57′22″的 中国陕西省中部，覆盖面积约1358平方公里。在该区域的中部和南部地区， 森林是主要的覆盖类型；北部地区则主要分布着农田、城市和裸地。除此之 外，辅助数据包括：相应区域的水体图像用于水体掩膜，SRTM提供的90米 分辨率的DEM数据用于进行传统地形纠正，以及相应区域的QuickBird高分 辨率影像用于地物类型的解译和分类精度的检验。由于DEM数据的空间分辨 率与TM影像并不匹配，因此需要先对DEM数据进行重采样，另外，所有的辅 助数据也需要与TM影像进行几何配准。

图4a为对图3所示的Landsat TM影像图片使用根据本发明的一个具体 实施例的一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法进行 阴影恢复后的图片；图4b为对图3所示的Landsat TM影像图片使用ATCOR模 型的地形纠正进行阴影恢复后的图片；图4c为图4a的局部放大图；图4d为 图4b的局部放大图。参加图3-4所示，通过利用根据本发明的一种Landsat  TM/ETM+图像中山体阴影区的光谱信息恢复方法的阴影纠正，图3所示的原始 影像中部分黑暗的阴影区的亮度被提升到了与周围非阴影区相近的水平，而 利用ATCOR模型的地形纠正进行阴影恢复的效果却较差，经ATCOR处理后， 山体阳坡与阴坡之间的亮度差异的确得到一定的纠正，但图像却变得很破碎， 同时仍然存在许多暗线和暗区域。

由上述图片实例可看出，根据本发明的一种Landsat TM/ETM+图像中山 体阴影区的光谱信息恢复方法可以很好的恢复阴影像元部位的地物光谱特征

本发明所提供的一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢 复方法，用于进行遥感图像阴影区的光谱信息恢复，该方法充分结合了阴影 像元本身的弱光谱信息及其邻近相似像元的光谱信息，摆脱了传统地形纠正 算法对DEM数据的依赖，当DEM数据的分辨率和精度不能满足地形纠正算法 得要求时，该算法可以提供一种行之有效的阴影处理方法。

本领域技术人员应当理解，虽然本发明是按照多个实施例的方式进行描 述的，但是并非每个实施例仅包含一个独立的技术方案。说明书中如此叙述 仅仅是为了清楚起见，本领域技术人员应当将说明书作为一个整体加以理解， 并将各实施例中所涉及的技术方案看作是可以相互组合成不同实施例的方式 来理解本发明的保护范围。

以上所述仅为本发明示意性的具体实施方式，并非用以限定本发明的范 围。任何本领域的技术人员，在不脱离本发明的构思和原则的前提下所作的 等同变化、修改与结合，均应属于本发明保护的范围。