CN104834814B - 遥感影像地形标准化方法 - Google Patents
遥感影像地形标准化方法 Download PDFInfo
- Publication number
- CN104834814B CN104834814B CN201510211271.6A CN201510211271A CN104834814B CN 104834814 B CN104834814 B CN 104834814B CN 201510211271 A CN201510211271 A CN 201510211271A CN 104834814 B CN104834814 B CN 104834814B
- Authority
- CN
- China
- Prior art keywords
- terrain
- radiation
- atmospheric
- reflectivity
- formula
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 13
- 239000000443 aerosol Substances 0.000 claims abstract description 27
- 230000006870 function Effects 0.000 claims abstract description 20
- 230000003287 optical effect Effects 0.000 claims abstract description 12
- 230000008859 change Effects 0.000 claims abstract description 9
- 230000005855 radiation Effects 0.000 claims description 68
- 238000002310 reflectometry Methods 0.000 claims description 47
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Chemical compound O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims description 37
- 238000012937 correction Methods 0.000 claims description 27
- 238000002834 transmittance Methods 0.000 claims description 24
- 230000000694 effects Effects 0.000 claims description 18
- 238000001228 spectrum Methods 0.000 claims description 18
- 238000010606 normalization Methods 0.000 claims description 15
- CBENFWSGALASAD-UHFFFAOYSA-N Ozone Chemical compound [O-][O+]=O CBENFWSGALASAD-UHFFFAOYSA-N 0.000 claims description 12
- 238000010521 absorption reaction Methods 0.000 claims description 10
- 238000009826 distribution Methods 0.000 claims description 10
- 238000011425 standardization method Methods 0.000 claims description 9
- 230000005540 biological transmission Effects 0.000 claims description 8
- 230000003595 spectral effect Effects 0.000 claims description 8
- 238000005286 illumination Methods 0.000 claims description 6
- 230000009471 action Effects 0.000 claims description 2
- 230000000873 masking effect Effects 0.000 abstract 1
- 238000004422 calculation algorithm Methods 0.000 description 8
- 238000004364 calculation method Methods 0.000 description 6
- 238000011160 research Methods 0.000 description 5
- FEIWNULTQYHCDN-UHFFFAOYSA-N mbba Chemical compound C1=CC(CCCC)=CC=C1N=CC1=CC=C(OC)C=C1 FEIWNULTQYHCDN-UHFFFAOYSA-N 0.000 description 4
- RWIUTHWKQHRQNP-NQLNTKRDSA-N Melinamide Chemical compound CCCCC\C=C/C\C=C/CCCCCCCC(=O)NC(C)C1=CC=CC=C1 RWIUTHWKQHRQNP-NQLNTKRDSA-N 0.000 description 3
- 238000010586 diagram Methods 0.000 description 3
- 238000012876 topography Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000008901 benefit Effects 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 244000025254 Cannabis sativa Species 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 230000008878 coupling Effects 0.000 description 1
- 238000010168 coupling process Methods 0.000 description 1
- 238000005859 coupling reaction Methods 0.000 description 1
- 238000009795 derivation Methods 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 230000008689 nuclear function Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000011158 quantitative evaluation Methods 0.000 description 1
- 230000005477 standard model Effects 0.000 description 1
- 238000011144 upstream manufacturing Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Landscapes
- Image Processing (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
一种遥感影像地形标准化方法，在Ross Thick‑Li Sparse线性核驱动模型中，将坡度坡向引入Ross体散射核函数与Li氏几何光学核函数，建立山区BRDF模型；在现有遥感影像地形标准化成果基础上，改变地表朗伯反射特征假设，引入山区BRDF模型，最大限度减少复杂地形区BRDF形状的改变对遥感数据的影响；将MODIS水汽和气溶胶大气产品作为地形标准化模型参数，获得卫星影像过境时刻的大气参数，消除大气光学因素的影响。该方法不仅去除了大气光学特性对遥感影像的辐射畸变，同时有效去除地形遮蔽等因子、地表BRDF特性引起的辐射畸变，从而获得地表真实反射率。
Description
技术领域
本发明属于遥感技术领域，涉及一种遥感影像地形标准化方法。
背景技术
地表反照率是影响地表短波净辐射的重要因素，高分辨率遥感影像成为反照率（Albedo）反演的重要数据源。然而，遥感观测不仅受水汽、气溶胶等大气光学因素影响，同时还受坡度、坡向、地形遮蔽、天空可视因子等地形因子以及地表二向反射分部函数（BRDF）或太阳-地表-传感器几何关系的影响，不能真实地反映地表状况。因此，遥感影像地形标准化已经成为地表Albedo精确反演的必要预处理手段，在高分辨率数字高程模型（DEM）基础上，用于同步消除或减小大气效应、地形效应和地表BRDF效应。
发明内容
本发明的目的是提供一种遥感影像地形标准化方法，消除或减小大气效应、地形效应和地表BRDF效应，真实反映地表状况。
为实现上述目的，本发明所采用的技术方案是：一种遥感影像地形标准化方法，其特征在于，该方法具体为：
1）卫星传感器获得的大气顶辐射亮度L TOP (λ)由路径辐射L p (λ)和来自地表反射的辐射L(λ)，即：
L TOP (λ)= L p (λ)+ L(λ)T(λ,θ ν ) (1)
在山区，来自地表反射辐射为：
L(λ)=(1/π)(E dir (λ)+ E aniso-dif (λ))×ρ T (λ)(i s ,i ν ,Φ s-ν )+ (1/π) (E iso-dif (λ)+E ref (λ))ρ Tsd (λ) (i ν ,Φ ν ) (2)
式中：E dir (λ)、E iso-dif (λ)、E aniso-dif (λ)和E ref (λ)分别表示太阳直接辐射、各向同性散射辐射、各向异性散射辐射和来自周围地形的反射辐射；ρ T (λ)和ρ Tsd (λ)分别表示坡地方向-方向反射率和半球-方向反射率；
对于相同地表覆盖类型而言，坡地与消除地形的平坦地表方向-方向反射率不同，其关系为：
ρ T (λ)( i s , i v ,Φ s-v )=Ω(λ) ( i s , i v ,Φ s-v ,θ s ,θ v ,φ s-v ) ×ρ H (λ) (θ s ,θ v ,φ s-v ) (5)
式中，Ω(λ) ( i s , i v ,Φ s-v ,θ s ,θ v ,φ s-v )为坡地BRDF归一化方向反射因子；由于半球-方向反射率比较复杂，同时坡地各向同性散射辐射和来自周围地形反射辐射的能量较小，因此不区分坡地和平坦地表半球-方向反射率，即：
ρ Tsd (λ)(i v ,Φ v )=ρ H (λ) (θ s ,θ v ,φ s-v ) (6)
因此，将公式（2）、公式（5）和公式（6）代入公式（1），地形标准化后的反射率为：
ρ H (λ) (θ s ,θ v ,φ s-v )=π(L TOP (λ)- L p (λ))/T(λ,θ v )〔(E dir (λ)+ E aniso-dif (λ)) Ω(λ)+E iso-dif (λ) + E ref (λ)〕 (7)
当Ω(λ)=1时，地表假设为朗伯体反射特征；存在两种情况，一种是对遥感影像只进行大气校正，公式可变为：
ρ H (λ) (θ s ,θ v ,φ s-v )=π(L TOP (λ)- L p (λ))/T(λ,θ v )E 0 (λ)D 0 T(λ,θ s ) (8)
另外一种情况是对遥感影像同时进行大气校正与地形校正，即地形标准化，但只是去除地形对太阳照射角的影响，则地表反射率反演公式为：
ρ H (λ) (θ s ,θ v ,φ s-v )=π(L TOP (λ)- L p (λ))/T(λ,θ v )E(λ) (9)
2）大气透过率受气溶胶、水汽含量的影响较大，同时也受瑞利散射和臭氧吸收的影响，即：
T(λ,θ s )= T o3 (λ,θ s )T w (λ,θ s ) T r (λ,θ s )T a (λ,θ s ) (10)
式中：T(λ,θ s )、T o3 (λ,θ s )、T w (λ,θ s )、T r (λ,θ s )和T a (λ,θ s )分别表示总的大气波谱透过率、臭氧透过率、水汽透过率、瑞利散射透过率和气溶胶透过率。由于臭氧吸收与瑞利散射随时空分布变化不大，臭氧光学厚度可采用多年平均的空间分布数据，瑞利散射可利用经验公式估算；然而在山区，气溶胶和水汽分布随时空变化较为强烈，水汽和气溶胶波谱透过率公式如下：
式中，k(λ)、m和w分别表示水汽有效吸收系数、大气质量及大气可降水厚度；α和β是angstrom波长指数与浑浊度；对于大气可降水厚度，MODIS水汽产品PW直接作为水汽透过率输入参数；AOD产品中的波长指数也作为气溶胶透过率的输入参数；β表示为：
Β=τ 0.55(0.55) a (14)
在山区，总的太阳波谱辐射由太阳直接辐射、太阳散射和周围地形反射辐射三部分组成，均为波长的函数；考虑到太阳方向对散射辐射的影响，将其又分解为各向同性散射辐射和各向异性散射辐射，即：
E(λ)= E dir (λ)+E iso-dif (λ)+E aniso-dif (λ)+E ref (λ) (15)
对于Landsat TM每个波段的太阳直接辐射与散射辐射同时受大气与地形作用的影响；
周围地形反射辐射估算方法如下：
式中，F ij 为地形结构因子，表示来自周围可见像元的反射辐射能量与可见像元自身总的太阳波谱辐射； L i (λ)是周围可见像元的辐射亮度，由Landsat TM影像测量而得；
3）坡地BRDF 归一化方向反射因子是遥感影像地形标准化的重要参数，是坡地与消除地形的平坦地表方向-方向反射率之比，而这两种方向反射率均遵循地表BRDF规律；半经验线性核驱动模型由一定物理意义核的线性组合来拟合地表二向反射特性，其公式如下：
ρ H (λ) (θ s ,θ v ,φ s-v )=f iso (λ)+ f vol (λ)K vol (θ s ,θ v ,φ s-v )+ f geo (λ)K geo (θ s ,θ v ,φ s-v ) (17)
式中，K vol (θ s ,θ v ,φ s-v )和 K geo (θ s ,θ v ,φ s-v )表示核函数，是太阳照射角与传感器观测角度的函数， f iso (λ)、f vol (λ) 和f geo (λ)表示BRDF模型各向同性、体散射与几何光学核函数的权重系数；
以坡元实际太阳照射角代替太阳天顶角，以坡元传感器实际观测角代替传感器天顶角，得到坡元在传感器过境时刻的反射率：
Ρ T (λ) (i s , i v ,Φ s-v )=f iso (λ)+ f vol (λ)K vol (i s , i v ,Φ s-v )+ f geo (λ)K geo (i s , i v ,Φ s-v ) (18)
其中，Φ s-v =∣Φ s -Φ v ∣，
式中：Φ s-v 为太阳与传感器相对方位角，Φ s 表示由于地形影响坡元上实际来自太阳的照射角，Φ v 表示由于地形影响坡元上实际来自传感器方向的观测角。
本发明地形标准化方法将坡度、坡向因子引入核函数，发展了山区地表BRDF模型，同时将MODIS大气可降水量PW、气溶胶AOD产品以及MODIS BRDF/albedo核系数产品引入遥感影像地形标准化模型。以DEM为基础，提取各地形因子；借助及地表坡度、坡向数据，利用Ross Thick-Li Sparse核驱动模型，获得坡面BRDF归一化方向反射因子；同时将MODIS水汽及气溶胶产品（0.55um处气溶胶光学厚度和波长指数）作为大气输入参数，在山区辐射传输模型基础上，同时消除大气、地形对遥感影像的影响。由于模型不仅考虑了地形引起的地表入射方向上太阳辐射能量的变化，同时也考虑了由于地表二向反射特性，在获得坡面BRDF归一化方向反射因子基础上反演地表反射率。实验表明，该模型不仅去除了大气光学特性对遥感影像的辐射畸变，同时有效去除地形遮蔽等因子、地表BRDF特性引起的辐射畸变，从而获得地表真实反射率。以黑河流域上游大野口流域为研究区，将DEM数据与MODIS核系数产品为模型输入参数，对Landsat TM影像进行了地形标准化校正，获得地表真实反射率。
附图说明
图1是大气顶传感器接收到的辐亮度。
图2是遥感影像地形标准化流程图。
图3是天空可视因子图。
图4是Landsat过境时刻考虑地形的方向反射率图。
图5是Landsat过境时刻不考虑地形的方向反射率图。
图6是坡地BRDF归一化方向反射因子图。
图7是研究地区遥感影像的原始影像图。
图8是用伯朗体假设对图7所示原始影像标准化后的影像图。
图9是用本发明方法对图7所示原始影像标准化后的影像图。
图10是太阳照射角余弦值用遥感影像大气校正算法进行大气校正后的地物反射率在TM 2、TM3和TM 4波段的散点图。
图11是太阳照射角余弦值用本发明标准化方法进行地形标准化校正后的地物反射率在TM 2、TM3和TM 4波段的散点图。
图12是水域像元TM影像进行MBAC大气校正的平坦区地物波谱曲线与FLAASH和Landsat CDR的比较图。
图13是农田像元TM影像进行MBAC大气校正的平坦区地物波谱曲线与FLAASH和Landsat CDR的比较图。
图14是林地像元（西北坡上）TM影像进行MBBA地形标准化的地物波谱曲线与MBLA和Landsat CDR的比较图。
图15是草地像元（西南坡上）TM影像进行MBBA地形标准化的地物波谱曲线与MBLA和Landsat CDR的比较图。
具体实施方式
下面结合附图和具体实施方式对本发明进行详细说明。
一、模型方法
1、遥感影像地形标准化原理
如图1所示，卫星传感器获得的大气顶（TOA）辐射亮度L TOP (λ)由两部分组成：路径辐射L p (λ)和来自地表反射的辐射L(λ)，即：
L TOP (λ)= L p (λ)+ L(λ)T(λ,θ ν ) (1)
式中：T(λ,θ ν )表示传感器观测方向的大气透过率，是传感器观测天顶角θ ν 和波长λ的函数。路径辐射主要由瑞利散射和气溶胶散射决定，具体计算可参考文献（Li et al.,2002）。而来自地表反射的辐射能量L(λ)由地表-传感器方向反射率与地表接收的太阳辐射决定。在山区，太阳辐射由四部分组成（Zhang et al., 2015a）。因此，来自地表反射辐射可以写为：
L(λ)=(1/π)(E dir (λ)+ E aniso-dif (λ))×ρ T (λ)(i s ,i ν ,Φ s-ν )+ (1/π) (E iso-dif (λ)+E ref (λ))ρ Tsd (λ) (i ν ,Φ ν ) (2)
式中：E dir (λ)、E iso-dif (λ)、E aniso-dif (λ)和E ref (λ)分别表示太阳直接辐射、各向同性散射辐射、各向异性散射辐射和来自周围地形的反射辐射。ρ T (λ)和ρ Tsd (λ)分别表示坡地方向-方向反射率和半球-方向反射率。在一个坡度、坡向（S,A）的像素上，来自太阳和传感器方向的太阳天顶角和方位角分别为（θ s ,φ s ）与（θ v ,φ v ）。然而，由于地形影响，坡元上实际来自太阳的照射角和传感器方向的观测角变为（i s ,Φ s ）和（i v ,Φ v ）。记Φ s-v 为太阳与传感器相对方位角，即：Φ s-v =∣Φ s -Φ v ∣。因此，坡地反射率成为新坐标系统下四个角度的函数（Li et al., 2012）。考虑到传感器观测角在一景Landsat影像中变化较小，可以认为影像中每个像元所对应的观测天顶角和方位角为零，因此传感器实际观测方向的天顶角和方位角可简化为：
cosi v =cosS (3)
tanΦ v =0 (4)
对于相同地表覆盖类型而言，坡地与消除地形的平坦地表方向-方向反射率不同，其关系为：
ρ T (λ)( i s , i v ,Φ s-v )=Ω(λ) ( i s , i v ,Φ s-v ,θ s ,θ v ,φ s-v ) ×ρ H (λ) (θ s ,θ v ,φ s-v ) (5)
式中，Ω(λ) ( i s , i v ,Φ s-v ,θ s ,θ v ,φ s-v )为坡地BRDF归一化方向反射因子（Wenet al., 2008; flood et al., 2013a）。由于半球-方向反射率比较复杂，同时坡地各向同性散射辐射和来自周围地形反射辐射的能量较小，因此不区分坡地和平坦地表半球-方向反射率，即：
ρ Tsd (λ)(i v ,Φ v )=ρ H (λ) (θ s ,θ v ,φ s-v ) (6)
因此，将公式（2）、（5）和（6）代入公式（1），地形标准化后的反射率可以写为：
ρ H (λ) (θ s ,θ v ,φ s-v )=π(L TOP (λ)- L p (λ))/T(λ,θ v )〔(E dir (λ)+ E aniso-dif (λ)) Ω(λ)+E iso-dif (λ) + E ref (λ)〕 (7)
当Ω(λ)=1时，地表假设为朗伯体反射特征。存在两种情况，一种是对遥感影像只进行大气校正，公式可变为：
ρ H (λ) (θ s ,θ v ,φ s-v )=π(L TOP (λ)- L p (λ))/T(λ,θ v )E 0 (λ)D 0 T(λ,θ s ) (8)
另外一种情况是对遥感影像同时进行大气校正与地形校正，即地形标准化，但只是去除地形对太阳照射角的影响，则地表反射率反演公式为：
ρ H (λ) (θ s ,θ v ,φ s-v )=π(L TOP (λ)- L p (λ))/T(λ,θ v )E(λ) (9)
从以上分析来看，遥感影像地形标准化由以下三类输入参数组成：1）在太阳-地表-传感器路径上，太阳辐射二次经过大气的衰减作用；2）山区地表接收的太阳波谱辐照度；3）坡地BRDF归一化方向反射因子。其流程图如图2所示。
2、几种地形标准化方法比较
为了验证本方法的优越性，增加其它四种算法，如表1所示。FLAASH是基于辐射传输方程MODTRAN的大气校正模型，其大气水汽与气溶胶参数直接由影像自身获得。LandsatCDR气候数据集提供了进行大气校正的地表反射率数据，由美国地质勘探局USGS发布（http://earthexplorer.usgs.gov/）。以MODIS大气产品作为输入，利用公式（8）进行遥感影像大气校正算法称为MBAC算法。基于公式（9），在朗伯体假设条件下进行地形标准化，且以MODIS大气产品作为输入，称为MBLA算法。而本发明提出的地形标准化方法用MBBA表示。
表1 四种大气与地形标准化模型.
二、模型输入参数描述
1、大气透过率
大气透过率受气溶胶、水汽含量的影响较大，同时也受瑞利散射和臭氧吸收的影响。Landsat TM影像大气波谱透过率计算公式参照李新估算方案（Li et al., 2002），即：
T(λ,θ s )= T o3 (λ,θ s )T w (λ,θ s ) T r (λ,θ s )T a (λ,θ s ) (10)
式中：T(λ,θ s )、T o3 (λ,θ s )、T w (λ,θ s )、T r (λ,θ s )和T a (λ,θ s )分别表示总的大气波谱透过率、臭氧透过率、水汽透过率、瑞利散射透过率和气溶胶透过率。由于臭氧吸收与瑞利散射随时空分布变化不大，臭氧光学厚度可采用多年平均的空间分布数据，瑞利散射可利用经验公式估算。然而在山区，气溶胶和水汽分布随时空变化较为强烈，若利用以往半经验公式在单点获取一个值代替整个研究区将会带来问题。随着美国NASA对地观测系统EOS的快速发展，许多卫星传感器携带有较可靠的大气产品，如MODIS气溶胶产品AOD与水汽产品PW。水汽和气溶胶波谱透过率公式如下（Li et al., 2002）：
式中，k(λ)、m和w分别表示水汽有效吸收系数、大气质量及大气可降水厚度；α和β是angstrom波长指数与浑浊度。对于Landsat TM各波段水汽有效吸收系数及大气质量计算公式参加文献（Li et al., 2002）。但对于大气可降水厚度，在我们的算法中，MODIS水汽产品PW可直接作为水汽透过率输入参数。AOD产品中的波长指数α也可作为气溶胶透过率的输入参数。但是浑浊度β不能直接输入，需要进行推导，具体推导方法可参看文献（Zhang etal., 2015a）。由于MODIS AOD产品提供了0.55微米处的气溶胶光学厚度τ 0.55，因此β可以表示为：
Β=τ 0.55(0.55) a (14)
2、太阳谱辐照度估算
在山区，总的太阳波谱辐射由太阳直接辐射、太阳散射和周围地形反射辐射三部分组成，均为波长的函数。考虑到太阳方向对散射辐射的影响，将其又分解为各向同性散射辐射和各向异性散射辐射，即：
E(λ)= E dir (λ)+E iso-dif (λ)+E aniso-dif (λ)+E ref (λ) (15)
对于Landsat TM每个波段的太阳直接辐射与散射辐射同时受大气与地形作用的影响，其详细的估算方法参见文献Zhang et al.（2015a），只是Landsat TM大气顶太阳波谱辐照度常数见其官方网站（http://landsat.usgs.gov/science_L5_cpf.php）。
周围地形反射辐射估算方法如下：
式中，F ij 为地形结构因子，表示来自周围可见像元的反射辐射能量与可见像元自身总的太阳波谱辐射，具体计算参见（Li et al., 2002）。L i (λ)是周围可见像元的辐射亮度，由Landsat TM影像测量而得。
3、BRDF 模型
坡地BRDF 归一化方向反射因子是遥感影像地形标准化的重要参数，是坡地与消除地形的平坦地表方向-方向反射率之比。而这两种方向反射率均遵循地表BRDF规律。以EOS-MODIS 采用的Ross Thick-Li Sparse算法（AMBRALS算法）为代表的半经验线性核驱动模型由一定物理意义核的线性组合来拟合地表二向反射特性，备受人们的关注，其公式如下：
ρ H (λ) (θ s ,θ v ,φ s-v )=f iso (λ)+ f vol (λ)K vol (θ s ,θ v ,φ s-v )+ f geo (λ)K geo (θ s ,θ v ,φ s-v ) (17)
式中，K vol (θ s ,θ v ,φ s-v )和 K geo (θ s ,θ v ,φ s-v )表示核函数，是太阳照射角与传感器观测角度的函数，其详细的计算参见文献Schaaf et al. (1994)。f iso (λ)、f vol (λ) 和f geo (λ)表示BRDF模型各向同性、体散射与几何光学核函数的权重系数，直接通过MODIS BRDF核系数产品MCD43A1作为输入。
地表BRDF是波长及4个角度的函数。在平坦地表某瞬间，太阳-地表-传感器几何位置关系是固定的，太阳天顶角与方位角、传感器观测天顶角和方位角恒定。但是在山区，太阳-地表-传感器三者几何相对位置关系随地形而改变，甚至每个像元各不同。由于坡度、坡向的影响，坡元太阳相对天顶角与方位角、坡元相对传感器观测天顶角和方位角四个角度均发生了变化，具体公式参加文献Li et al.（2012）。
MODIS BRDF/Albedo 没有直接地考虑地形坡度、坡向的影响，不能直接应用于计算坡地BRDF，需要引入DEM数据，利用前向模型对其进行修改， 获得考虑坡度、坡向的天顶处实际方向反射率。在Ross Thick-Li Sparse 线性核驱动模型中以坡元实际太阳照射角代替太阳天顶角，以坡元传感器实际观测角代替传感器天顶角，得到坡元在传感器过境时刻的反射率：
Ρ T (λ) (i s , i v ,Φ s-v )=f iso (λ)+ f vol (λ)K vol (i s , i v ,Φ s-v )+ f geo (λ)K geo (i s , i v ,Φ s-v ) (18)
其中，Φ s-v =∣Φ s -Φ v ∣。
三、模型精度验证
下面分别从目视效果与定量评价两方面检验本发明提供的地形标准化模型的精度：
在研究区，Landsat TM在当地过境的同一天与MODIS 上午星过境时刻一般相差1小时左右，因此MODIS大气产品成为TM影像地形标准化的大气输入参数。图3～图6是研究区天空可视因子及坡地BRDF归一化方向反射因子的空间分布图。图3是天空可视因子图，图4是不考虑地形的方向反射率图，图5是考虑地形的方向反射率图，图6是坡地BRDF归一化方向反射因子图。为了目视评价耦合MODIS大气产品与地表BRDF特性的地形标准化模型去除地形效应的优势，将研究区地形标准化结果部分放大显示，即选取图7所示的遥感影像原始影像，然后用朗伯体假设对图7的影像进行标准化，产生图8所示的影像图；用本发明方法对对图7的影像进行标准化，产生图9所示的影像图。从图7可以看出，由于地形崎岖使得原始影像上面向太阳方向的阳坡地物较亮，而背向太阳方向的阴坡地物较暗，从而产生地形起伏与阴暗效应明显，为遥感影像解译与定量遥感研究带来困难。如图8和图9所示，经过地形标准化后的遥感影像，地形平坦化，减弱了阴坡与阳坡辐射亮度的差异，同时也提高了阴坡与半阴坡处地物的反射率。去除地形效应后的反射率影像上，阴坡处更多地物细节信息凸现。然而，由于地表朗伯体反射特性的假设，往往在阴影处出现过校正问题，如图8中A、B所示地物。而本发明标准化方法考虑了地表BRDF反射特性，抑制了这种过校正效应，如图9中A、B所示地物所示。
由于一方面地表波谱反射率测量值缺乏，另一方面在西北干旱区，地形特征决定了阴坡阳坡分布不同地表覆盖类型地物，选择阴坡、阳坡同种地表覆盖类型验证地形标准化结果比较困难。因此，像元太阳照射角余弦值与地表反射率是否存在线性相关成为一种评判地形标准化效果的常用定量检测方法。在不同地形条件与地表覆盖类型中任选481个像元，制作太阳照射角余弦值与只做大气校正的MBAC及地形标准化MBBA后的地物分别在TM2、TM3和TM 4波段的反射率散点图，如图10和图11所示。未消除地形影响的地表反射率随地形坡度、坡向而发生变化，随太阳相对入射角余弦值的增大而增大，图10表明只进行大气校正的地表反射率随地形坡度、坡向而发生变化，随太阳相对入射角余弦值的增大而增大。然而，对遥感影像进行地形标准化处理后，由于去除了地形对遥感影像的影响，其反演的反射率与太阳相对入射角余弦值之间基本消除了这种线性依赖关系，从而恢复地表真实反射率（图11）。
另外，提取四种不同覆盖类型地物反射波谱曲线用于进一步比较本方法地形标准化优势的评价。首先，在平坦地表选择两种地物，分别比较其在FLAASH与Landsat CDR数据集的波谱反射率，用于分析本发明标准化方法进行大气校正（MBAC）的有效性。图12和图13所示的地物波谱曲线变化表明，基于MODIS水汽和气溶胶产品的大气校正模型能精确描述MODIS过境时刻大气状况，即MBAC模型能够有效去除水汽、气溶胶对遥感影像的影响，从而恢复了水体与农作物真实波谱反射率。例如，水体反射率值大气校正前比校正后较高，尤其是在近红外波段。这是因为大气校正去除了程辐射对遥感影像的影响，从而恢复了水体真实反射率，使得近红外波段几乎为零。同时，在坡度为34°的西北坡选择一个林地像元，在坡度为15°的西南坡选择一个草地像元。本发明的地形标准化方法在MODIS大气产品及地形效应基础上，构建了山区地表BRDF反射率模型。从图14和图15可以看出，与MBLA模型和Landsat CDR数据集相比，利用本方法地形标准化后，阴坡林地与半阴坡草地反射率值得到提高，恢复了阴影处地物的信息。
从以上对比分析来看，本标准化方法提供的遥感影像地形标准化模型校正效果较好，同步进行大气校正与地形校正后反演的地物波谱反射率值更真实客观。
本发明标准化方法从三个方面进行重点研究：一是在Ross Thick-Li Sparse线性核驱动模型中，将坡度、坡向引入Ross体散射核函数与Li氏几何光学核函数，建立山区BRDF模型；二是在现有遥感影像地形标准化成果（Li et al., 1999, 2002）基础上，改变地表朗伯反射特征假设，引入山区BRDF模型，最大限度减少复杂地形区BRDF 形状的改变对遥感数据的影响；三是，将MODIS水汽和气溶胶大气产品作为地形标准化模型参数，获得卫星影像过境时刻的大气参数，消除大气光学因素的影响。
Claims (1)
1.一种遥感影像地形标准化方法，其特征在于，该方法具体为：
1）卫星传感器获得的大气顶辐射亮度L TOP (λ)由路径辐射L p (λ)和来自地表反射的辐射L(λ)，即：
L TOP (λ)= L p (λ)+ L(λ)T(λ,θ ν ) (1)
在山区，来自地表反射辐射为：
L(λ)=(1/π)(E dir (λ)+ E aniso-dif (λ))×ρ T (λ)(i s ,i ν ,Φ s-ν )+ (1/π) (E iso-dif (λ)+ E ref (λ))ρ Tsd (λ) (i ν ,Φ ν ) (2)
式中：E dir (λ)、E iso-dif (λ)、E aniso-dif (λ)和E ref (λ)分别表示太阳直接辐射、各向同性散射辐射、各向异性散射辐射和来自周围地形的反射辐射；ρ T (λ)和ρ Tsd (λ)分别表示坡地方向-方向反射率和半球-方向反射率；
对于相同地表覆盖类型而言，坡地与消除地形的平坦地表方向-方向反射率不同，其关系为：
ρ T (λ)( i s , i v ,Φ s-v )=Ω(λ) ( i s , i v ,Φ s-v ,θ s ,θ v ,φ s-v ) ×ρ H (λ) (θ s ,θ v ,φ s-v ) (5)
式中，Ω(λ) ( i s , i v ,Φ s-v ,θ s ,θ v ,φ s-v )为坡地BRDF归一化方向反射因子；由于半球-方向反射率比较复杂，同时坡地各向同性散射辐射和来自周围地形反射辐射的能量较小，因此不区分坡地和平坦地表半球-方向反射率，即：
ρ Tsd (λ)(i v ,Φ v )=ρ H (λ) (θ s ,θ v ,φ s-v ) (6)
因此，将公式（2）、公式（5）和公式（6）代入公式（1），地形标准化后的反射率为：
ρ H (λ) (θ s ,θ v ,φ s-v )=π(L TOP (λ)- L p (λ))/T(λ,θ v )〔(E dir (λ)+ E aniso-dif (λ)) Ω(λ)+ E iso-dif (λ) + E ref (λ)〕 (7)
当Ω(λ)=1时，地表假设为朗伯体反射特征；存在两种情况，一种是对遥感影像只进行大气校正，公式可变为：
ρ H (λ) (θ s ,θ v ,φ s-v )=π(L TOP (λ)- L p (λ))/T(λ,θ v )E 0 (λ)D 0 T(λ,θ s ) (8)
另外一种情况是对遥感影像同时进行大气校正与地形校正，即地形标准化，但只是去除地形对太阳照射角的影响，则地表反射率反演公式为：
ρ H (λ) (θ s ,θ v ,φ s-v )=π(L TOP (λ)- L p (λ))/T(λ,θ v )E(λ) (9)
2）大气透过率受气溶胶、水汽含量的影响较大，同时也受瑞利散射和臭氧吸收的影响，即：
T(λ,θ s )= T o3 (λ,θ s )T w (λ,θ s ) T r (λ,θ s )T a (λ,θ s ) (10)
式中：T(λ,θ s )、T o3 (λ,θ s )、T w (λ,θ s )、T r (λ,θ s )和T a (λ,θ s )分别表示总的大气波谱透过率、臭氧透过率、水汽透过率、瑞利散射透过率和气溶胶透过率， 由于臭氧吸收与瑞利散射随时空分布变化不大，臭氧光学厚度可采用多年平均的空间分布数据，瑞利散射可利用经验公式估算；然而在山区，气溶胶和水汽分布随时空变化较为强烈，水汽和气溶胶波谱透过率公式如下：
式中，k(λ)、m和w分别表示水汽有效吸收系数、大气质量及大气可降水厚度；α和β是angstrom波长指数与浑浊度；对于大气可降水厚度，MODIS水汽产品PW直接作为水汽透过率输入参数；AOD产品中的波长指数α也作为气溶胶透过率的输入参数；β表示为：
Β=τ 0.55(0.55) a (14)
在山区，总的太阳波谱辐射由太阳直接辐射、太阳散射和周围地形反射辐射三部分组成，均为波长的函数；考虑到太阳方向对散射辐射的影响，将其又分解为各向同性散射辐射和各向异性散射辐射，即：
E(λ)= E dir (λ)+E iso-dif (λ)+E aniso-dif (λ)+E ref (λ) (15)
对于Landsat TM每个波段的太阳直接辐射与散射辐射同时受大气与地形作用的影响；
周围地形反射辐射估算方法如下：
式中，F ij 为地形结构因子，表示来自周围可见像元的反射辐射能量与可见像元自身总的太阳波谱辐射； L i (λ)是周围可见像元的辐射亮度，由Landsat TM影像测量而得；
3）坡地BRDF 归一化方向反射因子是遥感影像地形标准化的重要参数，是坡地与消除地形的平坦地表方向-方向反射率之比，而这两种方向反射率均遵循地表BRDF规律；半经验线性核驱动模型由一定物理意义核的线性组合来拟合地表二向反射特性，其公式如下：
ρ H (λ) (θ s ,θ v ,φ s-v )=f iso (λ)+ f vol (λ)K vol (θ s ,θ v ,φ s-v )+ f geo (λ)K geo (θ s ,θ v ,φ s-v ) (17)
式中，K vol (θ s ,θ v ,φ s-v )和K geo (θ s ,θ v ,φ s-v )表示核函数，是太阳照射角与传感器观测角度的函数， f iso (λ)、f vol (λ) 和f geo (λ)表示BRDF模型各向同性、体散射与几何光学核函数的权重系数；
以坡元实际太阳照射角代替太阳天顶角，以坡元传感器实际观测角代替传感器天顶角，得到坡元在传感器过境时刻的反射率：
Ρ T (λ) (i s , i v ,Φ s-v )=f iso (λ)+ f vol (λ)K vol (i s , i v ,Φ s-v )+ f geo (λ)K geo (i s , i v ,Φ s-v ) (18)
其中，Φ s-v =∣Φ s -Φ v ∣，
式中：Φ s-v 为太阳与传感器相对方位角，Φ s 表示由于地形影响坡元上实际来自太阳的照射角，Φ v 表示由于地形影响坡元上实际来自传感器方向的观测角。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510211271.6A CN104834814B (zh) | 2015-04-29 | 2015-04-29 | 遥感影像地形标准化方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201510211271.6A CN104834814B (zh) | 2015-04-29 | 2015-04-29 | 遥感影像地形标准化方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN104834814A CN104834814A (zh) | 2015-08-12 |
| CN104834814B true CN104834814B (zh) | 2018-04-10 |
Family
ID=53812697
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201510211271.6A Active CN104834814B (zh) | 2015-04-29 | 2015-04-29 | 遥感影像地形标准化方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN104834814B (zh) |
Families Citing this family (9)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN105975777B (zh) * | 2016-05-04 | 2021-01-26 | 中国科学院合肥物质科学研究院 | 顾及实际天空光分布影响的地表反照率遥感模型 |
| CN107247038B (zh) * | 2017-06-14 | 2020-02-14 | 电子科技大学 | 一种获取河流冰凌红外波段散射特性的方法 |
| CN107656289A (zh) * | 2017-08-23 | 2018-02-02 | 中国科学院光电研究院 | 基于地基辐亮度的星载光学载荷绝对辐射定标方法及系统 |
| CN108132220B (zh) * | 2017-12-25 | 2021-03-05 | 中国林业科学研究院资源信息研究所 | 林区机载推扫式高光谱影像的brdf归一化校正方法 |
| CN109934788B (zh) * | 2019-03-22 | 2022-12-30 | 鲁东大学 | 一种基于标准遥感图像的遥感图像缺失数据修复方法 |
| CN111198162B (zh) * | 2020-01-09 | 2022-10-18 | 胡德勇 | 一种城区表面反射率遥感反演方法 |
| CN111737958B (zh) * | 2020-06-05 | 2023-06-30 | 中国科学院空天信息创新研究院 | 一种遥感模型数据标准化处理方法 |
| CN113553551B (zh) * | 2021-07-28 | 2024-02-02 | 生态环境部华南环境科学研究所 | 一种耦合景观格局的臭氧浓度预测系统 |
| CN115524763B (zh) * | 2022-09-27 | 2023-08-11 | 中国科学院、水利部成都山地灾害与环境研究所 | 一种多时相高分辨率山地卫星影像地形辐射校正方法 |
Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102622738A (zh) * | 2012-03-08 | 2012-08-01 | 北京师范大学 | 一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法 |
| CN103760565A (zh) * | 2014-02-10 | 2014-04-30 | 中国科学院南京地理与湖泊研究所 | 一种区域尺度森林冠层高度遥感反演方法 |
| CN103940407A (zh) * | 2014-02-13 | 2014-07-23 | 鲁东大学 | 一种基于地形和遥感影像融合技术提取冲沟方法 |
| CN104166782A (zh) * | 2014-06-05 | 2014-11-26 | 刘健 | 一种林地土壤有机碳遥感估测的方法 |
-
2015
- 2015-04-29 CN CN201510211271.6A patent/CN104834814B/zh active Active
Patent Citations (4)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102622738A (zh) * | 2012-03-08 | 2012-08-01 | 北京师范大学 | 一种Landsat TM/ETM+图像中山体阴影区的光谱信息恢复方法 |
| CN103760565A (zh) * | 2014-02-10 | 2014-04-30 | 中国科学院南京地理与湖泊研究所 | 一种区域尺度森林冠层高度遥感反演方法 |
| CN103940407A (zh) * | 2014-02-13 | 2014-07-23 | 鲁东大学 | 一种基于地形和遥感影像融合技术提取冲沟方法 |
| CN104166782A (zh) * | 2014-06-05 | 2014-11-26 | 刘健 | 一种林地土壤有机碳遥感估测的方法 |
Non-Patent Citations (7)
| Title |
|---|
| an integrated approach to estimate shortwave solar radiation on clear-sky days in rugged terrain using MODIS atmospheric products;yanli Zhang等;《solar energy》;20150206;第2.1节、第2.2节 * |
| an operational scheme for deriving standardised surface reflectance from landsat TM/ETM+ and spot hrg imagery for eastern australia;neil flood等;《remote sensing》;20130104;第2.4节 * |
| Modeling the land surface relectance for optical remote sensing data in rugged terrain;Wen Jianguang等;《science in china series D:earth sciences》;20080815;第51卷(第8期);第1169-1178页 * |
| retrieval of snow reflectance from landsat data in rugged terrain;Li Xin等;《annals of glaciology》;20020101;第34卷(第1期);第31-37页 * |
| sensitivity of topographic correction to the DEM spatial scale;Yanli Zhang等;《IEEE geoscience and remote sensing letters》;20140630;第12卷(第1期);图1，第2节 * |
| topographic effects on bidirectional and hemispherical relectances calculated with a geometric-optical conopy model;Crystal Barker等;《IEEE transactions on geoscience and remote sensing》;20020806;第32卷(第6期);第1186-1193页 * |
| 利用遥感资料估算复杂地形条件下的净辐射;李净;《万方在线公开：http://d.wanfangdata.com.cn/Thesis/y1630081》;20110525;第3章 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN104834814A (zh) | 2015-08-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104834814B (zh) | 遥感影像地形标准化方法 | |
| Mondejar et al. | Near infrared band of Landsat 8 as water index: a case study around Cordova and Lapu-Lapu City, Cebu, Philippines | |
| Fisher et al. | Comparing Landsat water index methods for automated water classification in eastern Australia | |
| Harmel et al. | Sunglint correction of the Multi-Spectral Instrument (MSI)-SENTINEL-2 imagery over inland and sea waters from SWIR bands | |
| Snyder et al. | BRDF models to predict spectral reflectance and emissivity in the thermal infrared | |
| Li et al. | A physics-based atmospheric and BRDF correction for Landsat data over mountainous terrain | |
| Li et al. | Satellite-derived land surface temperature: Current status and perspectives | |
| Sterckx et al. | SIMilarity Environment Correction (SIMEC) applied to MERIS data over inland and coastal waters | |
| Bernstein et al. | Validation of the QUick atmospheric correction (QUAC) algorithm for VNIR-SWIR multi-and hyperspectral imagery | |
| Walther et al. | Implementation of the daytime cloud optical and microphysical properties algorithm (DCOMP) in PATMOS-x | |
| Maignan et al. | Polarized reflectances of natural surfaces: Spaceborne measurements and analytical modeling | |
| CN104406686B (zh) | 复杂地形条件下太阳短波入射辐射估算方法 | |
| Harmel et al. | Estimation of the sunglint radiance field from optical satellite imagery over open ocean: Multidirectional approach and polarization aspects | |
| Salama et al. | Two-stream remote sensing model for water quality mapping: 2SeaColor | |
| Dierssen | Hyperspectral measurements, parameterizations, and atmospheric correction of whitecaps and foam from visible to shortwave infrared for ocean color remote sensing | |
| Wu et al. | The definition of remotely sensed reflectance quantities suitable for rugged terrain | |
| Mishra et al. | Assessment of different topographic corrections in AWiFS satellite imagery of Himalaya terrain | |
| Zhang et al. | Sensitivity of topographic correction to the DEM spatial scale | |
| Cahalan et al. | Cloud characterization and clear-sky correction from Landsat-7 | |
| Rotta et al. | Atmospheric correction assessment of SPOT-6 image and its influence on models to estimate water column transparency in tropical reservoir | |
| Wang et al. | An adaptive atmospheric correction algorithm for the effective adjacency effect correction of submeter-scale spatial resolution optical satellite images: Application to a WorldView-3 panchromatic image | |
| Wu et al. | A comparison of illumination geometry-based methods for topographic correction of QuickBird images of an undulant area | |
| Al-Shehhi | Uncertainty in satellite sea surface temperature with respect to air temperature, dust level, wind speed and solar position | |
| Meng et al. | An analytical model for optical polarimetric imaging systems | |
| Hadjimitsis et al. | The use of an improved atmospheric correction algorithm for removing atmospheric effects from remotely sensed images using an atmosphere–surface simulation and meteorological data |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| C06 | Publication | ||
| PB01 | Publication | ||
| EXSB | Decision made by sipo to initiate substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20190201 Address after: 310004 21 Floor of Blue International Building 533 Dongxin Road, Hangzhou City, Zhejiang Province Patentee after: HANGZHOU DADI TECHNOLOGY CO., LTD. Address before: 730070 No. 967 Ning Ning Road, Anning District, Lanzhou, Gansu Patentee before: Northwest Normal University |