CN104133203A - 一种立体测绘卫星三线阵ccd相机的大气订正方法 - Google Patents

Abstract

一种立体测绘卫星三线阵CCD相机的大气订正方法，(1)计算每个像元的入瞳处辐射亮度；并确立每个像元的地表反射率计算模型；(2)建立不同气溶胶类型和气溶胶光学厚度下的立体测绘卫星三线阵CCD相机每个波段的反演参数查找表；(3)依据Rahman BRDF模型，建立立体测绘卫星前视、后视CCD相机与正视CCD相机地表反射率关系；并利用地表反射率计算模型计算三个相机不同气溶胶类型和光学厚度下的地表反射率；(4)建立最小化判定函数，确定函数的最优解，根据最优解确定三线阵CCD相机的气溶胶类型和光学厚度；(5)对步骤(4)的结果进行修正，获得三线阵CCD相机修正后的气溶胶光学厚度；(6)利用地表反射率计算模型计算得到最终的地表反射率。

Description

一种立体测绘卫星三线阵CCD相机的大气订正方法

技术领域

本发明涉及立体测绘卫星三线阵CCD相机的大气订正方法，特别涉及ZY-3立体测绘卫星三线阵CCD相机的大气订正方法。

背景技术

立体测绘卫星三线阵CCD采用3台线阵CCD(前视相机、正视相机、后视相机)在不同时刻对同一地物进行多角度成像，不仅可以获得目标立体信息，同时还可以得到多角度的定量化信息，但是要获得高精度的定量化信息，就必须详细考虑大气对于成像波段传输过程的定量影响，并对其影响做出准确的大气订正，实现地物目标辐射特性真实再现，以达到提高遥感信息质量及定量化应用水平目的。

目前国内外大气订正方面的成果较多，但由于大气状况瞬息万变，当时、当地的大气对遥感数据影响程度和重要性不同，再加上各卫星的相机波段设置不同，特性不同，研究目的和要求也不同，即使相同地区，大气状况也是复杂多变的，每景图像成像时的大气状况也存在明显差异，因此很难将某种大气订正方法直接应用于立体测绘卫星三线阵CCD相机。

发明内容

本发明的技术解决问题是：克服现有技术的不足，提供一种易于业务化运行的大气订正方法。

本发明的技术解决方案是：一种立体测绘卫星三线阵CCD相机的大气订正方法，所述的三线阵CCD相机包括前视、正视、后视三台CCD相机；步骤如下：

(1)根据立体测绘卫星三线阵CCD相机的遥感图像数据，计算每个像元的入瞳处辐射亮度；并依据太阳-大气-地表-大气-相机相互作用的大气辐射传输方程确立每个像元的地表反射率计算模型；

(2)依据步骤(1)中建立的地表反射率计算模型，设置气溶胶类型和气溶胶的光学厚度，建立不同气溶胶类型和气溶胶光学厚度下的立体测绘卫星三线阵CCD相机每个波段的反演参数查找表；

(3)依据Rahman BRDF模型，建立立体测绘卫星前视CCD相机与正视CCD相机地表反射率关系，以及后视CCD相机与正视CCD相机地表反射率关系；并利用上述反演参数查找表、每个像元入瞳处辐射亮度以及地表反射率计算模型计算三个相机不同气溶胶类型和光学厚度下的地表反射率；

(4)根据步骤(3)的结果，建立最小化判定函数，确定函数的最优解，根据最优解确定三线阵CCD相机的气溶胶类型和光学厚度；

(5)根据气溶胶光学厚度与波长的关系对步骤(4)的结果进行修正，获得三线阵CCD相机修正后的气溶胶光学厚度；

(6)根据步骤(4)中确定的气溶胶类型以及步骤(5)中获得的气溶胶光学厚度，代入上述反演参数查找表，利用地表反射率计算模型计算得到最终的地表反射率。

所述步骤(1)中每个像元的地表反射率计算模型如下：

${\mathrm{\ρ}}_{x,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)=\frac{L_{x,\mathrm{\λ}}^{\mathrm{TOA}}(i,j)-L_{x,\mathrm{\λ}}^{\mathrm{Path}}(i,j)}{T_{x,\mathrm{\λ}}(i,j)+S_{x,\mathrm{\λ}}(i,j)[L_{x,\mathrm{\λ}}^{\mathrm{TOA}}(i,j)-L_{x,\mathrm{\λ}}^{\mathrm{Path}}(i,j)]}$

其中，为x相机第λ波段像元(i,j)处的地表反射率；x＝f，n，b代表立体测绘卫星的前视、正视、后视相机；i，j为每台相机的像元位置；

为x相机第λ波段像元(i,j)处的入瞳处辐射亮度；

为x相机第λ波段像元(i,j)处的路径散射；

T(θ_x,s)为太阳到地面方向的大气漫射散射因子；T(θ_x,v)为地面到相机方向的大气漫射散射因子；θ_x,s为x相机观测时的太阳天顶角；E_S,x,λ为x相机第λ波段的大气外太阳光谱辐照度E_S；

S_x,λ(i,j)为x相机第λ波段像元(i,j)处的大气球面反照率。

所述步骤(4)中的最小化判定函数如下：

$\min \left|Y_{\mathrm{\λ}}(i,j)\right|=\underset{x=f,b}{\mathrm{\Σ}}{[{\mathrm{\ρ}}_{x,\mathrm{\λ}}^{R-\mathrm{Obj}}(i,j)-{\mathrm{\ρ}}_{x,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)]}^2$

其中，为利用步骤(3)中的地表反射率关系以及步骤(3)中计算的正视CCD相机不同气溶胶类型和光学厚度下的地表反射率，得到的x相机的地表反射率，x＝f，b分别代表立体测绘卫星的前视、后视相机；

为步骤(3)中计算的x相机不同气溶胶类型和光学厚度下的地表反射率。

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

(1)本发明利用立体观测卫星三线阵CCD相机的多角度观测特点，结合大气辐射传输模型和Rahman BRDF模型，实现了立体观测卫星三线阵CCD相机的大气订正。

(2)本发明提出的立体观测卫星三线阵CCD相机大气订正方法，实现了逐像元的大气订正，克服以往大气订正算法中多区域采用单一气溶胶光学厚度的缺点，有效地提高了大气订正的精度。

(3)本发明充分利用了立体观测卫星三线阵CCD相机的多角度观测特点，直接从多角度遥感图像本身出发进行大气订正，不要任何实际测量参数或其他辅助参数，实施简易、快速、经济使得工程化的业务化大气订正得以实现。

附图说明

图1为本发明方法流程图。

具体实施方式

下面结合附图对本发明做详细说明，如图1所示，本发明步骤具体如下：

(1)相机入瞳处辐射亮度计算：读取立体测绘卫星三线阵CCD相机的遥感图像数据，依据数据附带文件中的绝对辐射定标系数，按照公式①计算每个像元的入瞳处辐射亮度：

$L_{x,\mathrm{\λ}}^{\mathrm{TOA}}(i,j)={\mathrm{Gain}}_{x,\mathrm{\λ}}\×{\mathrm{DN}}_{x,\mathrm{\λ}}(i,j)+{\mathrm{Bias}}_{x,\mathrm{\λ}}$    ①

式中：x＝f，n，b代表立体测绘卫星的前视、正视、后视相机；

λ代表每台相机的不同波段；

i，j为每台相机的像元位置(i∈[1,K]，j∈[1,L])；

为x相机第λ波段像元(i,j)处的入瞳处辐射亮度；

DN_x,λ(i,j)为x相机第λ波段像元(i,j)处的数字DN值；

Gain_x,λ和Bias_x,λ为x相机第λ波段的绝对辐射定标系数。

(2)地表反射率计算模型确立：依据太阳-大气-地表-大气-相机相互作用的大气辐射传输方程，按照公式②确立每个像元的地表反射率计算模型：

$L_{x,\mathrm{\λ}}^{\mathrm{TOA}}(i,j)=L_{x,\mathrm{\λ}}^{\mathrm{path}}(i,j)+\frac{E_{S,x,\mathrm{\λ}}\cos \left({\mathrm{\θ}}_{x,s}\right)T\left({\mathrm{\θ}}_{x,s}\right)T\left({\mathrm{\θ}}_{x,v}\right)}{\mathrm{\π}}\×\frac{{\mathrm{\ρ}}_{x,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)}{1-{\mathrm{\ρ}}_{x,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)S_{x,\mathrm{\λ}}(i,j)}$    ② $=L_{x,\mathrm{\λ}}^{\mathrm{Path}}(i,j)+\frac{T_{x,\mathrm{\λ}}(i,j){\mathrm{\ρ}}_{x,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)}{1-{\mathrm{\ρ}}_{x,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)S_{x,\mathrm{\λ}}(i,j)}$

式中：为x相机第λ波段像元(i,j)处的路径散射；

E_S,x,λ为x相机第λ波段的大气外太阳光谱辐照度E_S；

θ_x,s为x相机观测时的太阳天顶角

θ_x,v为x相机观测时的相机观测天顶角；

φ_x,s为x相机观测时的太阳方位角；

φ_x,v为x相机观测时的相机观测方位角；

T(θ_x,s)为向下方向(太阳到地面)的大气漫射散射因子；

T(θ_x,v)为向上方向(地面到相机)的大气漫射散射因子；

为x相机第λ波段像元(i,j)处的地表反射率；

S_x,λ(i,j)为x相机第λ波段像元(i,j)处的大气球面反照率；

$T_{x,\mathrm{\λ}}(i,j)=\frac{E_{S,x,\mathrm{\λ}}\cos \left({\mathrm{\θ}}_{x,s}\right)T\left({\mathrm{\θ}}_{x,s}\right)T\left({\mathrm{\θ}}_{x,v}\right)}{\mathrm{\π}}.$

对②进行变化得到地表反射率计算模型：

${\mathrm{\ρ}}_{x,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)=\frac{L_{x,\mathrm{\λ}}^{\mathrm{TOA}}(i,j)-L_{x,\mathrm{\λ}}^{\mathrm{Path}}(i,j)}{T_{x,\mathrm{\λ}}(i,j)+S_{x,\mathrm{\λ}}(i,j)[L_{x,\mathrm{\λ}}^{\mathrm{TOA}}(i,j)-L_{x,\mathrm{\λ}}^{\mathrm{Path}}(i,j)]}$    ③

由上知已确定，故地表反射率由T_x,λ(i,j)和S_x,λ(i,j)参数确定；这3个反演参数是相机波段、太阳天顶角、卫星观测天顶角、太阳方位角、卫星观测方位角、气溶胶类型及气溶胶光学厚度τ的函数即有 $\{L_{x,\mathrm{\λ}}^{\mathrm{Path}}(i,j),T_{x,\mathrm{\λ}}(i,j),S_{x,\mathrm{\λ}}(i,j)\}=f(\mathrm{\λ},{\mathrm{\θ}}_{x,s},{\mathrm{\θ}}_{x,v},{\mathrm{\φ}}_{x,s},{\mathrm{\φ}}_{x,v},\mathrm{\τ}),$ 卫星成像时除气溶胶类型及光学厚度的参数外其他参数都唯一确定，那么在其他参数状态确定的前提下这3个反演参数由气溶胶类型及光学厚度确定。

(3)参数的查找表建立：依据步骤(2)中确立的参数模型，设置气溶胶类型和气溶胶光学厚度，通过大气辐射传输模型为立体测绘卫星每台线阵CCD相机的每个波段建立以(λ,θ_x,s,θ_x,v,φ_x,s,φ_x,v,τ)为索引的T_x,λ(i,j)和S_x,λ(i,j)在不同参数状态下的反演参数查找表值；

(4)气溶胶类型及气溶胶光学厚度确定和修正：依据Rahman BRDF模型，建立立体测绘卫星前后、正视、后视相机的地表反射率关系以确定判定函数，通过判定函数确定气溶胶类型及气溶胶光学厚度，最后利用公式修正气溶胶光学厚度值；

(4.1)基于Rahman BRDF模型计算地表反射率

依据Rahman BRDF模型，立体测绘卫星前视、正视、后视相机的地表反射率关系为：

${\mathrm{\ρ}}_{f,\mathrm{\λ}}^{R-\mathrm{Obj}}(i,j)={\mathrm{\ρ}}_{n,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)\frac{{\left(\cos {\mathrm{\θ}}_{f,s}\right)}^{k-1}{\left(\cos {\mathrm{\θ}}_{f,v}\right)}^{k-1}}{{(\cos {\mathrm{\θ}}_{f,s}+\cos {\mathrm{\θ}}_{f,v})}^{1-k}}{F}_f\left(\mathrm{\Θ}\right)[1+H_{f,\mathrm{\λ}}\left(G\right)]$

${\mathrm{\ρ}}_{b,\mathrm{\λ}}^{R-\mathrm{Obj}}(i,j)={\mathrm{\ρ}}_{n,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)\frac{{\left(\cos {\mathrm{\θ}}_{b,s}\right)}^{k-1}{\left(\cos {\mathrm{\θ}}_{b,v}\right)}^{k-1}}{{(\cos {\mathrm{\θ}}_{b,s}+\cos {\mathrm{\θ}}_{b,v})}^{1-k}}{F}_b\left(\mathrm{\Θ}\right)[1+H_{b,\mathrm{\λ}}\left(G\right)]$    ④

式中：k为常数k∈[0,1]，表征目标表面各向异性的特性；

F(Θ)为气溶胶散射特性的相函数；

$F\left({\mathrm{\Θ}}_f\right)=\frac{1-{g_f}^2}{{[1+{g_f}^2-{2g}_f\cos \left({\mathrm{\Θ}}_f\right)]}^{1.5}}$

$F\left({\mathrm{\Θ}}_b\right)=\frac{1-{g_b}^2}{{[1+{g_b}^2-{2g}_b\cos \left({\mathrm{\Θ}}_b\right)]}^{1.5}}$

Θ_x为相函数的相位角

$\cos \left({\mathrm{\Θ}}_f\right)=\cos \left({\mathrm{\θ}}_{f,s}\right)\cos \left({\mathrm{\θ}}_{f,v}\right)-\sin \left({\mathrm{\θ}}_{f,s}\right)\sin \left({\mathrm{\θ}}_{f,v}\right)\cos ({\mathrm{\φ}}_{f,s}-{\mathrm{\φ}}_{f,v})$

$\cos \left({\mathrm{\Θ}}_b\right)=-\cos \left({\mathrm{\θ}}_{b,s}\right)\cos \left({\mathrm{\θ}}_{b,v}\right)-\sin \left({\mathrm{\θ}}_{b,s}\right)\sin \left({\mathrm{\θ}}_{b,v}\right)\cos ({\mathrm{\φ}}_{b,s}-{\mathrm{\φ}}_{b,v})$

g_x为常数依据气溶胶类型设置，表征不对称因子；

g_f∈[0,1]，g_b∈[-1,0]

H_x,λ(G)为BRDF的热点效应；

$H_{x,\mathrm{\λ}}\left(G\right)=\frac{1-{\mathrm{\ρ}}_{n,\mathrm{\λ}}^{\mathrm{Obj}}}{1+G_x}$

G_x为几何因子；

$G_x=\sqrt{{\left[\tan \left({\mathrm{\θ}}_{x,s}\right)\right]}^2+{\left[\tan \right]\left({\mathrm{\θ}}_{x,v}\right)}^2-2\tan \left({\mathrm{\θ}}_{x,s}\right)\tan \left({\mathrm{\θ}}_{f,v}\right)\cos ({\mathrm{\φ}}_{x,s}-{\mathrm{\φ}}_{x,v})}$

(4.2)判定函数确立

利用Rahman BRDF模型确立的立体测绘卫星前后、正视、后视相机地表反射率关系得到的与大气辐射传输模型确立的地表反射率的等价关系，建立最小化判定函数为：

$\min \left|Y_{\mathrm{\λ}}(i,j)\right|=\underset{x=f,b}{\mathrm{\Σ}}\left|{[{\mathrm{\ρ}}_{x,\mathrm{\λ}}^{R-\mathrm{Obj}}(i,j)-{\mathrm{\ρ}}_{x,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)]}^2\right|$    ⑤

其中，为利用步骤(3)中的地表反射率关系以及步骤(3)中计算的正视CCD相机不同气溶胶类型和光学厚度下的地表反射率，得到的x相机的地表反射率，x＝f，b分别代表立体测绘卫星的前视、后视相机；为步骤(3)中计算的x相机不同气溶胶类型和光学厚度下的地表反射率。

(4.3)气溶胶类型及气溶胶光学厚度确立

利用立体测绘卫星成像时获得的θ_x,s,θ_x,v,φ_x,s,φ_x,v，通过步骤(3)反演参数查找表值获得以(λ,θ_x,s,θ_x,v,φ_x,s,φ_x,v)索引的不同气溶胶类型及光学厚度下的T_x,λ(i,j)和S_x,λ(i,j)，并根据大气辐射传输模型确立公式③计算出立体测绘卫星三线阵CCD相机对应的地表反射率再根据步骤(4)中确定的地表反射率关系④计算出不同气溶胶类型及光学厚度下正视相机地表反射率对应的前视、后视相机的地表反射率和将上述两种方式计算多种气溶胶类型和光学厚度条件下的立体测绘卫星三线阵CCD相机的地表反射率和代入方程⑤式获得各种气溶胶类型和光学厚度条件下的判定函数值，将满足判定函数最小值对应条件下立体测绘卫星三线阵CCD相机各地表反射率对应的气溶胶类型和光学厚度确定为函数最优解；

(4.4)气溶胶类型及气溶胶光学厚度确修正

根据步骤(4.3)确定的立体测绘卫星三线阵CCD相机对应的气溶胶光学厚度，通过气溶胶光学厚度与波长满足的公式τ(λ)＝αλ^-β进行最小二乘法拟合以确定公式系数α、β，最后通过已知参数α、β的公式获得修正后立体测绘卫星三线阵CCD相机不同波段的气溶胶光学厚度。

(5)大气订正处理：根据步骤(4)获得的气溶胶类型及修正的气溶胶光学厚度，结合立体观测卫星三线阵CCD相机的几何参数θ_x,s,θ_x,v,φ_x,s,φ_x,v，查找出以(λ,θ_x,s,θ_x,v,φ_x,s,φ_x,v)为索引的T_x,λ(i,j)和S_x,λ(i,j)反演参数值，即可通过方程③式获得经过大气订正后的立体测绘卫星三线阵CCD相机对应的地表反射率。

本发明未详细说明部分属于本领域技术人员公知常识。

Claims (3)

1.一种立体测绘卫星三线阵CCD相机的大气订正方法，所述的三线阵CCD相机包括前视、正视、后视三台CCD相机；其特征在于步骤如下：

(1)根据立体测绘卫星三线阵CCD相机的遥感图像数据，计算每个像元的入瞳处辐射亮度；并依据太阳-大气-地表-大气-相机相互作用的大气辐射传输方程确立每个像元的地表反射率计算模型；

(2)依据步骤(1)中建立的地表反射率计算模型，设置气溶胶类型和气溶胶的光学厚度，建立不同气溶胶类型和气溶胶光学厚度下的立体测绘卫星三线阵CCD相机每个波段的反演参数查找表；

(3)依据Rahman BRDF模型，建立立体测绘卫星前视CCD相机与正视CCD相机地表反射率关系，以及后视CCD相机与正视CCD相机地表反射率关系；并利用上述反演参数查找表、每个像元入瞳处辐射亮度以及地表反射率计算模型计算三个相机不同气溶胶类型和光学厚度下的地表反射率；

(4)根据步骤(3)的结果，建立最小化判定函数，确定函数的最优解，根据最优解确定三线阵CCD相机的气溶胶类型和光学厚度；

(5)根据气溶胶光学厚度与波长的关系对步骤(4)的结果进行修正，获得三线阵CCD相机修正后的气溶胶光学厚度；

(6)根据步骤(4)中确定的气溶胶类型以及步骤(5)中获得的气溶胶光学厚度，代入上述反演参数查找表，利用地表反射率计算模型计算得到最终的地表反射率。

2.根据权利要求1所述的一种立体测绘卫星三线阵CCD相机的大气订正方法，其特征在于：所述步骤(1)中每个像元的地表反射率计算模型如下：

${\mathrm{\ρ}}_{x,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)=\frac{L_{x,\mathrm{\λ}}^{\mathrm{TOA}}(i,j)-L_{x,\mathrm{\λ}}^{\mathrm{Path}}(i,j)}{T_{x,\mathrm{\λ}}(i,j)+S_{x,\mathrm{\λ}}(i,j)[L_{x,\mathrm{\λ}}^{\mathrm{TOA}}(i,j)-L_{x,\mathrm{\λ}}^{\mathrm{Path}}(i,j)]}$

其中，为x相机第λ波段像元(i,j)处的地表反射率；x＝f，n，b代表立体测绘卫星的前视、正视、后视相机；i，j为每台相机的像元位置；

为x相机第λ波段像元(i,j)处的入瞳处辐射亮度；

为x相机第λ波段像元(i,j)处的路径散射；

T(θ_x,s)为太阳到地面方向的大气漫射散射因子；T(θ_x,v)为地面到相机方向的大气漫射散射因子；θ_x,s为x相机观测时的太阳天顶角；E_S,x,λ为x相机第λ波段的大气外太阳光谱辐照度E_S；

S_x,λ(i,j)为x相机第λ波段像元(i,j)处的大气球面反照率。

3.根据权利要求1所述的一种立体测绘卫星三线阵CCD相机的大气订正方法，其特征在于：所述步骤(4)中的最小化判定函数如下：

$\min \left|Y_{\mathrm{\λ}}(i,j)\right|=\underset{x=f,b}{\mathrm{\Σ}}{[{\mathrm{\ρ}}_{x,\mathrm{\λ}}^{R-\mathrm{Obj}}(i,j)-{\mathrm{\ρ}}_{x,\mathrm{\λ}}^{\mathrm{Obj}}(i,j)]}^2$

其中，为利用步骤(3)中的地表反射率关系以及步骤(3)中计算的正视CCD相机不同气溶胶类型和光学厚度下的地表反射率，得到的x相机的地表反射率，x＝f，b分别代表立体测绘卫星的前视、后视相机；

为步骤(3)中计算的x相机不同气溶胶类型和光学厚度下的地表反射率。