综上所述，利用数字高程模型(DEM)或者数字表面模型(DSM)结合卫星传感器所接收的实际辐射观测值来获取地表真实参数是地形辐射校正的一般思路[6]。但是，目前几乎所有的模型及实验都主要针对于中低分辨率的影像，例如MODIS、TM等影像，对于高分辨率处理相对较少。目前，随着高分辨率遥感卫星发射日益增多，应用需求也越来越多，而地形效应在高分辨率影像的表现更加明显，对辐射信号影响更大。传统的地形辐射校正方法因对分辨率提高所带来的误差源控制以及模型建模等因素考虑不足，对于高分辨率影像适用性存在一定问题。究其原因，除了高精度、高分辨率的DEM(或DSM)不易获取是主要限制外，针对高分辨率影像的地形辐射校正模型和方法自身也存在一定问题。与此同时，相当数量遥感卫星针对全色影像并不提供绝对辐射定标系数，因此如何在没有定标系数的条件下实现对高分辨率全色影像进行地形辐射校正也是亟待解决的问题。
针对上述问题，本文重点利用辐射传输模型和更高分辨率(10 m)的DSM，建立了一种面向高分辨率遥感影像的地形辐射校正方法。通过利用北京北部山区的资源三号01星高分辨率多光谱和全色遥感影像进行地形辐射校正实验。实验结果表明，本文方法优于传统方法，更适合于高分辨率影像地形辐射校正。
1 地形辐射校正误差源分析 众多模型中，基于物理模型的地形辐射校正的实质是地形因子和大气参数同步解算，在建模过程中需综合考虑各种影响因素，考虑的影响因素越多，输入的参数(例如绝对辐射定标系数、大气参数等)越真实，中间计算的过程(例如DEM与影像分辨率越接近，两者配准越精确)越合理，最终的校正结果就越准确。因此，影响高分辨率影像地形辐射校正精度的主要因素可概括为2个：①校正模型；②输入数据及参数的误差。
1.1 校正模型分析 对于复杂地形而言，地面一点所接收的有效光照可分解为3项：①太阳的直射辐照度Ed；②天空漫散射辐照度Ef；③周围邻近山坡辐照度Ea。复杂地形地表总入射太阳辐照度可表示为
 | (1) |
如果模型中只考虑前一或两项辐射，计算较为简单，但效果并不好。原因在于当目标位于深谷或阴影区时，周围邻近山坡辐射项Ea的辐射量可能会在一点所接收的总辐射量中占有较大的比重，而模型却没有如实反映。除此之外，研究发现，当地表有雪覆盖或是近红外波段有植被覆盖时，邻近地形的反射辐射也不可忽略[7]。这些都充分说明地形辐射校正模型中仅考虑太阳直接辐射是不够的[8]。
本文在建模上全面考虑太阳直接辐射、天空漫散射辐射和周围邻近山坡辐射等所有辐射项，以真实反映辐射能量分布的实际情况。
1.2 输入数据及参数的误差分析 根据物理模型，输入的误差源主要表现为DEM(DSM)、大气辐射传输计算和地表特性3类。具体如下：
1) 计算地形因子的DEM或DSM
对于高分辨率遥感影像而言，DSM包含地物高程信息，比DEM更具备适应性，而且DSM(或DEM)应与高分辨率影像相当或稍低[9]，如果相差较大，则DSM无法准确描述真实地形起伏情况。同时，无论用什么方法，DSM与高分辨率影像的配准精度也非常重要，这些将决定后续计算地形因子的精度。
2) 大气辐射传输计算
大气参数估算是地形辐射校正中非常重要的一步，如果有条件最好能够针对实验区域同步(或准同步)实测大气参数[10]，减少估算误差，最大限度提高地形辐射校正的精度。
3) 地表特性
地表特性的假设是地形辐射校正中的一项重要误差源。很多地形辐射校正模型是基于朗伯体假设，把地表理想化为均匀、各向同性的表面，而实际地表是具有明显方向性的非朗伯体。可采用常数K简化描述地表二向性分布函数[11]。
以上3个误差源的分析表明，DSM的选取及配准精度影响最大；其次是地表特性描述准确程度；最后是大气参数精度。研究表明：大气参数误差对地形辐射校正影响相对较小，特别是天空漫散射和程辐射, 并且在不同大气模式下, 大气参数误差引起地表反射率的反演误差相差不大[12]。因此，针对高分辨率影像的地形辐射校正，需要从上述源头控制好误差，才能提升高分辨率影像的地形辐射校正精度。本文通过控制这些误差，研究了一种适用于高分辨率影像的地形辐射校正方法。
2 地形辐射校正方法 2.1 光照分解模型 大部分地物都是由3部分辐射综合作用经地表目标的方向反射进入传感器而被观测。根据Sandmeier和Itten[13]的模型提出，在太阳有效入射角is和天顶角θs、地表目标所处地形坡度S和坡向A下，山区地表目标接收的总入射辐射照度EinTotal为
 | (2) |
式中：b(x, y)为二进制阴影遮蔽因子，用来判断像元是否接收到直射辐射，当处于非阴影区时b=1，处于阴影区时b=0；Vf为天空开阔度因子，定义为从坡面像元观察的天空区域与在无遮挡的水平面上观察的天空区域之比，取值为0~1之间；等式右边第1项表示太阳直射辐照度，Edh为平坦地表接收的太阳直射辐照度；等式右边第2项和第3项分别表示环日各向异性散射辐照度与天空各向同性散射辐照度，两者综合作用为天空漫散射辐照度Ef，Efh为平坦地表接收的天空漫散射辐照度，K为各向异性指数，表示环日各向异性散射占天空漫散射的权重，用地面法线方向接收的太阳直射辐照度与大气层顶太阳辐照度之比计算；等式右边第4项表示周围地形的反射辐射Ea，针对高分辨率遥感影像，此项辐射不可忽略。
2.2 多光谱遥感影像地形辐射校正模型 根据光照分解模型，推导出多光谱遥感影像地形辐射校正模型如下：
 | (3) |
式中：ρ(x, y)为地物反射率；DN为卫星影像灰度值；DNba为程辐射对应的图像DN值，在无法准确估算的情况下，可以采用图像有效像元区域最小值来表示；c0和c1为影像绝对辐射定标系数；Eiter为临近像元的多次反射辐照度；模型中与大气相关的因子τv、τs、Ed和Edif分别为大气上行透射率、下行透射率、地面直射辐照度和地面散射辐照度，可直接由6SV软件计算得到。
模型中与地形相关的因子有遮蔽因子b(x, y)、天空开阔度因子Vf、太阳有效入射角is以及邻近地物反射的天空观测因子Vsky，所有这些地形相关因子均可由DSM计算的坡度、坡向等地形参数获取。其中太阳有效入射角is公式为
 | (4) |
式中：?s为太阳方位角。坡度S、坡向A是通过DEM或者DSM计算出来的。邻近像元的多次反射辐照度Eiter的计算公式如下：
 | (5) |
式中：Eg为水平表面接收的太阳的直射辐照度与漫散射辐照度之和；Vsky=0.5+0.5cos S; ρterrain(i-1)表示像元周边一定范围内(一般0.5 km)的平均反射率，属于中间结果。多次反射辐照度Eiter的计算通常迭代2~3次即可。
2.3 全色遥感影像地形辐射校正模型 全色遥感影像的地形辐射校正方法基于严格的物理模型，综合了辐射传输模型与经验系数校正的优点。如果全色遥感影像在有绝对辐射定标系数前提下，上述多光谱遥感影像地形辐射校正模型也适用于全色。但实际上很多型号卫星并不提供全色遥感影像的定标系数，因而无法得到绝对地表反射率。例如HJ-1[14]和CBERS-04等卫星，以及GF-2卫星的PSM2载荷等。因此，对于这些定标系数缺失的全色遥感影像，本文提出了如下校正方法计算相对反射率实现地形辐射校正：
 | (6) |
式中：DNba可选择图像有效像元区域最小值代表该数值；Vf由坡度计算得到。
令ξ=(1-Vf)/10，由于无法得到绝对地表反射率，难以进行多次迭代，可采用式(7)简化计算：
 | (7) |
本文针对全色遥感影像的地形辐射校正采用了3种方法，即绝对辐射定标系数法(即多光谱遥感影像地形辐射校正模型)、缺少定标系数的相对反射率法(即本节模型)以及相对反射率后带掩模优化法。其中，相对反射率后带掩模优化法主要考虑到影像分辨率与所采用DEM分辨率差异大，通过加掩模有利于补偿校正效果。
2.4 处理步骤及流程 步骤1??重投影、正射校正与配准。
按照影像的投影方式对DSM进行重新投影，确保两者投影一致。利用有理函数模型(RPC)参数及DSM数据对一级影像进行正射校正，获取正射影像。通过正射校正，确保影像和DSM数据之间配准精度。
步骤2??数据预处理。
1) 裁剪。对DSM和影像的重叠区域进行裁剪，获取对应相同范围的DSM数据和影像数据。
2) 辐射定标。将绝对辐射定标系数应用到多光谱及全色遥感影像中，得到辐亮度影像(单位为W·m-2·sr-1·μm-1)。
步骤3??地形辐射校正。
地形辐射校正主要分3步，可采用图 1所示流程图表示。
 |
| 图 1 地形辐射校正流程图 Fig. 1 Topographic correction flowchart |
| 图选项 |
1) 地形相关因子计算。通过DSM完成坡度、坡向、天空观测因子、遮蔽因子计算。
2) 大气参数计算。输入实测气象数据，通过6S模型计算大气辐射传输量输出项，分别是：到达地面直射辐照度、到达地面漫射辐射辐照度、日地平均距离处太阳辐照度、程辐射、上行透射率、半球反照率。
3) 地表反射率计算。基于辐射传输参量计算校正后的地表反射率。
3 地形辐射校正实验 3.1 实验数据
3.1.1 图像数据 实验区域位于北京北部以十三陵水库为中心的山区，该区域以山地为主，海拔800~1 000 m。实验数据为2016年5月16日成像的资源三号01星多光谱遥感影像及全色遥感影像共2景，以DSM数据的实际分布情况作为特征区域的选择依据，本实验在多光谱、全色2景遥感影像中分别选取10个特征区域作为样本影像。其中，多光谱遥感影像分辨率为5.8 m，全色遥感影像分辨率为2.1 m。影像中心经度为116.89°，中心纬度为40.03°。数据格式为Geotiff，影像产品级别为1级加上配套RPC参数文件，可自行处理生成2级或正射影像。配套的地形数据是资源三号01星立体像对数据通过摄影测量方法生成的10 m分辨率的DSM数据(见图 2)。
 |
| 图 2 实验图像示例 Fig. 2 Example of experimental images |
| 图选项 |
3.1.2 辅助数据 1) 绝对辐射定标系数
多光谱遥感影像的绝对辐射定标系数来自于资源卫星中心网站，由于影像数据获取时间为2016年5月16日10:58:43，采用最新2015年提供的绝对辐射定标系数；全色遥感影像2015年没有绝对辐射定标系数，但2014年提供了绝对辐射定标系数。全色与多光谱遥感影像的绝对辐射定标系数见表 1。
表 1 资源三号01星遥感影像绝对辐射定标系数 Table 1 Absolute radiometric calibration coefficients for ZY3-01 satellite images
| 载荷 | 波段 | 光谱范围/ μm | 定标系数增益项 | 定标系数偏移量 |
| 多光谱 | Band 1 | 0.45~0.52 | 0.233 0 | 0 |
| Band 2 | 0.52~0.59 | 0.216 2 | 0 | |
| Band 3 | 0.63~0.69 | 0.178 9 | 0 | |
| Band 4 | 0.77~0.89 | 0.194 9 | 0 | |
| 全色 | Pan | 0.50~0.80 | 0.170 8 | 0 |
表选项
利用绝对辐射定标系数将卫星影像DN值转换为辐亮度图像的公式为
 | (8) |
式中：Le(λe)为转换后辐亮度，W·m-2·sr-1·μm-1；G为定标斜率；Offset为绝对辐射定标系数偏移量，空缺值为0。
2) 光谱响应函数
在资源中心网站查询多光谱波段及全色数据的光谱响应函数。
3) 大气数据
影像区域含有2个AERONET观测站，分别为Beijing站(图 3中绿点所示)，经度116.381°，纬度39.977°，以及Beijing_CAMS站(图 3中红点所示)，经度116.381°，纬度39.933°。
 |
| 图 3 观测站位置示意图 Fig. 3 Schematic of location of observation station |
| 图选项 |
取成像时刻前后15 min平均值作为图像大气数据。根据Angstrom公式，分别计算2个站点550 nm气溶胶光学厚度，然后求均值得到AOT550=0.156，水汽含量为0.912 g/cm2，可以发现影像获取时刻区域内能见度较高，气溶胶影响相对较小。
3.2 处理结果及评价
3.2.1 处理结果对比 1) 多光谱遥感影像处理前后辐射对比(见图 4)。
 |
| 图 4 多光谱遥感影像处理结果辐射对比 Fig. 4 Comparison of radiation of multispectral image experimental process results |
| 图选项 |
2) 多光谱遥感影像处理前后细节对比(见图 5)。
 |
| 图 5 多光谱遥感影像处理细节对比 Fig. 5 Comparison of process results for details of multispectral images |
| 图选项 |
3) 全色遥感影像处理前后辐射对比(见图 6~图 8)。
 |
| 图 6 全色遥感影像处理结果辐射对比 Fig. 6 Comparison of radiation of panchromatic image process results |
| 图选项 |
 |
| 图 7 针对全色遥感影像的不同方法结果对比 Fig. 7 Comparison of different methods for panchromatic images |
| 图选项 |
 |
| 图 8 不同方法的全色遥感影像校正结果对比 Fig. 8 Comparison of correction results for panchromatic images among different methods |
| 图选项 |
3.2.2 结果评价与分析 针对地形辐射校正的结果评价一般分为主观评价和客观评价。主观评价以目视评价为主；客观评价以校正前后客观指标值变化进行量化评判，主要包含2部分，一是统计特征评价，例如直方图、标准差等，二是光谱一致性评价，例如相关系数、变异系数等。
1) 目视评价
从处理前后不同区域的辐射对比图(见图 4~图 8)中可见，校正后影像的立体感变弱，阴坡区域辐亮度提升，说明本文方法较好地消除了地形起伏带来的像元辐亮度畸变或地物反射率畸变。其中对于高分辨率多光谱遥感影像，位于阴、阳坡的同一类地物光谱差异明显变小，同时较好地保留了影像细节。对于全色遥感影像，本文分别采用3种方法处理，绝对辐射定标系数法采用的是2014年度提供的定标系数。通过3种方法校正结果对比，绝对辐射定标系数法在某些阴影区域信息并未恢复出来，原因可能是2014年的绝对辐射定标系数不能适用于全色遥感影像。相对反射率法部分恢复阴影信息，相对反射率后带掩模优化法则较好地恢复了阴影区域信息。
实验中同时还采用了传统的C校正和商业软件ATCOR针对相同的数据源进行处理并对比。其中C校正结果稍差，体现为细节缺失；ATCOR算法与本文方法结果基本相当，但是ATCOR算法不能处理缺少绝对辐射定标系数的图像，且对于图像分辨率远高于DSM分辨率的情况处理效果不甚理想。本文方法具有一定的普适性。
2) 统计特征评价
进行统计特征评价前，首先需要统一对比基准，本实验采用地表未经地形辐射校正的反射率影像作为处理前影像，与经地形辐射校正后的反射率结果进行对比，统计特征的实验区域如图 9所示。对比的统计特征值分别为直方图、均值、标准差、变异系数以及X、Y方向剖面图。
 |
| 图 9 统计特征区域示意图 Fig. 9 Schematic of statistical characteristic area |
| 图选项 |
由图 10可以看出，校正前存在不止一个波峰，说明存在坡度、坡向影响，校正后呈高斯分布，说明与自然界地物随机性一致。
 |
| 图 10 多光谱遥感影像处理前后直方图统计对比 Fig. 10 Statistic comparison of histogram for multispectral image process results |
| 图选项 |
由表 2可以看出，校正后标准差较校正前有所降低，整体均值没有明显变化且略有提升，说明校正了地形效应影响，同时较好保持了影像细节信息。这里同时引用变异系数[15]，定义为标准差与均值的比值，一般如果处理后有效果，该值应该变小，从统计结果也可看到，处理后每个波段的变异系数均有所降低。
表 2 多光谱遥感影像校正前后统计特征值对比 Table 2 Comparison of statistical eigenvalues before and after correction for multispectral images
| 波段 | 校正前 | 校正后 | |||||
| 均值/ 10-2 | 标准差/ 10-2 | 变异系数 | 均值/ 10-2 | 标准差/ 10-2 | 变异系数 | ||
| 1 | 0.407 | 0.420 | 1.032 | 0.505 | 0.420 | 0.832 | |
| 2 | 2.339 | 0.832 | 0.356 | 2.529 | 0.757 | 0.299 | |
| 3 | 2.400 | 0.655 | 0.273 | 2.630 | 0.579 | 0.220 | |
| 4 | 29.049 | 5.722 | 0.197 | 31.469 | 5.107 | 0.162 | |
表选项
由表 3可以看出，与多光谱遥感影像类似，全色遥感影像经过不同方法处理后，标准差及变异系数都有所降低，整体均值都略有提升，符合地形效应影响减弱或消除的特征，说明不同方法针对全色遥感影像都有效。其中从统计特征中的变异系数项看，本文方法结果稍优于C校正和ATCOR算法。需要说明的是，全色遥感影像对比的统一基准是相对反射率影像。
表 3 全色遥感影像校正前后统计特征值对比 Table 3 Comparison of statistical eigenvalues before and after correction for panchromatic images
| 统计量 | 原始图像 | C校正法 | ATCOR算法 | 本文方法 |
| 均值/10-2 | 7.59 | 10.07 | 8.28 | 9.55 |
| 标准差/10-2 | 4.36 | 2.98 | 2.03 | 2.10 |
| 变异系数 | 0.574 | 0.296 | 0.245 | 0.220 |
表选项
分析结果如图 11所示, 图中白色曲线为校正前剖面曲线，红色曲线为校正后剖面曲线。由X和Y方向剖面线数值变化可看出，经过地形辐射校正后的区域数值幅度明显高于校正前数值，这说明本文方法的有效性。
 |
| 图 11 多光谱遥感影像处理前后X与Y方向剖面对比 Fig. 11 Comparison of X and Y direction profiles before and after correction for multispectral image |
| 图选项 |
3) 光谱一致性评价
对于高分辨率多光谱遥感影像，可采用相似度即相关系数指标对实验区影像进行统计，以分析影像波谱响应的相似程度。具体方法是：校正前后影像中，先选取若干个(5~10个)阴坡、阳坡区域，各自构成区域对(见图 12(a))，分别计算校正前后相邻阴坡、阳坡区域对中各自点集光谱曲线的相关系数。除此之外，也对处理前后相同位置的阴阳坡点对的光谱线离散程度进行对比。
 |
| 图 12 光谱一致性评价图 Fig. 12 Evaluation chart of spectral consistency |
| 图选项 |
从图 12(b)、(c)的阴阳坡光谱曲线图可看到，校正后的曲线离散度明显小于校正前。从平均相关系数来看，尽管校正前后都接近于1，但校正后仍有0.3%的提升，说明地形辐射校正减弱了地形对地物信息遮挡的效应，从整体上改善了图像质量。
4 结论 本文利用资源三号01星全色及多光谱遥感影像，通过相应DSM数据、实时观测大气数据以及绝对辐射定标系数等配套输入数据，反演地表反射率或相对反射率，探讨针对高分辨率影像的地形辐射校正方法。
通过对实验结果评价分析可看出，本文方法实现了高分辨率影像的地形辐射精确校正，影像的纹理细节保真效果较好。同时，本文也针对缺少绝对辐射定标系数的高分辨率全色影像的地形辐射校正提出了新的方法，并进行了不同方法效果的对比，实验结果表明该方法有效。
参考文献
| [1] | COLBY J D. Topographic normalization in rugged terrain[J].Photogrammetric Engineering and Remote Sensing, 1991, 57: 531–537. |
| [2] | 陈志明. 遥感影像地形辐射校正方法研究与系统实现[D]. 福州: 福建师范大学, 2009: 4. CHEN Z M. Research on topographic correction methods for remote sensing image and system implementation[D]. Fuzhou: Fujian Normal University, 2009: 4(in Chinese).http://cdmd.cnki.com.cn/Article/CDMD-10394-2010035006.htm |
| [3] | 武瑞东. 卫星遥感影像数据的地形影响校正[J].遥感信息, 2005(4): 31–34. WU R D. Topographic correction of satellite remote sensing image data[J].Remote Sensing Information, 2005(4): 31–34.(in Chinese) |
| [4] | SHEPHERD J D, DYMOND J R. Correcting satellite imagery for the variance of reflectance and illumination with topography[J].International Journal of Remote Sensing, 2003, 24(17): 3503–3514.DOI:10.1080/01431160210154029 |
| [5] | RICHTER R. Correction of atmospheric and topographic effects for high spatial resolution satellite imagery[J].International Journal of Remote Sensing, 1997, 18(5): 1099–1111.DOI:10.1080/014311697218593 |
| [6] | 段四波, 阎广建. 山区遥感图像地形辐射校正模型研究综述[J].北京师范大学学报(自然科学版), 2007, 43(3): 362–366. DUAN S B, YAN G J. A review of models for topographic correction of remotely sensed images in mountainous area[J].Journal of Beijing Normal University(Natural Science), 2007, 43(3): 362–366.(in Chinese) |
| [7] | PROY C, TANRE D, DESCHAMPS P Y. Evaluation of topographic effects in remotely sensed data[J].Remote Sensing of Environment, 1989, 30(1): 21–32.DOI:10.1016/0034-4257(89)90044-8 |
| [8] | 高永年, 张万昌. 遥感影像地形辐射校正研究进展及其比较试验[J].地理研究, 2008, 27(2): 467–468. GAO Y N, ZHANG W C. Comparison test and research progress of topographic correction on remotely sensed data[J].Geographical Research, 2008, 27(2): 467–468.(in Chinese) |
| [9] | ZHANG Y, YAN G, BAI Y. Sensitivity of topographic correction to the DEM spatial scale[J].IEEE Geoscience Remote Sensing Letters, 2015, 12(1): 53–57.DOI:10.1109/LGRS.2014.2326000 |
| [10] | SOLA I, GONZáLEZ-AUDíCANA M, áLVAREZ-MOZOS J. Validation of a simplified model to generate multispectral synthetic images[J].Remote Sensing, 2015, 7(3): 2942–2951.DOI:10.3390/rs70302942 |
| [11] | VANONCKELEN S, LHERMITTE S, VAN ROMPAEY A. The effect of atmospheric and topographic correction on pixel-based image composites:Improved forest cover detection in mountain environments[J].International Journal of Applied Earth Observation and Geoinformaiton, 2015, 35B: 320–328. |
| [12] | 段四波, 阎广建, 穆西晗, 等. 基于DEM的山区遥感图像地形辐射校正方法[J].地理与地理信息科学, 2007, 23(6): 23–24. DUAN S B, YAN G J, MU X H, et al. DEM based remotely sensed imagery topographic correction method in mountainous areas[J].Geography and Geo-Information Science, 2007, 23(6): 23–24.(in Chinese) |
| [13] | SANDMEIER S, ITTEN K I. A physically-based model to correct atmospheric and illumination effects in optical satellite data of rugged terrain[J].IEEE Transactions on Geoscience and Remote Sensing, 1997, 35(3): 708–717.DOI:10.1109/36.581991 |
| [14] | JIANG B, LIANG S, TOWNSHEND J R, et al. Assessment of the radiometric performance of Chinese HJ-1 satellite CCD instruments[J].IEEE Journal of Selected Topics in Applied Earth Observation & Remote Sensing, 2013, 6(2): 840–850. |
| [15] | 罗良清, 魏和清. 统计学[M].北京: 中国财政经济出版社, 2011: 49. LUO L Q, WEI H Q. Statistics[M].Beijing: China Financial and Economic Publishing House, 2011: 49.(in Chinese) |
