CN105403201B - 一种基于像元分解的遥感图像大气程辐射获取方法 - Google Patents

Abstract

本发明公开了一种基于像元分解的遥感图像大气程辐射获取方法，它是在像元分解技术的基础上，通过光谱端元概括的方法，在像元尺度上，逐点获得卫星遥感影像的大气程辐射。本发明的方法仅仅依靠遥感影像本身即可逐点获得遥感影像大气程辐射，不需要大气参数等其它辅助数据，且最后得到的结果可以很好的反映城市市区尺度的大气灰霾浓度的空间差异，从而可用于实现遥感监测城市气溶胶污染地区差异的目标。

Description

一种基于像元分解的遥感图像大气程辐射获取方法

技术领域

本发明涉及可以应用在市区或者乡镇范围尺度的城乡建设、农业、林业、气象、生态环境等行业部门的一种基于像元分解的遥感图像大气程辐射获取方法，属于遥感成像技术领域。

背景技术

电磁波在大气传输过程中，因为大气中的气体分子、气溶胶等物质存在容易发生散射、透过等衰减作用，降低传感器接收的电磁波辐射信号强度(如图1所示)。尤其是因为散射造成地物边界模糊、对比度降低等而影响图像质量，该进入传感器的部分称为大气程辐射，对于传统的遥感专题分类研究而言是一个非常大的“噪声”，需要大气校正来提高遥感图像质量。

反向思维，大气程辐射是一个仅反映大气状况、而与地面信息无关的辐射能量信号，对于遥感监测大气气溶胶浓度的空间分布格局有着重要的意义。

传统的去除大气程辐射的大气校正方法，主要为两类方法，其一为：暗目标法、波段归一化比值运算法、直方图匹配法等，这些方法虽然较为简单，但是整幅图像只有一个程辐射值，运用在城市区域显然难以模拟城市大气气溶胶污染的空间差异；而另一种方法为：裂窗法、6S模型、LOWTRAN模型、MODTRAN模型、ATCOR模型等，这些模型校正结果较为精确，但是需要相关大气参数。

张兆明等人提出的发明专利“一种逐像元计算卫星遥感影像大气程辐射的物理方法”，它是基于遥感成像物理原理，在MODIS光学厚度相关波段信息基础上进行逐点计算，是一个较为实用的大气程辐射精确值计算的方法。但是该方法基于MODIS数据，计算结果最后呈现出来的效果对于研究城市市区尺度的大气灰霾浓度的空间差异而言，适用性受到很大限制。难以应用到城市规划、产业布局规划的管理决策中。

李先华、徐丽华等人提出的“基于土地利用分类图和同期星下点反射率实测方法，逐点计算大气程辐射遥感值”的方法，也是基于遥感成像物理原理，利用TM和MODIS数据，逐点计算大气程辐射精确值的方法。但是该方法基于研究区土地利用分类图和同期、同位置的地物波谱实测，计算结果最后呈现出来的效果可以很好的反映城市市区尺度的大气灰霾浓度的空间差异，但是按照该方法每次计算都需要进行地物波谱实测和获取研究区土地利用现状图，使得后期推广受限。

本发明从遥感卫星影像的大气程辐射计算的理论基础出发，在前面的基础上，基于研究区的地物波谱库基础数据库的建立，提出了一种基于像元分解的遥感影像大气程辐射的逐点计算方法，该方法完全基于遥感影像本身，不需要辅助数据，较容易掌握和理解，对于监测城市灰霾天气有着重要意义。

发明内容

本发明的目的在于，提供一种基于像元分解的遥感图像大气程辐射获取方法。本发明的方法仅仅依靠遥感影像本身即可逐点获得遥感影像大气程辐射，不需要大气参数等其它辅助数据，且最后得到的结果可以很好的反映城市市区尺度的大气灰霾浓度的空间差异，从而可用于实现遥感监测城市气溶胶污染地区差异的目标。

本发明的技术方案：一种基于像元分解的遥感图像大气程辐射获取方法，其特点是,包括以下步骤：

第一步、对原始遥感影像(针对研究区的)进行几何校正。

第二步、将采集的地面波谱数据进行归纳为低照度不透水面、高照度不透水面和植被三种类型的光谱端元，并在遥感影像的波段区间获得各光谱端元在波段区间内的平均反射率；其中的不透水面包括干燥裸露的土壤；各光谱端元在各波段区间的平均反射率_ρi可按公式

(1)计算： ${\mathrm{\ρ}}_i=\frac{{\∫}_a^{b}r\left(x\right)d_x}{(b-a)}---\left(1\right)$

在公式(1)中，r(x)为实测的波谱数据向量，[a,b]为影像波段区间范围。

第三步、基于像元分解方法，获取原始遥感影像中植被、高照度不透水面、低照度不透水面三个光谱端元比例分布图；计算是采用基于约束性的混合像元分解方法(下面的公式2)，得到研究区的植被、高照度不透水面、低照度不透水面的比例图像。

${\mathrm{\ρ}}_i=\underset{i=1}{\overset{n}{\Σ}}\left(a_{ij}{x}_j\right)+e_i---\left(2\right)$

且0≤x_j≤1 (3)

$RMS=\sqrt{(\underset{i=1}{\overset{m}{\Σ}}{e}_i^2)/m}---\left(4\right)$

其中ρ_i为包含一个或多个组分的像元在第i个波段上的平均光谱反射率；α_ij为像元内第i个波段上第j个组分的反射率；x_j为像元内第j个组分所占的比例；e_i为第i个光谱波段的误差项；m(i＝1,2,,m)为传感器系统的波段数；RMS≤0.05为最小均方根误差。

第四步、结合第二部和第三部，得到原始遥感影像的地面反射率图像；

经过第三步，获取植被、低照度、高照度三个光谱端元的比例分配遥感图。

在第二步基础上，基于像元尺度，逐点计算第i个波段的地面反射率值R_i。

R_i＝ρ_i×x_j (5)

⑤在遥感成像的原理下，根据步骤④得到的地面反射率图像，对像元进行逐点分解，得到每个像元的大气程辐射遥感值。

上述的基于像元分解的遥感图像大气程辐射获取方法中，在进行步骤⑤的大气程辐射遥感值的计算前，根据步骤④获得的地面反射率，划定不同地物类型的边界阈值，获得原始遥感影像所示的研究区的地物类型边界线，从而将整幅图分为同一图斑和不同图斑边界线两个区域分别计算大气程辐射遥感值。

图斑边界线上的大气程辐射遥感值根据公式(6)计算。

DNa＝DN₁-r₁×(DN₁-DN₂)/(r₁-r₂) (6)

其中DN₁、DN₂为相邻的不同地物类型的像元遥感值，r₁、r₂为不同遥感像元对应地面的反射率。

而同一图斑内的大气程辐射遥感值是根据地物类型边界线上所有像元点的原始值与大气程辐射遥感值、地面反射率值之间的相关关系利用最小二乘法计算得到的。将图斑边界线上和同一图斑内的大气程辐射遥感值相加，可得到研究区的大气程辐射遥感遥感图像。

前述的基于像元分解的遥感图像大气程辐射获取方法中，所述植被三种类型的光谱端元根据时间分别按照春、夏、秋、冬的波谱测量时间进行分类。

与现有技术相比，本发明的方法完全基于遥感影像本身，从像元成像的角度去实现大气程辐射的逐点获取，从而可以用于实现城市市区范围内的大气气溶胶浓度的遥感监测。本发明的方法具有以下具体优势：

1、不需要相关的大气参数，可以应用在任何一个区域范围；

2、基于像元的计算尺度，可以应用在一个城市市区范围内，直观反映出城市的大气气溶胶污染空间分布差异；

3、研究区地物波谱库的建立，可以作为多年科研项目的积累，推广到其他科研项目应用。

4、本发明得到的大气程辐射的在蓝光波段效果尤为显著，最能反映大气气溶胶的空间分异特征。

5、本发明的地物反射率是基于光谱实测得到的，并可结合网络下载的波谱库进行归纳总结，便于推广应用。

6、本发明所使用的像元分解算法技术目前较为成熟，是准确获取地面反射率遥感图的关键。

7、本方法在用于中等尺度多波段遥感影像时效果更优，如ASTER、Landsat TM、环境卫星等。

附图说明

图1为电磁波辐射传输原理示意图

图2为光谱端元选择结果示意图

图3为用本发明实施例中得到的研究区地面反射率分布图；

图4为实施例中根据本发明方法得到的杭州市区大气程辐射遥感图像(蓝光波段为例)。

具体实施方式

下面结合实施例对本发明作进一步的说明，但并不作为对本发明限制的依据。

实施例。大气程辐射是瑞利散射和气溶胶米氏散射综合作用的结果，到达传感器的总的大气程辐射可以看成地表-传感器路径上大气上行辐射的积分(如图1所示)，程辐射的理论计算公式(7)。当电磁波进入大气层，大气层对电磁波有衰减作用(吸收作用、散射作用)，经地面反射，再次经大气层的衰减作用，最终到达传感器。传感器接收到的能量值——DN(单波段)表示为：

DN_ij＝K×E_ij×r_ij×τ_1ij×τ_2ij/π+K×Na_ij＝K×E_ij×r_ij×τ_ij/π+DNa_ij (7)

式(7)中，K为卫星传感器所固有的增益系数(可以从影像头文件获取，通常为Gains系数)；E_ij为像元的地面辐射照度；r_ij为像元的地面反射率；τ_1ij为像元的大气光谱下行透过率；τ_2ij为像元的大气光谱上行透过率；则τ_ij为像元的大气光谱上行透过率与下行透过率之积；Na_ij是大气程辐射值；DN_ij、DNa_ij分别为像元的遥感值和大气程辐射遥感值(汤定元等，1979；李先华，1993)。大气程辐射遥感值DNa_ij为大气亮度，其物理意义为大气总的反向散射，是仅与大气状况有关而与地面目标无关的大气向上散射，可以表征大气组分特征。从原始的遥感影像像元值中“分离”出大气程辐射遥感值，生成一幅新的遥感影像——大气程辐射遥感影像，以此来研究大气环境质量。大气程辐射的产生的物理意义是本专利申请的能产生社会效应的依据。

对于公式(7)，逐点计算DNa_ij，需要获取每个像元的E_ij、r_ij、τ_ij。本实施例的研究区定为杭州市，r_ij通过杭州典型地面波谱测量以及混合像元分解方法来获取。现在需要求出E_ij、τ_ij，按照常规计算方法，需要获取成像当天的气象条件，显然对计算要求较高，不适合技术推广。本次计算采用假设推算法。

假设在1km²范围内，大气状况均一，在此范围内，τ_ij、E_ij不变。因此我们可以将一幅30m*30m空间分辨率的影像，按照1km*1km大小重新划定网格，每个网格内可以逐点计算。其方法是：

对于1km²的网格内，可以计算两个相邻的像元点(A₁和A₂)的大气程辐射遥感值DNa_A1、DNa_A2。根据假设，该两点位置中，应存在以下等式：

E_A1＝E_A2＝E，DNa_A1＝DNa_A2＝DNa，τ_A1＝τ_A2＝τ (8)

根据公式(8)和公式(7)，可得：

DNa＝DN_A1-r_A1×(DN_A1-DN_A2)/(r_A1-r_A2) (9)

式(9)中，DN_A1、DN_A2为相邻的像元A₁、A₂的像元遥感值。从结果可以看出，若想等式成立或者能计算得到结果，r_A1、r_A2必须不相等，即A₁、A₂像元代表的地物必须是不同的地物类型，即同一地物类型的区域需要其他计算方法。对于城市区域而言，30m的像元尺度很少有同质像元，相邻的同质像元更少。因此我们在计算的时候，只计算相邻像元不相等情况下的DNA值；当相邻像元相等时，先不进行计算。

研究区的未计算的相邻像元、像元值相等的区域，计算方法如下。

根据大气辐射过程对像元点的原始值、大气程辐射遥感值、地面反射率值之间的相关关系进行理论推导。

大气辐射过程中传感器接受的信号可表示为:

L＝τL₀+DNa (10)

式中L为传感器接受的信号；DNa为大气呈辐射(太阳辐射通过大气分子与气溶胶的单次散射和多次散射)；L₀为直接太阳辐射和太阳辐射通过大气分子与气溶胶的单次散射及多次散射经地面目标反射到传感器方向的辐射；τ为大气光谱透过率。

根据图像构成，传感器输出的遥感图像的辐照度，即通常本领域所认为的DN值)可表示为：

DN＝HL+Z (11)

式中：H代表传感器对不同波段的输入响应，以便在不同的光照条件下得到合理的读数范围；Z为传感器的零输入响应。

在理想的漫反射(即目标为朗伯体)，目标上的辐照度为E,目标上的反射率为r则

对传感器的某一波段λ而言，存在下列关系式：

${\mathrm{DN}}_{\mathrm{\λ}}=\frac{{\mathrm{\τH}}_{\mathrm{\λ}}{E}_{\mathrm{\λ}}}{\mathrm{\π}}{r}_{\mathrm{\λ}}+H_{\mathrm{\λ}}{\mathrm{DNa}}_{\mathrm{\λ}}+Z_{\mathrm{\λ}}---\left(12\right)$

由式(12)可解出DNa_λ：

${\mathrm{DNa}}_{\mathrm{\λ}}=\frac{1}{H_{\mathrm{\λ}}}{\mathrm{DN}}_{\mathrm{\λ}}-\frac{{\mathrm{\τE}}_{\mathrm{\λ}}}{\mathrm{\π}}{r}_{\mathrm{\λ}}-\frac{1}{H_{\mathrm{\λ}}}{Z}_{\mathrm{\λ}}={\mathrm{aDN}}_{\mathrm{\λ}}+{\mathrm{br}}_{\mathrm{\λ}}+c---\left(13\right)$

其中： $a=\frac{1}{H_{\mathrm{\λ}}},b=-\frac{{\mathrm{\τE}}_{\mathrm{\λ}}}{\mathrm{\π}},c=-\frac{Z_{\mathrm{\λ}}}{H_{\mathrm{\λ}}}$

公式(13)表明某一波段像元的程辐射值与地物的DN值、反射率值呈回归关系。

图斑内部同质像元的遥感值计算可以利用图斑周边像元已知大气呈辐射遥感值、地物的DN值、反射率值，用最小二乘法来计算式(13)中回归系数a,b,c；然后利用图斑同质像元的地物的DN值和反射率值计算得到同质像元的大气程辐射值。

r_ij的计算方法。与原来依靠土地利用现状图的方法不同的是，本次依靠混合线性像元分解方法来获取。线性混合像元分解中，在任意一个波段上，任意一个像元的响应是该像元内各组分响应的线性和。因此，像元第i波段的反射率ρ可以表示为：

${\mathrm{\ρ}}_i=\underset{i=1}{\overset{n}{\Σ}}\left(a_{ij}{x}_j\right)+e_i---\left(14\right)$

式(14)中，ρ_i为包含一个或多个组分的像元在第i个波段上的平均光谱反射率；α_ij为像元内第i个波段上第j个组分的反射率；x_j为像元内第j个组分所占的比例；e_i为第i个光谱波段的误差项；假设一像元内有n个组分(j＝1,2,,n)，传感器系统的波段数为m(i＝1,2,,m)。此处需对方程加以约束条件，即：且0≤x_j≤1 (15)

$RMS=\sqrt{\left(\underset{i=1}{\overset{m}{\Σ}}{e}_i^2\right)/m}---\left(16\right)$

当RMS满足规定要求，光谱分解结果可信。依靠MNF(MinimumNoise Fraction Transform)方法，进行光谱端元选择，实现线性光谱分解。根据城市地物特征，本次专利选择的光谱端元为植被端元、高照度不透水面端元(高反照率)、低照度不透水面端元(低反照率)，具体的端元选择过程如图2所示。

(10)式中的a_ij，通过实测研究区地物波谱，建立城市典型地物波谱库，并在光谱端元概况的基础上获取。表1为本次研究的三个端元在TM8的各个波段上的地物反射率概括值(a_ij值)。

表1基于Landsat 8光谱端元的地物反射率概括值

有地理坐标的大气程辐射遥感影像的最终生成。按照上述方法生成的大气程辐射遥感值，只是矩阵形式，而在空间分析过程中，需要对影像赋予地理坐标。参与计算的遥感影像与地面反射率图均是经过精纠正，有相同的地理坐标系统和投影信息，因此生成的大气程辐射遥感影像只需要对矩阵(ascii文件格式)赋予有地理坐标描述和投影信息的头文件，便可使用ARCGIS、ERDAS、ENVI等专业软件，转为有地理坐标和投影信息的影像文件。该地理坐标描述和投影，这是将研究结果进行空间分析的必要前提。

综合上所述，通过实测地物波谱，在混合线性分解模型基础上，概况出端元反射率值，便可逐点计算反演大气程辐射遥感值，从而最终生成有精确坐标的大气程辐射遥感影像。

本次实施例中杭州市Landsat TM8第二波段(蓝光波段)反射率遥感图如图3所示。而杭州市Landsat TM8第二波段(蓝光波段)大气程辐射遥感图如图4所示。

Claims (3)

1.一种基于像元分解的遥感图像大气程辐射获取方法，其特征在于,包括以下步骤：

①对原始遥感影像进行几何校正；

②将采集的地面波谱数据进行归纳为低照度不透水面、高照度不透水面和植被三种类型的光谱端元，并在遥感影像的波段区间获得各光谱端元在波段区间内的平均反射率；

③基于像元分解方法，获取原始遥感影像中植被、高照度不透水面、低照度不透水面三个光谱端元比例分布图；

④结合步骤②和③，得到原始遥感影像的地面反射率图像；

⑤在遥感成像的原理下，根据步骤④得到的地面反射率图像，对像元进行逐点分解，得到每个像元的大气程辐射遥感值；

在进行步骤⑤的大气程辐射遥感值的计算前，根据步骤④获得的地面反射率，划定不同地物类型的边界阈值，获得原始遥感影像所示的研究区的地物类型边界线，从而将整幅图分为同一图斑和不同图斑边界线两个区域分别计算大气程辐射遥感值；同一图斑内的大气程辐射遥感值是根据地物类型边界线上所有像元点的原始值与大气程辐射遥感值、地面反射率值之间的相关关系利用最小二乘法计算得到的。

2.根据权利要求1所述的基于像元分解的遥感图像大气程辐射获取方法，其特征在于：图斑边界线上的大气程辐射遥感值根据公式DNa＝DN₁-r₁×(DN₁-DN₂)/(r₁-r₂)进行计算，其中DN₁、DN₂为相邻的不同地物类型的像元遥感值，r₁、r₂为不同遥感像元对应地面的反射率。

3.根据权利要求1所述的基于像元分解的遥感图像大气程辐射获取方法，其特征在于：所述低照度不透水面、高照度不透水面和植被这三种类型的光谱端元根据时间分别按照春、夏、秋、冬的波谱测量时间进行分类。