山地区域的遥感影像会因地形引起目标辐射特性差异，产生同物异谱或同谱异物现象，影响地物解译、分类与反演的结果^[1]。地形辐射校正旨在削弱甚至消除地形因素对辐射值产生的影响。

阴影是地形效应的最主要形式，其对地表辐射值的减弱降低了影像数据的利用效率，所以在地形辐射校正前先进行阴影检测与修复，有利于遥感影像的分析与应用。阴影检测方法主要有模型法与特征法两种方法^[2]，单一的方法总是存在局限性，结合高空间、高光谱分辨率数据，以及高精度的DEM数据，可以同时利用两种方法的优势，提高阴影检测的精度^[3]。

对于检测出的阴影区进行修复时，修复方法主要有阴影增强^[4]、基于整体变分阴影补偿^[5]和纹理匹配阴影补偿等^[6]，也有重建光源模型^[7]，利用阴影的恒常特性估计标准光照条件下阴影区域地物应有的信息的补偿方法，以及用图像增强方法，如线性拉伸修正、分形插值、同态滤波^[8-10]等，都可以恢复阴影区域的信息。由于山地阴影的特殊性，阴影与非阴影之间不能找到合适的过渡区，难以直接利用阴影修复方法对山区阴影的像素进行拉伸与修复，并且阴影修复没有考虑向阳面过亮的问题，不能很好地达到地形效应的去除效果。

传统地形辐射校正结合了DEM的信息，可以同时处理向阳面过亮以及背阴面过暗的问题，常用的模型有C、SCS、cosin、Minnaert、VECA、三因子改正模型等^[11]，还有针对C校正的改进模型^[12-13]、Minnaert的改进模型^[14-15]，基于变分Retinex的校正法^[16]等。但是地形辐射校正物理模型只依据DEM与太阳之间的关系，逐像元地进行校正，山体间的相互关系考虑不充分，出现过校正和欠校正的情况。

本文提出一种基于DEM的山体阴影检测与地形辐射校正算法，充分利用了影像数据的多波段信息、DEM以及山体间的相互关系，将阴影按照定义进行更为精确的检测分类，并结合阴影修复与地形辐射校正模型，依据阴影的产生原因，对不同的地形区域进行不同方法的校正，对阴影较深的区域进行更好地修复，从而对传统物理模型的地形辐射校正效果有所改进。

1 山体阴影检测

阴影从构成原理上分为本影、投影以及半影区，这3种阴影在阴影检测中都被分为阴影区^[17]。由于产生阴影的原因不同(如图 1所示)，在恢复不同阴影处信息的时候，就需要不同的校正方法。同时，有的山区背阴面在产生原理上不是阴影，但是在遥感影像中往往亮度较暗，在一般阴影检测中被判为阴影，这些都会影响地形辐射校正的精度。

图 1中θ_S为太阳天顶角；S₁—S₄为不同坡面的坡度角；i为太阳入射角。

本文提出的阴影检测方法融合了模型法与特征法，将阴影检测分为3个独立的步骤，对山区阴影进行分类处理，然后综合分析。

首先用多波段信息进行基于色彩空间的阴影检测。将RGB波段转到HSI空间，引入蓝通道归一化信息，并用近红外波段的数据去除云和水^[18-19]，最终用阈值法得到阴影区和非阴影区。利用特征P，如式(1)所示，增强低亮度和高色调的对比，使阈值更加明显，并改善阴影检测^[2]

式中，H与I分别是HSI色彩空间的色调与亮度。根据式(1)得到最终的阴影检测结果

式中，set₁为第1步阴影检测结果；set_1S为阴影像素集合；set_1L为非阴影像素集合。

然后利用DEM和太阳角度信息计算得到每个像素位置处的太阳入射角信息^[20]，如式(3)所示

式中，i为太阳入射角；θ_S为太阳天顶角；S为坡度角；ϕ_S为太阳方位角；A为坡向。用阈值法区分出水平面、背阴面以及向阳面。根据式(3)可得

式中，set₂为第2步阴影检测结果；set_2S为背阴面像素集合；set_2L为非背阴面像素集合。

接着进行本影以及投影检测，依据地形的相关关系与太阳高度角进行计算^[11]。本影的产生条件是在太阳入射方向上的坡度角θ_sn大于或等于太阳高度角θ_S。其中θ_sn如式(5)所示

式中，θ_n为坡度；ϕ_s为太阳方位角；ϕ_n为坡向；太阳方位角ϕ_s与坡向ϕ_n以正北为起始方向，顺时针为正。

若θ_sn≥θ_S，则为本影，反之不产生本影。投影的检测为在太阳入射方向按间隔取采样点，分别用高程计算地形遮蔽角θ_j。若θ_j≥θ_S，则视为遮蔽，并计算在该条射线上遮蔽点所占的比例，就可以得出该点受遮蔽的程度，从而进行相对应的校正，其中θ_j的计算如式(6)所示

式中，H_j为各采样点高程；H_O为地面单元O点的高程；θ_j是地面单元O点至j点间的仰角；ΔL为采样点间的步长。最终得到

式中，set₃为第3步阴影检测结果；set_3S为投影与本影像素集合；set_3L为非投影影像素集合。

2 阴影补偿

由于山区阴影产生原因的复杂性，单一的地形辐射校正或者阴影补偿方法不能对各种情况做到合理的处理。本文综合阴影修复与地形辐射校正模型的长处，并结合对cos i的调整，对山区阴影进行补偿与校正。

阴影检测的结果中，对阴影进行了细致的定义与划分，所以在阴影修复中，对不同的阴影进行区别处理。3种结果综合，就组合成了8种不同的地形与阴影定义，整理后见表 1。

表 1中详细展示了阴影检测结果中各种情况的出现原因，并针对每种情况作出对应的校正，将3种阴影检测结果分别赋值于红、绿、蓝波段，得到综合的结果图，其中红色为set_1S，绿色为set_2S，蓝色区为set_3S，非阴影区及向阳面等赋值为0，故为黑色。

在校正方法中主要有传统校正法与阴影补偿法，并且依照经验模型对太阳入射角余弦值进行调整。传统地形辐射校正法，是逐像素的按照几何物理模型校正，本文分别采用的C校正^[23]、SCS+C^[24]校正与VECA^[25]校正，其模型分别如式(8)—(10)所示

式中，L_H为水平面的辐亮度值；L_T为倾斜表面的辐射值；θ_S为太阳天顶角；S为坡度角；i为太阳入射角；C_(λ)为分波段计算得到的统计值；L_a为校正前像素辐亮度平均值；m、b分别为辐射值与cos i线性回归方程的斜率与截距。

阴影补偿中，首先采用小波变换的方法，将影像分为高频与低频，在对低频部分用一种利用局部均值与方差进行补偿的映射模型进行增强处理^[22]，如式(11)所示

式中，LC(i, j)与LC′(i, j)分别为原始与校正后的影像低频系数；m_L、m_S分别为该区域理想情况下以及阴影情况下低频系数的均值；σ_L与σ_S分别为理想情况下与阴影情况下低频系数的均方差。

在这里，阴影补偿的目的是将投影区的亮度调节成为不受投影影响的情况，这就需要在该区域附近寻找地形相似，且地物类型相近的非投影区，用这个作为理想情况下的亮度，对投影区进行亮度的拉伸，所以也叫作投影补偿。投影补偿后的影像亮度，仍旧是倾斜地表的亮度，只是去除了遮挡关系的影响，因此投影补偿后的区域还是要结合地形辐射校正模型，才能得到水平地表的亮度值。

在投影区域附近寻找不受投影影响的区域时，首先用HSV变换得到亮度图，用于将少数高亮的地面点视为粗差去除，然后在该像素周围邻域内计算cos i、坡度以及高程的距离，如式(12)所示

式中，i_j与i_o分别为非投影与投影区太阳入射角；hight_j与hight_o分别为非投影与投影区高度；slope_j与slope_o分别为非投影与投影区坡度角。对归一化后的3个量进行计算，并用dist_t选出地形相近的像素集，然后计算该像素集中的光谱角θ与NDVI，并计算欧氏距离dist_s，如式(13)所示，用阈值法找到与原像素地物类型也相近的像素集

式中，θ_j与θ_o分别为像素集中待匹配像素与原像素的光谱角；NDVI_j与NDVI_o分别为像素集中待匹配像素与原像素的归一化植被指数。光谱角计算^[26]如式(14)所示

式中，θ为光谱角；X与Y分别是两像素处不同波段辐射亮度值组成的一维向量。

对于选定的地形与地物都相近的像素集，在求均值与方差的时候，用式(15)作为权值，以保证与原像素点地形和地物类型越相似的，对原像素的影响越大，避免出现不同区域植被不同，以及大气情况复杂等诸多因素引入不必要的误差，从而带来的精度降低。最后对投影区校正后的影像进行地形辐射校正，得到结果

式中，dist_t表示地形的相似度；dist_s表示光谱的相似度。

对于影像检测结果为阴影，但是却是向阳面的部分，采用调整太阳入射角余弦值cos i的方法进行校正^[21]。这些区域在传统地形辐射校正模型中不能得到很好的校正，是因为该区域受环境光以及邻近像元区域的影响，计算所得cos i与真实情况不相符，使得辐射值与cos i并不满足线性关系，所以用式(16)对cos i进行调整，得到cos i′作为改正后的太阳入射角余弦值，代入校正模型进行校正

式中，cos i与cos i′分别为原始于调整后的太阳入射余弦值角；L_a为校正前像素辐亮度平均值；m、b分别为辐射值与cos i线性回归方程的斜率与截距。

3 试验数据与结果分析

为对本文方法进行验证，选用Landsat 8湖北双峰山区域的影像进行试验，影像的空间分辨率为30 m。DEM为5 m分辨率的数据，由资源三号影像密集匹配得到。为了保证各种情况的阴影都存在，选取的测试区域高程差为800 m，平均高程为160 m，并选用太阳高度角较低的影像进行试验。本文分别对两组数据采用3种传统地形辐射校正模型进行试验，并将本文方法结合3种方法，得到对比的结果，试验结果如图 2-7所示。

DEM数据与Landsat 8影像的空间对应精度会影响试验的结果^[27]，本试验中Landsat 8为正射影像，与DEM在空间上可以对应，不在此讨论空间对应精度对试验结果的影响。

3.1 阴影检测结果

第1组影像的太阳方位角为151.501 9°，太阳高度角为48.177 3°，第2组影像的太阳方位角为150.117 2°，太阳高度角为50.269 6°。

图 2(a)和图 3(a)为同一地区不同时间的Landsat 8影像。图 2(d)和图 3(d)是3种阴影检测结果的综合分析，其中：

红色部分是利用多波段数据检测的阴影结果，就是原影像中较暗，所以信息量减少的区域，即图 2(b)和图 3(b)的检测结果。

绿色部分是利用DEM与太阳角度信息计算所得cos i值得到的背阴面，cos i结果如图 2(c)和图 3(c)所示，用阈值法就可以得到cos i大的山体背阴面。

蓝色区域是利用山体间的相关关系所得落影检测结果，黑色区域为无影区。

详细的阴影分类分析见表 1，这样就实现了从原理上对不同的阴影进行区分，使得校正更加精确有效。图 2和图 3的3种阴影检测结果对比中，可以看出单一的阴影检测结果都不能完整地包含所有可能产生阴影区域，对阴影并不能做出很好的判定。由于每种结果都将影像分为阴影与非阴影两部分，将3种结果综合分析便得到8种阴影的分类。表 1对8种结果作了详细说明，可以对地形和阴影有更全面的把握。

3.2 地形辐射校正结果

根据表 1对不同地形及阴影的分类，得到了多波段阴影检测结果、背阴面及投影区的综合分析结果，结合地形辐射校正模型C、SCS+C及VECA，分别对两张影像用不同的方法进行组合校正，得到的地形辐射校正结果与原校正模型的校正结果进行对比，如图 4—7所示。

该区域地表覆盖类型主要是植被，同时也有水、裸地及其他地物类型。可以看出本文方法进行的阴影校正可以更加有效提高阴影区的亮度，有效地降低了波段辐射值与地形的关系，并且对无须校正的水与裸地区域，以及水平的植被区，都能保持原辐亮度值，没有改变地物类型。

对得到的影像结果进行统计特征分析，得到各方法校正前后的统计值，均值体现整张影像的亮度，均值越大表明整体亮度越亮，方差体现色彩的均匀性，在地形辐射校正中就说明背阴面与向阳面的亮度差异越小，地形辐射校正的效果也就越好。表 2中对整幅影像的全部像素进行了统计值的计算，并对比校正前后的均值与方差值。可以看出本文方法对影像均值有所提升，方差有所减小，以波段3为例，校正前的均值为35.882 5，C校正后均值为36.289 8，结合本文方法后，均值为36.562 2，同样，加入本文方法后，原方差7.819 8也减小为7.454 78。这是因为本文方法对较暗的投影面与自影面，都在保证原地物类型的前提下，进行了有效的恢复。

表 3对比了校正前后斜率值。斜率体现了各个波段的辐亮度值与地形之间的相关关系，是衡量地形辐射校正结果的最重要参数，斜率越小，则辐亮度值于地形的线性相关越小，即影像受地形影响越小。如表 3所示，以波段4为例，校正前的斜率为13.395 9，经C校正后减小为2.858 45，结合本文方法后，则减小到1.885 01，可以看出本文方法在各个波段都能对斜率有所减小，对传统的地形辐射校正模型有改进的效果。

为了评价地形辐射校正的效果，进行进一步统计定量分析，将影像的辐射亮度值与cos i结合，画出散点图并进行线性回归拟合。如图 8所示，图中纵坐标为影像三波段的辐射亮度值，横坐标为太阳入射角余弦值cos i，r为相关系数，体现二者的相关性，mean为均值，SD为方差，图中点颜色越浅表示越密集。

从计算结果可以看出校正前辐射亮度值与cos i具有较强的相关关系，斜率为14.783 1，相关系数为0.241 7，经C校正后斜率为3.362 04，相关系数为0.036 6，已经有了一定改进，再结合本方法后，斜率减小为2.190 62，相关系数减小为0.013，已经非常接近于0，说明辐射亮度值与cos i的相关性已经基本消除，再结合均值与方差，可以看出本文方法较未改进前有明显的改善。

4 结论

传统地形辐射校正主要是像素级的校正，针对每个像素的地理位置以及与太阳的入射角度关系进行辐射值的校正，虽然这个过程是依据物理模型的，可以得到较理想的结果，但还是存在如过校正或欠校正等问题。本文在阴影检测过程中引入投影的计算，考虑了山体间的相互关系，并且用多波段阴影检测结果进行辅助，即对地形类别和阴影进行了更加精细的分类，考虑到每种阴影产生的原因以及校正方法，并针对每种地形和阴影分别进行校正，从而得到更加理想的结果。这个过程可以理解为对象级的地形辐射校正，从原理上改进了校正的过程与结果。

同时，3种阴影检测结果之间可以相互约束，并且增强了阈值法阴影检测中对阈值选择正确性的宽容度，从而增强了算法的稳定性。

本文试验还有一定的局限性。由于缺少临近像素之间辐亮度的相互约束，会出现个别像素亮度变化不是很平缓，可加入超像素分割法以解决，即用光谱与纹理的统计值在超像素间进行阴影区和非阴影区的匹配与阴影去除，增强了相邻同质像素间的相互约束，从而改善这个问题。