专利摘要

本发明公开了气溶胶和地表参数联合反演方法，通过所提取的两个光谱基础向量来构造0.45~0.52μm、0.52~0.59μm、0.63~0.69μm和0.77~0.89μm这4个波段的地表反射率，平均绝对误差小于0.01。本发明气溶胶和地表参数联合反演方法通过卫星遥感多光谱地表反射率的重建方法，可有效实现不同地物的多光谱地表反射率的精确重建，针对GF‑1卫星遥感数据的地溶胶和地表参数联合反演方法。该方法不需要事先准备好的地表反射率，也不需要利用时间序列的卫星数据来约束反演，只需要单景的GF‑1卫星数据即得到对应的气溶胶光学厚度和地表反射率等结果。

权利要求

1.一种气溶胶和地表参数联合反演方法，其特征在于，包括以下步骤：

1)根据公式(5)获得无云像元近似的地表反射率；

其中L(λ)表示GF-1卫星实测的归一化辐亮度值，LRay(λ)为从瑞利散射查找表中得到的归一化的瑞利散射辐亮度结果,μ0为太阳天顶角θ0的余弦值；

2)根据公式(6)获得近似的多光谱地表反射率权重系数；

和为近似的多光谱地表反射率权重系数；

3)以公式(7)为价值函数，进行大气校正仪标量观测的地气解耦和；

其中下标j表示所采用GF-1卫星多光谱观测的4个波段，上标i表示对应的参数组合的标号，Lmeas和Lsimu分别为实际观测和辐射传输标量模拟的归一化辐亮度结果，参数向量为所设的气溶胶微物理参数V0和FMFV和地表多光谱重建权重系数a1和a2的网格点，进而可以计算对应的价值函数并进行从小到大排序；通过遍历查找中先进行太阳天顶角、观测天顶角和相对方位角观测几何的匹配，最终查找到minJ对应的和即为反演得到的V0和FMFV，对应的和即为多光谱地表反射率重建的权重系数组合；a1下限为上限为a2的下限为上限为

4)基于所反演得到的气溶胶柱体积浓度的结果，利用所建立的气溶胶微物理参数V0和FMFV到光学参数转换的查找表，索引和插值得到对应的光谱气溶胶光学厚度τa(λ)结果；根据公式(8)得到反演结果；

τa(λ)＝LUT(V0,FMFV)(8)

其中τa(λ)表示对应的AOD。

说明书

技术领域

本发明涉及大气污染治理领域，尤其是一种大气气溶胶观测技术。

背景技术

大气中的气溶胶对全球辐射平衡、气候变化、人类健康等起着直接或间接的作用。气溶胶光学厚度作为气溶胶重要的光学参数，表示气溶胶粒子的消光能力，被广泛应用于各个相关领域。目前气溶胶的常规监测主要有地基、星载两种方式，地基观测具有精度高的特点。相比于地基观测，星载监测气溶胶具有范围广、成本低、受地形条件约束等特点，但是基于卫星传感器反演气溶胶难度较大，另外传统的基于卫星传感器数据反演的气溶胶光学厚度产品存在空间分辨率低的缺点，例如POLarization and Directionality of theEarth’s Reflectances (POLDER)的官方二级产品空间分辨率为18km，ModerateResolution Imaging Spectroradiometer(MODIS)的官方产品初始空间分辨率为10km。这些产品远远不能满足如城市等小尺度区域的气溶胶研究。

为了发展高空间分辨率气溶胶光学厚度产品，已有大量学者展开了相关的研究，并取得了一些成果，例如王中挺等基于高分一号传感器数据，利用经典的暗目标方法中红蓝波段的线性关系反演气溶胶光学厚度，但是卫星反演结果整体偏高，且该方法只适用于植被等暗地表地区。Bao等对经典暗目标法中红蓝波段关系进行了改进，统计了不同的Normalized Difference Vegetation Index(NDVI)阈值区间的红蓝波段地表反射率线性关系，并将此方法应用于高分一号传感器数据反演气溶胶光学厚度，但是反演结果还是存在整体偏高的问题。 Sun等利用深蓝算法，并基于高分一号传感器数据反演气溶胶光学厚度，但是反演结果还是存在整体偏高的问题，而且该方法需要基于长时间序列数据制备地表反射率库。She等利用 MODIS地表反射率产品，统计可见光与近红外波段地表反射率线性关系，并将该方法应用高分四号传感器数据反演气溶胶光学厚度，并与地基站点结果验证，相关系数达到0.92，但是该线性关系不适用于某些地表。Chen等基于图像时序的方法，利用“大气变化快，地表变化慢”的特点来对高分四号卫星传感器数据进行地气解耦合反演气溶胶光学厚度，反演了800m 空间分辨率的气溶胶光学厚度产品，其反演结果与地基站点验证显示相关性为R＝0.79。该方法类似于深蓝算法，需要事先制备地表反射率库，但是其可以根据反演的地表反射率更新地表反射率库。

由上述描述可见，目前基于高空间分辨率卫星遥感数据反演气溶胶光学厚度还存在各种问题，例如暗目标法只能应用于植被等暗地表地区，且波段间的线性关系对于特定地表不适合；地表反射率库方法需要事先制备地表反射率库，对于沙漠等地表变化不大的地区比较合适，但是对于植被等地表反射率变化大的地区往往是不适合的。图像序列方法虽然可以更新地表反射率库，但是基于自身反演的地表反射率，这样会导致误差累积增大，另外也存在图像配准不准等问题。

发明内容

本发明的目的是提供一种针对GF-1卫星遥感数据的地溶胶和地表参数联合反演方法。该方法不需要事先准备好的地表反射率，也不需要利用时间序列的卫星数据来约束反演，只需要单景的GF-1卫星数据即得到对应的气溶胶光学厚度和地表反射率等结果；具体技术方案为：

一种气溶胶和地表参数联合反演方法，包括以下步骤：

1)根据公式(5)获得无云像元近似的地表反射率；

其中L(λ)表示GF-1卫星实测的归一化辐亮度值，LRay(λ)为从瑞利散射查找表中得到的归一化的瑞利散射辐亮度结果,μ0为太阳天顶角θ0的余弦值；

2)根据公式(6)获得近似的多光谱地表反射率权重系数；

和 为近似的多光谱地表反射率权重系数；

3)根据公式(7)进行大气校正仪标量观测的地气解耦和；

其中下标j表示所采用GF-1卫星多光谱观测的4个波段，上标i表示对应的参数组合的标号， Lmeas和Lsimu分别为实际观测和辐射传输标量模拟的归一化辐亮度结果，a1和a2为多光谱地表反射率权重系数；

a1下限为 上限为 a2的下限为 上限为

4)根据公式(8)得到反演结果；

τa(λ)＝LUT(V0,FMFV)(8)

其中τa(λ)表示对应的AOD。

本发明气溶胶和地表参数联合反演方法通过卫星遥感多光谱地表反射率的重建方法，可有效实现不同地物的多光谱地表反射率的精确重建，针对GF-1卫星遥感数据的地溶胶和地表参数联合反演方法。该方法不需要事先准备好的地表反射率，也不需要利用时间序列的卫星数据来约束反演，只需要单景的GF-1卫星数据即得到对应的气溶胶光学厚度和地表反射率等结果。

附图说明

图1为对典型的地物光谱库4波段光谱重建结果散射图；

图2为本发明气溶胶和地表参数联合反演方法实施例1与现有技术的结果对比图；

图3为本发明气溶胶和地表参数联合反演方法实施例2与现有技术的结果对比图；

图4为本发明气溶胶和地表参数联合反演方法实施例3与现有技术的结果对比图；

图5为本发明气溶胶和地表参数联合反演方法实施例4与现有技术的结果对比图；

图6为本发明气溶胶和地表参数联合反演方法实施例5与现有技术的结果对比图；

图7为本发明气溶胶和地表参数联合反演方法实施例1-5与现有技术的数据结果对比图。

具体实施方式

下面利用实施例对本发明进行更全面的说明。本发明可以体现为多种不同形式，并不应理解为局限于这里叙述的示例性实施例。

GF-1卫星有效载荷技术指标如表4.1所示，包含全色多光谱相机和多光谱相机两个载荷。本实施例综合采用0.45～0.52μm、0.52～0.59μm、0.63～0.69μm和0.77～0.89μm在这4个多光谱谱段进行气溶胶和地表参数的联合反演。

数据处理的步骤如下：

(1)定标处理：从遥感影像像元亮度值DN值计算得到辐亮度值。

(2)重采样处理：将像元空间分辨率由8m或16m重采样为1km，减少观测噪声。

(3)读取观测几何信息：读从xml文件中读取对应的太阳天顶角、太阳方位角、观测天顶角和观测方位角信息。

(4)输出为txt文件，用于后续的反演处理。

表1 GF-1卫星有效载荷技术指标

对应GF-1卫星0.45～0.52μm、0.52～0.59μm、0.63～0.69μm和0.77～0.89μm这4多光谱波段，采用的地表反射率重建公式为

其中，向量 和 为4个波段对应的2个光谱基础向量；a1和a2为对应的线性组合权重系数，a1≥0，a2≥0， 也就是说，任意的多光谱地表反射率均可通过调整权重系数a1和a2的值来近似得到。在本方法中，我们通过对两个光谱基础向量权重系数a1和a2的反演来代替对4个波段地表反射率的反演，减少了需要反演的参数。

当前国际上比较常用的地物光谱库包括USGS光谱库和ASTER光谱库。其中，USGS光谱库(http://speclab.cr.usgs.gov/spectral.lib06/)是美国地质勘探局(UnitedStates Geological Survey，USGS)光谱实验室所建立的波长范围在0.2～3.0μm之间的光谱库，主要的光谱可分为矿物、混合物、涂料、挥发物、人造材料和植被等6大类。ASTER光谱库 (http://speclib.jpl.nasa.gov)是美国宇航局先进的空间热辐射反射测量仪(Advanced Space-borne Thermal Emission and Reflection Radiometer,ASTER)研究计划的一部分，主要包括矿物、岩石、土壤、人造材料、水和冰雪等物质的光谱，许多光谱的波长覆盖范围可达 0.4～14μm。

参照GF-1卫星多光谱波段设置，通过对这两个典型的光谱库中的光谱数据进行重采样，得到对应0.45～0.52μm、0.52～0.59μm、0.63～0.69μm和0.77～0.89μm的地表反射率光谱数据集，对地表反射率多光谱重建效果进行验证(包括植被、裸土、牧场、混凝土、岩石、水体及混合地表等)，如图1所示，光谱重建的平均绝对误差小于0.01。

根据AERONER和SONET地基实测数据及文献资料，对中国典型区域的气溶胶微物理参数进行统计，包括粗细模态的体积柱浓度、粒子有效半径、有效方差及对应的复折射指数等，建立典型区域的动态气溶胶模型，用于辐射传输计算和对应查找表的建立。其中，所建立的典型气溶胶模型气溶胶粒子谱分布的形式为

其中V表示体积柱浓度随粒子半径r的函数，r表示粒子半径，V0表示具体的体积柱浓度，rV表示体积谱对应的几何中值半径(单位μm)，σg为对应的方差，上标f和c分别表示细模态和粗模态的气溶胶。

本发明所采用的辐射传输模式为国际上近几年发布的UNL-VRTM(UnifiedLinearized- Vector Radiative Transfer Model)。UNL-VRTM模式在VLIDORT(VectorLinearized Discrete Ordinate Radiative Transfer)的基础上进行开发，增加了线性化的Mie散射、瑞利散射、大气廓线、气体吸收和方便的参数输入输出等模块，是一个开源易用的矢量辐射传输模型 (www.unl-vrtm.org)。

对应所设定的气溶胶模型，利用UNL-VRTM模式建立气溶胶微物理参数与光学参数之间的转换关系，并建立气溶胶微物理与光学参数查找表，以及索引和插值规范。在确定选用的体积柱浓度V0和体积柱浓度细模态占比FMFV的基础上，对应的细模态和粗模态的体积柱浓度为

可进一步得到光谱的气溶胶光学厚度(AOD)τa(λ)的形式为

其中，τa(λ)表示对应的AOD，reff为气溶胶粒子有效半径，Qext为对应的消光效率因子，可由辐射传输模式直接计算得到，上标f和c分别表示细模态和粗模态的气溶胶。

通过设定陆地上空对应的多光谱波段、观测几何、气溶胶模型参数、地表反射率参数，进行辐射传输计算，模拟得到对应的GF-1卫星多光谱观测结果，并建立对应的查找表。

表2利用辐射传输模式所建立查找表的参数范围

载入所建立的GF-1卫星多光谱观测查找表，基于标量查找表的维度规则，根据对应的太阳天顶角θ0、观测天顶角θv和相对方位角φ对原始的查找表进行插值和索引，建立观测几何 {θ0,θv,φ}对应的包括体积柱浓度V0、体积柱浓度细模态占比FMFV、地表反射率和归一化表观辐亮度L维度的子查找表。根据包括体积柱浓度V0体积柱浓度细模态占比FMFV所对应的AOD 结果，对所建立子查找表按照更为细化的AOD间隔(0.025)进行插值，得到插值后的子查找表，用于后续的地气解耦和反演处理。

表3对应观测几何{θ0,θv,φ}索引和插值得到的子查找表的参数范围

通过设定对应的观测几何和地面高程(对应气压)，进行瑞利散射在多光谱波段的散射贡献，建立瑞利散射查找表，用于后续的瑞利散射校正，以获取同步的地表反射率先验知识。

表4利用辐射传输模式所建立瑞利散射查找表的参数范围

基于GF-1卫星实际观测的多光谱辐亮度信息，通过瑞利散射校正，获取近似的多光谱地表形状曲线。瑞利校正后得到的近似的地表反射率 的形式为

其中L(λ)表示GF-1卫星实测的归一化辐亮度值(即由卫星遥感数据提供)，LRay(λ)为从瑞利散射查找表中得到的归一化的瑞利散射辐亮度结果,μ0为太阳天顶角θ0的余弦值。在此基础上，获取近似的多光谱地表反射率权重系数 和 对应的计算公式为

那么，就可以在 和 的上下边界范围内搜索实际GF-1卫星观测的多光谱地表反射率a1和a2的值，有效地缩小了搜索范围。

表5多光谱地表反射率光谱重建的权重系数范围

基于所提取的多光谱地表反射率先验知识和约束范围，对地表反射率进行多光谱重建，并根据观测几何所索引和插值得到的结果，建立包含气溶胶参数和权重系数的子查找表，通过参数遍历来计算标量观测的代价函数结果，进行多光谱观测的地气解耦。

为了进行大气校正仪标量观测的地气解耦和，我们定义对应的价值函数J的形式为

其中下标j表示所采用GF-1卫星多光谱观测的4个波段，上标i表示对应的参数组合的标号， Lmeas和Lsimu分别为实际观测和辐射传输标量模拟的归一化辐亮度结果。参数向量 为所设的气溶胶微物理参数(V0和FMFV)和地表多光谱重建权重系数(a1和 a2)的网格点，进而可以计算对应的价值函数并进行从小到大排序。

通过遍历查找中先进行太阳天顶角、观测天顶角和相对方位角等观测几何的匹配，最终查找到minJ对应的 和 即为反演得到的V0和FMFV，对应的 和 即为多光谱地表反射率重建的权重系数组合。

基于所反演得到的气溶胶柱体积浓度的结果，利用所建立的气溶胶微物理参数V0和FMFV到光学参数转换的查找表，索引和插值得到对应的光谱气溶胶光学厚度τa(λ)结果，对应的转化公式为

τa(λ)＝LUT(V0,FMFV) (8)

进而得到对应的气溶胶光学厚度反演结果。

实施例1：复杂地表低气溶胶情况下的高空间分辨率AOD精确反演

图2(a)为甘肃民勤2014年2月19日GF-1彩色合成图，(b)为本发明反演得到1km 空间分辨率的气溶胶光学厚度(AOD)结果空间分布图。该研究区域的地表类型非常复杂，包括沙漠、地山丘陵和平原(城市、植被、裸土、小型水体)等多种地表，从本发明反演的 AOD结果图上可以看出，在去除少部分超亮地表的基础上，气溶胶光学空间分布比较平滑，其与地观测结果的550nm AOD结果(AOD＝0.12)具有很好的一致性，反演精度在复杂亮地表和低气溶胶情况依然可为满足±(20％*AOD(地基实测)+0.05)的精度。

实施例2：复杂地表轻度污染情况下的高空间分辨率AOD精确反演

图3(a)为陕西西安2015年5月6日GF-1彩色合成图，(b)为本发明反演得到1km空间分辨率的AOD空间分布图。从GF-1卫星彩色合成图上可以看出，该研究区域包括城市、河流和山地等复杂地表类型，本发明反演的AOD空间分布明显呈现从人类生活密集的城市到活动较少的郊区扩散递减趋势，1km的AOD空间分辨率相比MODIS两种算法10km空间分辨率具有明显的高空间分辨率优势。同时，在市区轻度污染的情况下，反演可与地基观测结果(AOD＝0.44)比较一致，反演误差在(20％*AOD(地基实测)+0.05)精度之内。

实施例3：复杂地表轻度污染情况下的高空间分辨率AOD精确反演

图4(a)为陕西西安2014年12月30日GF-1彩色合成图，(b)为本发明反演得到1km空间分辨率的AOD空间分布图。从本发明反演的AOD空间分布图上可以看出，气溶胶光学厚度的分布同样从城市到郊区扩散递减趋势，且反演的AOD具有高空间分辨率的优势。在市区轻度污染的情况下，反演可与地基观测结果(AOD＝0.39)保持一致，反演误差小于 (20％*AOD(地基实测)+0.05)的精度。

实施例4：复杂地表中度污染情况下的高空间分辨率AOD精确反演

图5(a)为中度污染情况下陕西西安2015年1月3日GF-1彩色合成图，(b)为本发明反演得到1km空间分辨率的AOD空间分布图。从本发明反演的AOD空间分布图上可以看出，AOD的分布从城市中心向四周扩散递减，且展示了AOD的在城市的高空间分辨率分布细节。在市区中度污染的情况下，反演结果与地基观测结果(AOD＝0.77)较为吻合，反演误差小于(20％*AOD(地基实测)+0.05)的精度。

实施例5：复杂地表中度污染情况下的高空间分辨率AOD精确反演

图6(a)为中度污染情况下陕西西安2015年2月1日GF-1彩色合成图，(b)为本发明反演得到1km空间分辨率的气溶胶光学厚度(AOD)结果空间分布图。从本发明反演的 AOD空间分布图上可以看出，AOD的分布同样从城市中心向周边的郊区和山区扩散递减， AOD反演结果可与地基观测结果(AOD＝0.58)基本相当，误差在(20％*AOD(地基实测)+0.05) 的精度范围内。

图7(a)为对应实施例1至实施例5的550nm AOD反演结果，本发明基于GF-1反演的AOD结果与地基观测网SONET实测的AOD结果对比散点图，其中拟合线为y＝0.85x+0.03，相关系数科大R2＝0.93，且AOD误差均在(20％*AOD(地基实测)+0.05)的精度范围内；(b) 为相近时刻MODIS深蓝(DB)算法反演的AOD结果与地基观测网SONET实测的AOD结果对比散点图，其中拟合线为y＝0.46x+0.17，相关系数R2＝0.86。

通过以上对比验证可以看出，本发明得到的AOD结果空间分辨率和反演精度均优于 MODIS卫星AOD官方产品，且可适用与各种复杂亮地表的情况。

上述示例只是用于说明本发明，除此之外，还有多种不同的实施方式，而这些实施方式都是本领域技术人员在领悟本发明思想后能够想到的，故，在此不再一一列举。

一种气溶胶和地表参数联合反演方法专利购买费用说明

Q：办理专利转让的流程及所需资料

A：专利权人变更需要办理著录项目变更手续，有代理机构的，变更手续应当由代理机构办理。

1：专利变更应当使用专利局统一制作的“著录项目变更申报书”提出。

2：按规定缴纳著录项目变更手续费。

3：同时提交相关证明文件原件。

4：专利权转移的，变更后的专利权人委托新专利代理机构的，应当提交变更后的全体专利申请人签字或者盖章的委托书。

Q:专利著录项目变更费用如何缴交

A：（1）直接到国家知识产权局受理大厅收费窗口缴纳,（2）通过代办处缴纳，（3）通过邮局或者银行汇款,更多缴纳方式

Q：专利转让变更，多久能出结果

A：著录项目变更请求书递交后，一般1-2个月左右就会收到通知，国家知识产权局会下达《转让手续合格通知书》。