CN105261026A

CN105261026A - 一种星载多光谱相机的大气校正处理方法

Info

Publication number: CN105261026A
Application number: CN201510703415.XA
Authority: CN
Inventors: 杨磊; 傅俏燕; 张学文; 隋正伟; 高青山
Original assignee: China Center for Resource Satellite Data and Applications CRESDA
Current assignee: China Center for Resource Satellite Data and Applications CRESDA
Priority date: 2015-10-26
Filing date: 2015-10-26
Publication date: 2016-01-20
Anticipated expiration: 2035-10-26
Also published as: CN105261026B

Abstract

本发明公开了一种星载多光谱相机的大气校正处理方法，步骤如下：计算表观辐亮度；建立大气校正计算模型；测量并计算气溶胶光学厚度；大气校正快速计算。本发明对所在区域进行气溶胶光学厚度进行精确测量，获取得到高精度的气溶胶光学厚度参数，实现星载相机的大气校正，消除大气对遥感卫星影像的影响，增强卫星影像中地表反射率信息。

Description

一种星载多光谱相机的大气校正处理方法

技术领域

本发明涉及一种星载相机多光谱大气校正方法，可以实现星载相机的地表反射率生产，适用于在星载相机的实时大气校正处理。

背景技术

资源三号卫星(以下简称ZY-3卫星)是我国第一颗高分辨测绘卫星，能获取多光谱数据(MUX，5.8米)，载荷具体参数如表1所示。ZY-3卫星的高分辨多光谱数据在资源调查和监测等领域广泛应用，为防灾减灾、农林水利、环境保护、城市规划与建设、交通、国家重大工程等领域提供动态信息服务。开展ZY-3卫星数据大气校正工作并评价其反演地表参数的能力，为ZY-3卫星的业务应用提供高精度的大气校正产品，能够进一步提高ZY-3卫星的应用水平。

表1ZY-3多光谱相机技术参数

大气校正的方法主要有基于影像本身的暗目标法、基于地面的线性经验模型法和辐射传输模型法等3种方法。ZY-3卫星数据业务化应用要求从卫星数据本身来反演气溶胶光学厚度，进而完成大气校正处理，这就需要针对大气辐射特性设置用于反演大气参数的波段，而ZY-3卫星缺乏用于反演大气参数的波段，给ZY-3卫星数据的大气校正带来困难。

发明内容

本发明的技术解决问题是：针对现有技术的不足，提供了一种星载相机多光谱大气校正方法，本发明克服ZY-3卫星多光谱数据由于缺少短波红外数据而不能采用暗目标法计算气溶胶光学厚度的缺陷，通过建立卫星多光谱数据大气校正计算模型以及利用气溶胶光学厚度等参数计算地表反射率，实现在轨卫星数据高精度大气校正处理，提高了计算效率。

本发明的技术解决方案是：

一种星载多光谱相机的大气校正处理方法包括步骤如下：

(1)计算表观辐亮度：

L＝Gain×DN+Bias(1)

其中，L为表观辐亮度；Gain和Bias为星载相机增益和偏差值；DN为图像灰度值；

(2)建立大气校正计算模型：

L = L_{0} + \frac{ρ}{1 - s ρ} \frac{{TF}_{d}}{π} - - - (2)

其中，L为表观辐亮度；L₀为地表反射等于0时大气引起的程辐射；T为地表到传感器的大气透过率；F_d为到达地表的下行辐射通量密度；s为大气底层球面反照率；ρ为地表反射率；

(3)计算气溶胶光学厚度τ_aeraol(λ)：

τ_aeraol(λ)＝τ_total(λ)-τ_r(λ)-τ_oz(λ)

其中，τ_oz(λ)为臭氧的光学厚度；τ_r(λ)为大气分子瑞利散射光学厚度，

τ_{r} (λ) = 0.008569 λ^{- 4} (1 + 0.0113 λ^{- 2} + 0.0113 λ^{- 4}) * \frac{p}{1013.25} e^{- 0.125 h},

h为观测点的海拔，λ为波长；P为观测点在观测期间的大气气压；

τ_{t o t a l} (λ) = \frac{1}{m} ({lnE}_{0} (λ) + {lnd}_{s} - \ln E (λ)) = \frac{1}{m} l n \frac{E_{0} (λ) d_{s}}{E (λ)},

E₀(λ)为日地平均距离处地球大气上界的太阳单色直射辐照度，d_s为日地距离修正因子，E(λ)为通过地球到达地面的太阳单色直射辐照度；(4)大气校正计算

(4a)根据真实地表反射率的值域范围，设定3个地表反射率特征值，分别是ρ_a，ρ_b，ρ_c，将太阳-相机几何角度和步骤(3)得到的气溶胶光学厚度以及设定的地表反射率特征值代入MODTRAN辐射传输模型，计算得到表观辐亮度L_a，L_b，L_c；

(4b)将步骤(4a)设定的地表反射率特征值ρ_a，ρ_b，ρ_c以及计算得到的表观辐亮度L_a，L_b，L_c带入式(2)，求出大气系数s、L₀、T的值；

(4c)将步骤(1)中计算得到的实际观测点的表观辐亮度L和步骤(4b)求得大气系数s、L₀、T的代入(1b)中求得真实的地表反射率，完成对图像的逐像元进行大气校正计算。

步骤(3)中计算气溶胶光学厚度τ_aeraol(λ)的具体方式如下：

(3a)计算日地距离修正因子d_s：

d_{s} = {(\frac{r_{0}}{r})}^{2} = Σ_{n = 0}^{2} (a_{n} \cos n t + b_{n} \sin n t)

其中，r_o为日地平均距离；r为观测日期日地距离；n为正整数；a_n和b_n为日地距离修正因子；

(3b)计算太阳赤纬δ：

δ = Σ_{n = 0}^{3} (c_{n} \cos n t + d_{n} \sin n t)

其中，c_n和d_n为太阳赤纬计算系数，n为正整数；

(3c)计算观测时刻太阳时角ω：

ω = (μ t + \frac{φ}{15} + e - 12) \times 15

其中，μt为观测时刻的格林威治时间；φ为观测点经度，其每度对应的时间为4分钟，东经为正，西经为负；e为真太阳时与地方平均太阳时的修正系数，

e = \frac{24 (0.000075 + 0.001868 \cos t - 0.032077 \sin t - 0.014615 \cos 2 t - 0.04089 \sin 2 t)}{2 π},

d为积日，d＝0为1月1日，d＝364代表12月31日；

(3d)计算太阳天顶角θ_s：

其中，为观测地点的地理纬度；

(3e)计算观测点大气质量m：

m = \frac{1}{{cosθ}_{S} + 0.15 {(93.885 - θ_{S})}^{- 1.253}} * \frac{p}{1013.25};

(3f)计算吸收气体的透过率t_g(λ)：

t_g(λ)＝t_oz(λ)＝exp(-m*τ_oz(λ))

其中，t_oz(λ)为臭氧的透过率，τ_oz(λ)为臭氧的光学厚度，，α_oz(λ)为臭氧吸收系数，U为臭氧含量(DU)；

(3g)计算通过地球到达地面的太阳单色直射辐照度E(λ)：

E(λ)＝E₀(λ)d_sexp(-mτ_atm(λ))t_g(λ)

其中，E₀(λ)为日地平均距离处地球大气上界的太阳单色直射辐照度；d_s为日地距离修正因子；m为大气质量；τ_atm(λ)为大气光学厚度，主要是大气分子瑞利散射光学厚度和气溶胶光学厚度；

(3h)根据步骤(2g)和(2f)的公式可以得到下式：

E(λ)＝E₀(λ)d_sexp(-m*(τ_atm(λ)+τ_oz(λ)))

对上式两边取自然对数可得：

lnE(λ)＝lnE₀(λ)+lnd_s-m*(τ_atm(λ)+τ_oz(λ))

τ_total(λ)＝τ_atm(λ)+τ_oz(λ)＝τ_r(λ)+τ_aeraol(λ)+τ_oz(λ)

由上面两式得：

τ_{t o t a l} (λ) = \frac{1}{m} ({lnE}_{0} (λ) + {lnd}_{s} - \ln E (λ)) = \frac{1}{m} l n \frac{E_{0} (λ) d_{s}}{E (λ)}

太阳分光光度计输出的电压值(DN)与太阳辐照度之间是线性的：

\frac{E_{0} (λ)}{E (λ)} = \frac{{DN}_{0} (λ)}{D N (λ)}

τ_{t o t a l} (λ) = \frac{1}{m} l n \frac{{DN}_{0} (λ) d_{s}}{D N (λ)}

气溶胶的光学厚度：

τ_aeraol(λ)＝τ_total(λ)-τ_r(λ)-τ_oz(λ)。

本发明与现有技术相比的有益效果是：

(1)本发明能够大规模业务化生产遥感卫星数据的地表反射率产品，确定了采用太阳分光度计测量气溶胶光学厚度值和精确的大气计算模型，为精确计算地表反射率产品提供了条件，提高了计算效率，增强了工程适用性，通用性大大增强，在国内外属于首例。

(2)本发明中的气溶胶光学厚度采用先进的地基测量和计算，避免采用误差较大的气溶胶光学厚度值影响大气校正的精度，为最后的地表反射率的校正提供了准确的支持。

附图说明

图1为本发明流程图。

具体实施方式

下面结合附图对本发明的具体实施方式进行进一步的详细描述。

如图1所述，本发明一种星载多光谱相机的大气校正处理方法包括步骤如下：

(1)计算表观辐亮度：

L＝Gain×DN+Bias(1)

(2)建立大气校正计算模型：

L = L_{0} + \frac{ρ}{1 - s ρ} \frac{{TF}_{d}}{π} - - - (2)

(3)计算气溶胶光学厚度τ_aeraol(λ)：

τ_aeraol(λ)＝τ_total(λ)-τ_r(λ)-τ_oz(λ)

τ_{r} (λ) = 0.008569 λ^{- 4} (1 + 0.0113 λ^{- 2} + 0.0113 λ^{- 4}) * \frac{p}{1013.25} e^{- 0.125 h},

h为观测点的海拔(km)，λ为波长(μm)；P为观测点在观测期间的大气气压(Hpa)；

τ_{t o t a l} (λ) = \frac{1}{m} ({lnE}_{0} (λ) + {lnd}_{s} - \ln E (λ)) = \frac{1}{m} l n \frac{E_{0} (λ) d_{s}}{E (λ)},

E₀(λ)为日地平均距离处地球大气上界的太阳单色直射辐照度，d_s为日地距离修正因子，E(λ)为通过地球到达地面的太阳单色直射辐照度；计算气溶胶光学厚度τ_aeraol(λ)的具体方式如下：

(3a)计算日地距离修正因子d_s：

d_{s} = {(\frac{r_{0}}{r})}^{2} = Σ_{n = 0}^{2} (a_{n} \cos n t + b_{n} \sin n t)

其中，r_o为日地平均距离；r为观测日期日地距离；n为正整数；a_n和b_n为日地距离修正因子，具体取值如表1所示：

表1日地距离修正因子取值表

(3b)计算太阳赤纬δ：

δ = Σ_{n = 0}^{3} (c_{n} \cos n t + d_{n} \sin n t)

其中，c_n和d_n为太阳赤纬计算系数，n为正整数；

表2计算太阳赤纬的系数

(3c)计算观测时刻太阳时角ω：

ω = (μ t + \frac{φ}{15} + e - 12) \times 15

e = \frac{24 (0.000075 + 0.001868 \cos t - 0.032077 \sin t - 0.014615 \cos 2 t - 0.04089 \sin 2 t)}{2 π},

d为积日，d＝0为1月1日，d＝364代表12月31日；

(3d)计算太阳天顶角θ_s：

其中，为观测地点的地理纬度；

(3e)计算观测点大气质量m：

m = \frac{1}{{cosθ}_{S} + 0.15 {(93.885 - θ_{S})}^{- 1.253}} * \frac{p}{1013.25};

(3f)计算吸收气体的透过率t_g(λ)：

t_g(λ)＝t_oz(λ)＝exp(-m*τ_oz(λ))

(3g)计算通过地球到达地面的太阳单色直射辐照度E(λ)(W/m²)：

E(λ)＝E₀(λ)d_sexp(-mτ_atm(λ))t_g(λ)(2a)

(3h)根据步骤(2g)和(2f)的公式可以得到下式：

E(λ)＝E₀(λ)d_sexp(-m*(τ_atm(λ)+τ_oz(λ)))

对上式两边取自然对数可得：

lnE(λ)＝lnE₀(λ)+lnd_s-m*(τ_atm(λ)+τ_oz(λ))

τ_total(λ)＝τ_atm(λ)+τ_oz(λ)＝τ_r(λ)+τ_aeraol(λ)+τ_oz(λ)

由上面两式得：

τ_{t o t a l} (λ) = \frac{1}{m} ({lnE}_{0} (λ) + {lnd}_{s} - \ln E (λ)) = \frac{1}{m} l n \frac{E_{0} (λ) d_{s}}{E (λ)}

\frac{E_{0} (λ)}{E (λ)} = \frac{{DN}_{0} (λ)}{D N (λ)}

τ_{t o t a l} (λ) = \frac{1}{m} l n \frac{{DN}_{0} (λ) d_{s}}{D N (λ)}

气溶胶的光学厚度：

τ_aeraol(λ)＝τ_total(λ)-τ_r(λ)-τ_oz(λ)。(4)大气校正计算

(4b)将步骤(4a)设定的地表反射率特征值ρ_a，ρ_b，ρ_c以及计算得到的表观辐亮度L_a，L_b，L_c带入式(1b)，求出大气系数s、L₀、T的值；

卫星采集时太阳-相机的几何角度从卫星数据的辅助数据获取。

太阳-相机几何角度可以从卫星标准数据产品的辅助数据获取，分别是：

(1)太阳方位角(SolarAzimuth)

(2)太阳天顶角(SolarZenith)

(3)卫星方位角(SatelliteAzimuth)

(4)卫星天顶角(SatelliteZenith)

(5)卫星相机侧摆角(YawSatelliteAngle)

选取利用敦煌地区的ZY-3卫星多光谱数据，采用本方法进行大气校正。计算结果表明，未经大气校正的影像的光谱曲线变化平缓，经过大气校正后，多光谱相机影像的蓝波段(band1)的反射率大大降低，红波段(band3)和红外波段(band4)的反射率提高明显，且与地面实测的曲线一致。结果表明多光谱相机的大气校正效果较好。

表1ZY-3大气校正效果

本发明说明书中未作详细描述的内容属于本领域专业技术人员的公知技术。

Claims

1.一种星载多光谱相机的大气校正处理方法，其特征在于步骤如下：

(1)计算表观辐亮度：

L＝Gain×DN+Bias(1)

(2)建立大气校正计算模型：

L = L_{0} + \frac{ρ}{1 - s ρ} \frac{{TF}_{d}}{π} - - - (2)

(3)计算气溶胶光学厚度τ_aeraol(λ)：

τ_aeraol(λ)＝τ_total(λ)-τ_r(λ)-τ_oz(λ)

τ_{r} (λ) = 0.008569 λ^{- 4} (1 + 0.0113 λ^{- 2} + 0.0113 λ^{- 4}) * \frac{p}{1013.25} e^{- 0.125 h},

τ_{t o t a l} (λ) = \frac{1}{m} (\ln E_{0} (λ) + \ln d_{s} - \ln E (λ)) = \frac{1}{m} l n \frac{E_{0} (λ) d_{s}}{E (λ)},

E₀(λ)为日地平均距离处地球大气上界的太阳单色直射辐照度，d_s为日地距离修正因子，E(λ)为通过地球到达地面的太阳单色直射辐照度；

(4)大气校正计算

(4c)将步骤(1)中计算得到的实际观测点的表观辐亮度L和步骤(4b)求得大气系数s、L₀、T的代入(2)中求得真实的地表反射率，完成对图像的逐像元进行大气校正计算。

2.根据权利要求1所述的一种星载多光谱相机的大气校正处理方法，其特征在于：步骤(3)中计算气溶胶光学厚度τ_aeraol(λ)的具体方式如下：

(3a)计算日地距离修正因子d_s：

d_{s} = {(\frac{r_{0}}{r})}^{2} = Σ_{n = 0}^{2} (a_{n} \cos n t + b_{n} \sin n t)

(3b)计算太阳赤纬δ：

δ = Σ_{n = 0}^{3} (c_{n} \cos n t + d_{n} \sin n t)

其中，c_n和d_n为太阳赤纬计算系数，n为正整数；

(3c)计算观测时刻太阳时角ω：

ω = (μ t + \frac{φ}{15} + e - 12) \times 15

e = \frac{24 (0.000075 + 0.001868 \cos t - 0.032077 \sin t - 0.014615 \cos 2 t - 0.04089 \sin 2 t)}{2 π},

d为积日，d＝0为1月1日，d＝364代表12月31日；

(3d)计算太阳天顶角θ_s：

其中，为观测地点的地理纬度；

(3e)计算观测点大气质量m：

m = \frac{1}{{cosθ}_{s} + 0.15 (93.885 - θ_{S}) - 1.253} * \frac{p}{1013.25};

(3f)计算吸收气体的透过率t_g(λ)：

t_g(λ)＝t_oz(λ)＝exp(-m*τ_oz(λ))

其中，t_oz(λ)为臭氧的透过率，τ_oz(λ)为臭氧的光学厚度， α_oz(λ)为臭氧吸收系数，U为臭氧含量(DU)；

(3g)计算通过地球到达地面的太阳单色直射辐照度E(λ)：

E(λ)＝E₀(λ)d_sexp(-mτ_atm(λ))t_g(λ)(2a)

(3h)根据步骤(2g)和(2f)的公式可以得到下式：

E(λ)＝E₀(λ)d_sexp(-m*(τ_atm(λ)+τ_oz(λ)))

对上式两边取自然对数可得：

lnE(λ)＝lnE₀(λ)+lnd_s-m*(τ_atm(λ)+τ_oz(λ))

τ_total(λ)＝τ_atm(λ)+τ_oz(λ)＝τ_r(λ)+τ_aeraol(λ)+τ_oz(λ)

由上面两式得：

τ_{t o t a l} (λ) = \frac{1}{m} (\ln E_{0} (λ) + \ln d_{s} - \ln E (λ)) = \frac{1}{m} l n \frac{E_{0} (λ) d_{s}}{E (λ)}

\frac{E_{0} (λ)}{E (λ)} = \frac{{DN}_{0} (λ)}{D N (λ)}

τ_{t o t a l} (λ) = \frac{1}{m} l n \frac{{DN}_{0} (λ) d_{s}}{D N (λ)}

气溶胶的光学厚度：

τ_aeraol(λ)＝τ_total(λ)-τ_r(λ)-τ_oz(λ)。

3.根据权利要求1所述的一种星载多光谱相机的大气校正处理方法，其特征在于：步骤(4a)中的卫星采集时太阳-相机的几何角度从卫星数据的辅助数据获取。