CN104951656B

CN104951656B - 宽视域卫星影像地表反射率反演方法

Info

Publication number: CN104951656B
Application number: CN201510345432.0A
Authority: CN
Inventors: 张兆明; 何国金; 王猛猛; 龙腾飞; 王桂周; 张晓美
Original assignee: Institute of Remote Sensing and Digital Earth of CAS
Current assignee: Institute of Remote Sensing and Digital Earth of CAS
Priority date: 2015-06-23
Filing date: 2015-06-23
Publication date: 2018-01-12
Anticipated expiration: 2035-06-23
Also published as: CN104951656A

Abstract

太阳照射天顶角和卫星观测天顶角是卫星遥感影像大气校正和地表反射率反演中很重要的输入参数。针对视场角较大的影像，本发明逐像元计算太阳照射天顶角和卫星观测天顶角，然后反演地表反射率，提出一种宽视域陆地卫星遥感影像地表反射率获取的方法。

Description

宽视域卫星影像地表反射率反演方法

技术领域

本发明涉及一种宽视域卫星遥感影像地表反射率反演的方法，能够应用在林业、农业、气象、生态环境等行业部门。本发明提出的方法适用于星下点成像的陆地卫星影像。

背景技术

太阳照射方向和卫星观测方向(天顶角和方位角)信息对于遥感数据的处理和应用具有重要的作用。太阳照射天顶角和卫星观测天顶角的不同会导致辐射传输中大气路径的长短，从而影响太阳照射方向和卫星观测方向的大气透过率，其次，在给定的太阳入射条件下，传感器所接收到的反射辐射能量还受到卫星观测角度的影响；同样的太阳照射条件，卫星观测角度的变化会产生不同的地表反射率反演值，即BRDF(BidirectionalReflectance Distribution Function，BRDF，双向反射分布函数)效应。因此太阳照射天顶角和卫星观测天顶角是卫星遥感影像大气校正和地表反射率反演中很重要的输入参数。然而现有的陆地卫星遥感影像地表反射率获取方法中，整景影像通常只利用一个太阳照射天顶角和卫星观测天顶角值(景中心点的)。对于视场角较小的卫星影像(如ASTER影像)，一景影像范围内太阳照射角度和卫星观测角度的变化较小，整景影像只利用一组角度信息通常不会带来大的误差。对于视场角较大的影像(如GF-1宽幅影像)，则需要针对影像内不同的区域提供更多的角度信息。事实上，对于高纬地区，即使对于ASTER之类视场角较小的卫星影像，一景给定一组角度值也是不够的，需要精确获取逐像元的太阳照射天顶角和卫星观测天顶角。针对现有技术的不足，本发明针对宽视域陆地卫星遥感影像提出一种地表反射率获取的方法。

发明内容

本发明的目的在于提供一种宽视域卫星遥感影像地表反射率获取技术，以克服现有方法的不足，从而促进相关领域遥感应用的发展。

为实现上述目的，本发明提出的方法包括以下步骤：

第一步、对卫星影像进行几何校正；

第二步、辐射定标，将影像像元亮度值转换为星上辐射亮度L_λ；

L_λ＝Gain×DN+Bias

L_λ为波段λ的星上辐射亮度，DN为波段λ的像元值，Gain和Bias分别为波段λ的增益和偏置。

第三步、获取大气程辐射L_p

L_p＝Gain×DN_min+Bias-0.01E₀cosθ_zT_zT_v/πd²

第四步、逐像元计算太阳天顶角

(4-1)计算太阳时：

其中t是太阳时(小时数，带小数位)，t_s是标准时(小时数，带小数位)，SM是该时区标准经线的经度(弧度)，L是该像元点的经度(弧度)，J是儒略日(即影像获取日期距离1月1日的天数)。

(4-2)计算太阳赤纬：

其中δ是太阳赤纬(弧度)，J是儒略日。

(4-3)计算太阳天顶角

其中θ_z是太阳天顶角(弧度)，l是该像元点的纬度(弧度)，δ是太阳赤纬(弧度)，t是太阳时。

第五步、在陆地卫星影像星下点成像情况下，按照以下步骤逐像元计算卫星观测天顶角：

(5-1)获取任一像元点与影像中心点的行距和列距；

(5-2)由行距和列距计算像元点与影像中心点的距离；

(5-3)由像元点与影像中心点的距离和卫星高度之间的三角关系获取该像元的卫星观测天顶角。

第六步、计算大气光学厚度τ

(6-1)根据陆地卫星遥感影像的经纬度和影像成像时间，获取同一地理范围、同步过境的MODIS气溶胶光学厚度(0.47μm和0.66μm)、总可降水汽、地表气压和臭氧浓度影像，并将各影像做几何校正，转换为和陆地卫星遥感影像一致的投影方式和分辨率，且影像行列数一致；

(6-2)由MODIS总可降水汽w计算大气水蒸汽光学厚度；

其中，τ_w是大气水蒸汽光学厚度，a_wλ是水汽吸收系数，w是总可降水汽，M是相对大气量；

(6-3)计算臭氧吸收光学厚度

τ_o＝C_ozone×A_ozone(λ)

τ_o是臭氧吸收光学厚度，C_ozone为MODIS臭氧浓度(单位为Dobson)，A_ozone(λ)为波段λ的臭氧吸收系数；

(6-4)由MODIS 0.47μm和0.66μm的气溶胶光学厚度，利用下式计算波长指数(α)及大气浑浊度系数(β)

β＝τ_a(λ₁)λ₁ ^α＝τ_α(λ₂)λ₂ ^α

λ₁和λ₂分别为0.47μm和0.66μm，τ_a(λ₁)和τ_a(λ₂)分别为0.47μm和0.66μm的气溶胶光学厚度；

在得到α和β后，可以利用下式得到任意波长λ的气溶胶光学厚度τ_a(λ)：

τ_a(λ)＝βλ^-α

其中，λ为波长(μm)，τ_a(λ)是该波长的气溶胶光学厚度；

(6-5)计算瑞利散射光学厚度；

其中，τ_r为瑞利散射光学厚度，λ是影像各波段的中心波长(μm)；p为MODIS的地表气压产品；

(6-6)τ＝τ_w+τ_o+τ_a+τ_τ

第七步、计算太阳照射方向和卫星观测方向的大气透过率；

T_z＝exp(-τ/cosθ_z)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_z}

T_v＝exp(-τ/cosθ_v)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_v}

其中，T_z和T_v分别是太阳照射方向和卫星观测方向的大气透过率，τ是大气光学厚度，θ_z是太阳天顶角，θ_v是卫星观测天顶角，τ_r、τ_a、τ_o和τ_w分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光学厚度和大气水蒸汽光学厚度；

第八步、计算地表反射率ρ

d＝1+0.0167sin[2π(J-93.5)/365]

E₀是大气层外相应波长的太阳光谱辐照度，d是日地距离(天文单位)，J是儒略日。

附图说明

图1是逐像元计算卫星观测天顶角示意图

具体实施方式

遥感影像地表反射率反演是以辐射传输方程为基础的，在假定地表均匀、朗伯面反射、大气性质均一，忽略大气多次散射的前提下，星上辐射亮度和地表反射率的关系如式1所示：

(式1)

式1中，L_λ是星上辐射亮度，L_p是程辐射，F_d是地表接收到的辐照度，T_v是卫星观测方向的大气透过率，ρ是地表反射率，d是日地距离，s是大气下界的半球反射率，因为s值很小，通常可以忽略，所以由式1可以得到地表反射率的计算公式(式2)：

(式2)

在不考虑大气下行漫射辐射的情况下，F_d可近似等于E₀cosθ_zT_z，其中E₀是大气层外相应波长的太阳光谱辐照度，θ_z是太阳天顶角，T_z是太阳照射方向的大气透过率。

L_λ由传感器辐射定标得到

L_λ＝Gain×DN+Bias (式3)

程辐射L_p利用暗目标方法来确定，即假定遥感影像中存在暗目标，暗目标的地表反射率应为0，但是由于程辐射的影响，在遥感影像上暗目标显示一定的DN值(DN_min)，可以认为该部分辐射即为程辐射，利用暗目标的DN值来确定程辐射值。影像中的暗目标可以选择山体阴影、洁净的水体或者浓密植被区域。

DN_min值确定后，利用下式计算L_p：

L_p＝Gain×DN_min+Bias (式4)

理想情况下，暗目标的反射率为0，实际的暗目标并非“全黑”，而是具有一个非常小的反射率，可以假定该反射率为0.01，因此需要将该部分辐射从上式中减掉，修正后的程辐射计算公式为：

L_p＝Gain×DN_min+Bias-0.01E₀cosθ_zT_zT_v/πd² (式5)

T_v和T_z按式6-7获得：

T_v＝exp(-τ/cosθ_v)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_v} (式6)

T_z＝exp(-τ/cosθ_z)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_z} (式7)

其中，τ是大气光学厚度，τ＝τ_r+τ_a+τ_o+τ_w，τ_r、τ_a、τ_o和τ_w分别是瑞利散射光学厚度、气溶胶光学厚度、臭氧吸收光学厚度和大气水蒸汽光学厚度，θ_v是卫星观测天顶角，θ_z是太阳天顶角；

(式8)

由MODIS 0.47μm和0.66μm的气溶胶光学厚度，利用式9-10计算波长指数(α)及大气浑浊度系数(β)

(式9)

β＝τ_a(λ₁)λ₁ ^α＝τ_α(λ₂)λ₂ ^α (式10)

在得到α和β后，可以利用式11得到任意波长λ的气溶胶光学厚度τ_a(λ)：

τ_a(λ)＝βλ^-α (式11)

其中，λ为波长(μm)，τ_a(λ)是该波长的气溶胶光学厚度；

臭氧吸收光学厚度τ_o按式12获得：

τ_o＝C_ozone×A_ozone(λ) (式12)

τ_o是臭氧吸收光学厚度，C_ozone为MODIS臭氧浓度(单位为Dobson)，A_ozone(λ)为波段λ的臭氧吸收系数，可以根据文献Aerosol Optical Depth Value-Added Product(Koontz Aet al.2013，https：//www.arm.gov/publications/tech_reports/doe-sc-arm-tr-129.pdf)得到。

大气水蒸汽光学厚度τ_w由总可降水汽w计算：

(式13)

其中，τ_w是大气水蒸汽光学厚度，a_wλ是水汽吸收系数，各波长的a_wλ值可以查文献得到[Bird R E，Riordan C.Simple solar spectral model for direct and diffuseirradiance on horizontal and tilted planes at the earth′s surface forcloudless atmospheres.Journal of climate and applied meteorology，1986，25：87-97.]。w是MODIS总可降水汽(cm)，M是相对大气量，由式14获得：

M＝[cosθ_z+0.15(93.885-θ_z)^-1.253]^-1 (式14)

θ_z是太阳天顶角。

日地距离d按式15计算：

d＝1+0.0167sin[2π(J-93.5)/365] (式15)

其中，J是儒略日，即影像获取日期距离1月1日的天数；

按式16-18计算逐像元的太阳天顶角θ_z

(式16)

其中t是太阳时(小时数，带小数位)，t_s是标准时(小时数，带小数位)，SM是该时区标准经线的经度(弧度)，L是该像元点的经度(弧度)，J是儒略日。

(式17)

其中δ是太阳赤纬(弧度)，J是儒略日。

(式18)

在陆地卫星影像星下点成像情况下，按照以下步骤逐像元计算卫星观测天顶角θ_v

1.获取任一像元点E与影像中心点O的行距(Δx)和列距(Δy)；

2.计算点E与点O的距离s

3.由s和卫星高度h之间的三角关系获取该像元点的卫星观测天顶角θ_v

θ_v＝arctan(s/h) 。

Claims

1.一种宽视域卫星影像地表反射率反演的方法，包括如下步骤：

第一步、对卫星影像进行几何校正；

L_λ＝Gain×DN+Bias

L_λ为波段λ的星上辐射亮度，DN为波段λ的像元值，Gain和Bias分别为波段λ的增益和偏置；

第三步、获取大气程辐射L_p

L_p＝Gain×DN_min+Bias-0.01E₀cosθ_zT_zT_v/πd²

L_p是大气程辐射，DN_min是暗目标DN值，E₀是大气层外相应波长的太阳光谱辐照度，θ_z是太阳天顶角，T_z是太阳照射方向的大气透过率，T_v是卫星观测方向的大气透过率，d是日地距离；

第四步、逐像元计算太阳天顶角

(4-1)计算太阳时：

其中t是太阳时，t_s是标准时，SM是该时区标准经线的经度，L是该像元点的经度，J是儒略日；

(4-2)计算太阳赤纬：

<mrow> <mi>&delta;</mi> <mo>=</mo> <mn>0.4093</mn> <mi>sin</mi> <mrow> <mo>(</mo> <mfrac> <mrow> <mn>2</mn> <mi>&pi;</mi> <mrow> <mo>(</mo> <mi>J</mi> <mo>-</mo> <mn>81</mn> <mo>)</mo> </mrow> </mrow> <mn>368</mn> </mfrac> <mo>)</mo> </mrow> </mrow>

其中δ是太阳赤纬，J是儒略日；

(4-3)计算太阳天顶角

<mrow> <msub> <mi>&theta;</mi> <mi>z</mi> </msub> <mo>=</mo> <mfrac> <mi>&pi;</mi> <mn>2</mn> </mfrac> <mo>-</mo> <mi>arcsin</mi> <mrow> <mo>(</mo> <mi>sin</mi> <mi>l</mi> <mi>sin</mi> <mi>&delta;</mi> <mo>+</mo> <mi>cos</mi> <mi>l</mi> <mi>cos</mi> <mi>&delta;</mi> <mi>cos</mi> <mfrac> <mrow> <mi>&pi;</mi> <mi>t</mi> </mrow> <mn>12</mn> </mfrac> <mo>)</mo> </mrow> </mrow>

其中θ_z是太阳天顶角，l是该像元点的纬度，δ是太阳赤纬，t是太阳时；

(5-1)获取任一像元点E与影像中心点O的行距Δx和列距Δy；

(5-2)由行距和列距计算像元点E与影像中心点O的距离s

<mrow> <mi>s</mi> <mo>=</mo> <msqrt> <mrow> <msup> <mrow> <mo>(</mo> <mi>&Delta;</mi> <mi>x</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> <mo>+</mo> <msup> <mrow> <mo>(</mo> <mi>&Delta;</mi> <mi>y</mi> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mo>;</mo> </mrow>

(5-3)由s和卫星高度h之间的三角关系获取该像元点的卫星观测天顶角θ_v

θ_v＝arctan(s/h)；

第六步、计算大气光学厚度τ

(6-1)根据陆地卫星遥感影像的经纬度和影像成像时间，获取同一地理范围、同步过境的MODIS 0.47μm和0.66μm的气溶胶光学厚度、总可降水汽、地表气压和臭氧浓度影像，并将各影像做几何校正，转换为和陆地卫星遥感影像一致的投影方式和分辨率，且影像行列数一致；

(6-2)由MODIS总可降水汽w计算大气水蒸汽光学厚度；

<mrow> <msub> <mi>&tau;</mi> <mi>w</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mn>0.2385</mn> <msub> <mi>a</mi> <mrow> <mi>w</mi> <mi>&lambda;</mi> </mrow> </msub> <mi>w</mi> <mi>M</mi> </mrow> <msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <mn>20.07</mn> <msub> <mi>a</mi> <mrow> <mi>w</mi> <mi>&lambda;</mi> </mrow> </msub> <mi>w</mi> <mi>M</mi> <mo>)</mo> </mrow> <mn>0.45</mn> </msup> </mfrac> </mrow>

(6-3)计算臭氧吸收光学厚度

τ_o＝C_ozone×A_ozone(λ)

τ_o是臭氧吸收光学厚度，C_ozone为MODIS臭氧浓度，A_ozone(λ)为波段λ的臭氧吸收系数；

(6-4)由MODIS 0.47μm和0.66μm的气溶胶光学厚度，利用下式计算波长指数α及大气浑浊度系数β

<mrow> <mi>&alpha;</mi> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <msub> <mi>ln&tau;</mi> <mi>a</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&lambda;</mi> <mn>1</mn> </msub> <mo>)</mo> </mrow> <mo>-</mo> <msub> <mi>ln&tau;</mi> <mi>a</mi> </msub> <mrow> <mo>(</mo> <msub> <mi>&lambda;</mi> <mn>2</mn> </msub> <mo>)</mo> </mrow> </mrow> <mrow> <msub> <mi>ln&lambda;</mi> <mn>1</mn> </msub> <mo>-</mo> <msub> <mi>ln&lambda;</mi> <mn>2</mn> </msub> </mrow> </mfrac> </mrow>

β＝τ_a(λ₁)λ₁ ^α＝τ_a(λ₂)λ₂ ^α

τ_a为气溶胶光学厚度，λ₁和λ₂分别为0.47μm和0.66μm，τ_a(λ₁)和τ_a(λ₂)分别为0.47μm和0.66μm的气溶胶光学厚度；

τ_a(λ)＝βλ^-α

其中，λ为波长，τ_a(λ)是波长λ的气溶胶光学厚度；

(6-5)计算瑞利散射光学厚度；

<mrow> <msub> <mi>&tau;</mi> <mi>r</mi> </msub> <mo>=</mo> <mfrac> <mi>p</mi> <mn>1013.25</mn> </mfrac> <mn>0.008569</mn> <msup> <mi>&lambda;</mi> <mrow> <mo>-</mo> <mn>4</mn> </mrow> </msup> <mrow> <mo>(</mo> <mn>1</mn> <mo>+</mo> <mn>0.0133</mn> <msup> <mi>&lambda;</mi> <mrow> <mo>-</mo> <mn>2</mn> </mrow> </msup> <mo>+</mo> <mn>0.00013</mn> <msup> <mi>&lambda;</mi> <mrow> <mo>-</mo> <mn>4</mn> </mrow> </msup> <mo>)</mo> </mrow> </mrow>

其中，τ_r为瑞利散射光学厚度，λ为波长；p为MODIS的地表气压产品；

(6-6)τ＝τ_w+τ_o+τ_a+τ_r

τ是大气光学厚度；

第七步、计算太阳照射方向和卫星观测方向的大气透过率；

T_z＝exp(-τ/cosθ_z)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_z}

T_v＝exp(-τ/cosθ_v)＝exp{(-τ_r-τ_a-τ_o-τ_w)/cosθ_v}

第八步、计算地表反射率ρ

<mrow> <mi>&rho;</mi> <mo>=</mo> <mfrac> <mrow> <mi>&pi;</mi> <mrow> <mo>(</mo> <msub> <mi>L</mi> <mi>&lambda;</mi> </msub> <mo>-</mo> <msub> <mi>L</mi> <mi>p</mi> </msub> <mo>)</mo> </mrow> <msup> <mi>d</mi> <mn>2</mn> </msup> </mrow> <mrow> <msub> <mi>E</mi> <mn>0</mn> </msub> <msub> <mi>cos&theta;</mi> <mi>z</mi> </msub> <msub> <mi>T</mi> <mi>z</mi> </msub> <msub> <mi>T</mi> <mi>v</mi> </msub> </mrow> </mfrac> </mrow>

d＝1+0.0167sin[2π(J-93.5)/365]

E₀是大气层外相应波长的太阳光谱辐照度，d是日地距离，J是儒略日。