高空间分辨率遥感影像(高分影像)以其精细的地物分辨能力，在资源环境遥感中发挥着重要作用。城市作为人类聚集区，其土地资源、生态环境与人类生活息息相关。利用高分影像进行城市资源环境监测是城市资源环境规划、保护与管理的重要手段和数据源(Myint等，2011;陈基伟和程之牧，2004)。然而，由于城区建筑阴影的广泛存在，城市高分影像中常常存在大量阴影，这严重制约了高分遥感数据的高效利用(Dare，2005)。

城市阴影在遥感影像上表现为低值区，并掩盖了阴影区地物本身辐射特性。当利用高分影像进行城市土地利用分类、资源调查时，这些低值区将极大地降低结果精度，因而降低数据使用效率(Zhan等，2005)。同时，城市高分影像获取成本通常较高，较低数据利用效率势必增加数据使用成本(陈锬，2014)。此外，由于城市阴影是固定的建筑等地物对太阳直接辐射遮挡产生的，而太阳方位角和高度角的变化规律也是固定的，这可能使得一些阴影区域在不同时期的遥感影像上一直存在，从而形成城市遥感监测的永久盲区(王卫国，2008)。因而，开展城市阴影校正(去除)并恢复阴影区有效地表辐射信息是城市遥感数据预处理的重要内容(Zhou等，2009)。

目前，国内外已开展了大量有关阴影校正算法方面的研究(Chen等，2007;Dare，2005;Jin等，2013;Richter和Müller，2005;Zhou等，2014)，其中经典方法有直方图匹配法(Histogram Matching)，伽玛(γ)校正(Gamma Correction)以及均值方差变换法MVT(Mean and Variance Transformation)(Shahtahmassebi等，2013)等。直方图匹配法是利用阴影区与临近非阴影区存在的相同地物进行直方图匹配，而后将该匹配参数应用于整个阴影区，从而实现阴影校正(Dare，2005;Sarab and i等，2004)。伽玛校正采用单一的γ因子用于整个阴影区，忽略了不同的阴影区背景地物的差异性特征(Nakajima等，2002)。均值方差变换法(MVT)则假设阴影区和非阴影区存在线性关系，通过提取各自区域的统计变量(均值与方差)构建该线性关系，从而实现阴影区校正(Chen等，2007;Zhan等，2005)。相关研究表明，MVT法以其较高的校正精度和操作简单方便的特点，相对于其他阴影校正方法具有明显的优越性，被广泛使用(Chen等，2007;Yamazaki等，2009;Liu和Yamazaki，2012)。

综观这些阴影校正方法，其核心都利用了阴影区内被记录的微弱辐射信号，并在此基础上利用阴影区和非阴影区之间的某种特殊统计特征关系实现阴影校正的(Shahtahmassebi等，2013)。事实上，阴影区与非阴影区除了存在这些统计特征关系之外，还广泛存在另一种物理特征关系：即在朗伯体表面假设下，同一地物在阴影区及其临近非阴影区具有相同的反射率(Pedrotti和Pedrotti，1993)。目前，这种关系在已有阴影校正方法体系中尚未被充分发掘(wen等，2015)。本研究将以此关系为基本出发点，发展一种基于辐射传输过程的城市高分遥感影像阴影区辐射校正方法。

由于大气的散射作用，太阳辐射在经过大气层达到地表的过程中分为直接辐射(I_d)和散射辐射两部分(图 1)。直接辐射在受到地表附近不透明物体(如建筑)阻挡时，形成阴影区。散射辐射可进一步分为向上和向下传输部分，其中向上传输部分沿传感器观测方向的叫做路径程辐射(I_p)，向下传输达到地表的部分叫做大气向下散射辐射(I_s)。同时，达到地表的各种辐射又将受到进一步各种复杂环境反射，一部分回馈给地表，这部分统称为环境反射辐射(I_e)。因此对于阴影区，到达区内的辐射主要为I_s和I_e；而对于非阴影区，到达区内的辐射除了I_s和I_e外，还有最主要的直接辐射I_d(图 1)。这些辐射到达地表后，被进一步反射，并再次通过大气层被传感器接收。

点击放大

下载原图

图 1 城市建筑阴影区及其附近非阴影区

Figure 1Main radiance transfer processes of an urban shadow and its nearby sunlit area

利用遥感数据中记录的辐射信号反演地表真实反射率是定量遥感应用过程中的重要前提，这个过程也被地表反射率定标，而大气校正是其中的重要内容。该过程可通过式(1)表示(Liang等，2001)：

${\rho _{sun}} = \frac{{\pi \left( {{L_{sun}} - {L_P}} \right)}}{{{T_V}\left( {{I_d} + {I_s} + {I_e}} \right)}} = \frac{{\pi \left( {{L_{sun}} - {L_P}} \right)}}{{{T_V}{I_{sun}}}}$ (1)

式中，ρ_sun表示非阴影区地表反射率，L_sun表示非阴影区表观辐亮度，I_p表示大气程辐射，T_v表示沿卫星观测方向大气辐射透过率，I_sun表示达到非阴影区地表总辐射，即I_d+I_s+I_e。

同理，对于阴影地表(无直接辐射I_d到达)反射率定标而言，亦可借鉴式(1)基本原理，即：

${\rho _{shw}} = \frac{{\pi \left( {{L_{shw}} - {L_P}} \right)}}{{{T_V}\left( {{I_s} + {I_e}} \right)}} = \frac{{\pi \left( {{L_{shw}} - {L_P}} \right)}}{{{T_V}{I_{shw}}}}$ (2)

式中，ρ_shw表示阴影区地表反射率，L_shw表示阴影区表观辐亮度，I_shw表示达到阴影区地表总辐射，即I_s+I_e。

若阴影区和临近非阴影区内同时包含某种相同地物，如城市绿化草坪、城市路面等(表示为“*”)，则在朗伯体地表假设下，该地物(*)在阴影区内的反射率与其在临近非阴影区的反射率相等(即$\rho _{_{shw}}^*$=$\rho _{_{sun}}^*$)，据此可得到方程组式(3)。

$\left\{ \begin{array}{l} \rho _{shw}^* = \frac{{\pi \left( {L_{shw}^* - {L_P}} \right)}}{{{T_V}{I_{shw}}}}\\ \rho _{sun}^*\frac{{\pi \left( {L_{sun}^* - {L_P}} \right)}}{{{T_V}{I_{sun}}}}\\ \rho _{shw}^* = \rho _{sun}^* \end{array} \right.$ (3)

式中，$\rho _{_{sun}}^*$和${L_{_{sun}}^*}$分别表示地物(*)在非阴影区像元的地表反射率和表观辐射亮度，$\rho _{_{shw}}^*$和${L_{_{shw}}^*}$分别表示该地物在阴影区像元的地表反射率和表观辐射亮度。通过对式(3)进一步整理可得：

$\frac{1}{{{I_{shw}}}} = \frac{{L_{sun}^* - {L_P}}}{{{I_{sun}}\left( {L_{shw}^* - {L_P}} \right)}}$ (4)

将式(4)代入式(2)即得到：

${\rho _{shw}} = \frac{{\pi fc\left({{L_{shw}} - {L_P}} \right)}}{{{T_V}{I_{sun}}}}$ (5)

式中，f_c表示为：

${f_c} = \frac{{L_{sun}^* - {L_P}}}{{L_{shw}^* - {L_P}}}$ (6)

由此，可通过式(5)和式(6)反演得到阴影区地表反射率，实现阴影校正。尽管如此，这种直接恢复阴影区地表反射率的阴影校正策略仍存在不足：需要的未知参数较多，包含f_c、I_p、T_v及I_sun。这些参数中，I_p、T_v及I_sun亦是大气校正过程需要的主要参数，往往很难精确估计。虽然通过辐射传输软件包如6S可以模拟计算，但是这种模拟方法需要更多其他同步获取参数，不利用推广使用。为此，本研究将尝试在只利用影像自身信息的情况下，实现阴影校正。我们采用的策略是：首先在尽量较少未知参数的前提下，将阴影校正后的变量变为表观辐射亮度，即可达到阴影校正的目的。

对式(5)进一步变形得到：

${I_{sun}}{\rho _{shw}}{T_v} = \pi fc\left({{L_{shw}} - {L_P}} \right)= L_{shw}^{hyp}$ (7)

式(7)的左边部分表示了如下虚拟过程：非阴影区直射辐射(I_sun)达到阴影区，并被阴影区地表反射(ρ_shw)，经过大气(T_v)被传感器接收($L_{shw}^{hyp}$)。众所周知，阴影区在现实中是不可能有直射辐射(光)到达的；若有，则证明阴影区已不存在，即阴影被校正(去除)。因此，这里的$L_{shw}^{hyp}$应理解为传感器接收到来自校正后的阴影区的辐射值。同时，传感器端接收的总辐射(即表观辐亮度)，还应包括大气程辐射(I_p)部分：

$L_{shw}^{rst} = \pi fc\left({{L_{shw}} - {L_P}} \right)+ {L_P}$ (8)

式中，$L_{shw}^{cor}$表示阴影区校正后的表观辐亮度。

采用式(6)与式(8)实现阴影校正，阴影校正结果变量为表观辐亮度，校正过程中未知变量只有I_p与f_c。由于该方法是建立在阴影区与非阴影区同一地物反射率相等这一基本假设上的，我们用RERB(Reflectance Equality Relationship Based)表示该方法。