一种雾霾条件下地表太阳短波辐射计算方法与流程

本申请涉及卫星遥感反演领域，尤其涉及一种雾霾条件下地表太阳短波辐射计算方法。

背景技术：

太阳辐射提供了地球各个圈层运行和生产的能量。而地表辐射收支过程与能量循环是地表过程的重要组成部分，是地区系统演化发展的重要驱动力，也是影响气候变化和地气相互作用的关键过程。太阳辐射穿透大气到达地表，需要经过大气的吸收和散射过程。

大气气溶胶(雾霾)是悬浮在大气中的固体颗粒和液体微粒与气体载体共同组成的多相体系，通常存在于对流层。气溶胶对太阳辐射的散射和吸收强烈影响地球辐射收支平衡，研究雾霾条件下地表太阳短波辐射对正确地评估气溶胶的气候效应有十分重要的作用。

所述地表太阳短波辐射计算方法包含：经验统计模型法、辐射传输方法、参数化方法和查找表方法等。其中，所述参数化方法(本申请文件中或称参数法)是基于大气辐射传输模型、模拟数据或机载、星载及观测数据等，建立关键大气、地表变量与下行短波辐射的参数关系式。所述查找表方法基于原理是大气光学特性和地物物理特性决定了大气顶太阳短波辐射和地表太阳短波辐射的变化，因此特定大气和地表参数会对应唯一的地表太阳短波辐射，所以可以预先设定不同的大气和地表参数，然后采用辐射传输模型计算输入参数情形下对应的地表太阳短波辐射，得到查找表，最后就可以根据不同输入图像的具体情况选择不同的查找表得到地表太阳短波辐射值数据结果。

经验统计方法和参数化方法基于局地观测资料或经验系数，精度高但普适性较差；辐射传输方法计算精度高但计算时间过长；而查找表方法由于输入参数和输出参数已经预先计算好，这样可以节约大量时间，但是计算精度和准确性不够。

技术实现要素：

本申请提出一种雾霾条件下地表太阳短波辐射计算方法，解决了现有技术只采用单一方法反演地表太阳短波辐射，无法提高精度和计算速度的问题。

本申请实施例提出一种雾霾条件下地表太阳短波辐射计算方法，包含以下步骤：

提取卫星的观测像元，去掉有云的像元图像；

用第一查找表计算地表太阳短波辐射，作为第一辐射值；所述第一查找表，包含太阳天顶角、气溶胶光学厚度、气溶胶类型和对应的地表太阳短波辐射的数据；

读取卫星的地表太阳短波辐射数据，为第二辐射值；与所述第一辐射值对比；若误差大于第一阈值，则该像元为校正像元；所述第二辐射值，为卫星产品自带的地表太阳短波辐射数据。

对所述校正像元用参数化方法计算地表太阳短波辐射，为第三辐射值。

本申请优化的实施例，用第一查找表之前，还包含以下步骤：计算像元的地表反射率，用第二查找表计算出气溶胶厚度和气溶胶类型；所述第二查找表，包含表观反射率均值、地表反射率和对应的气溶胶厚度、气溶胶类型。

本申请优化的实施例，在对校正像元用参数化方法计算第三辐射值后，还包含以下步骤：

提取校正像元的第三辐射值和第一辐射值进行对比得到调整值，用所述调整值对第一辐射值数据进行校准。

优选地，在本申请任意一项实施例中，对第一观测时刻的校正像元，计算第二观测时刻的第一辐射值，如果误差仍大于第一阈值，则以第一观测时刻所述校正像元的第三辐射值为基础，根据第一查找表，得出更新的太阳天顶角、气溶胶厚度和气溶胶类型；第一观测时刻和第二观测时刻之间的间隔小于第二阈值。

优选地，在本申请任意一项实施例中，对第一观测时刻的校正像元，计算第二观测时刻的第一辐射值，如果误差仍大于第一阈值，则以第一观测时刻所述校正像元的第三辐射值作为第二观测时刻所述校正像元的地表太阳短波辐射；第一观测时刻和第二观测时刻之间的间隔小于第二阈值。

优选地，通过静止卫星葵花-8采集观测数据。

优选地，所述参数化方法采用era-interim数据，参数化方法的计算参数包含：地表反射率、大气压、总臭氧含量、可降水量、气溶胶厚度和气溶胶类型。

优选地，所述第一查找表是基于大气辐射传输模式rstar构建的查找表。

进一步地，所述第一阈值为30w/m2。

优选地，所述太阳天顶角的参数范围为1-90度，间隔1度；

所述气溶胶光学厚度的参数范围为0-3，间隔为0.5。

优选地，通过葵花-8的二级云识别产品判断观测区域是晴空还是有云。

本申请实施例采用的上述至少一个技术方案能够达到以下有益效果：

用查找表法计算地表太阳短波辐射，用参数化方法校正结果，计算的速度快，精度高。本申请的实施例中，可采用静止轨道卫星手段进行观测，观测范围可覆盖全球，相比地基站点观测范围广，相比极轨卫星可以对区域进行持续监测，有更高时间分辨率，进一步提高速度和精度。本专利基于静止卫星高时空分辨率全球数据，数据量巨大，需要算法计算时间短、普适性高。

附图说明

此处所说明的附图用来提供对本申请的进一步理解，构成本申请的一部分，本申请的示意性实施例及其说明用于解释本申请，并不构成对本申请的不当限定。在附图中：

图1为本发明地表太阳短波辐射计算方法实施例流程图；

图2为本发明包含气溶胶信息计算步骤的方法实施例流程图；

图3为本发明包含查表法校准步骤的实施例流程图；

图4为本发明包含第二观测时刻气溶胶参数校正的实施例流程图；

图5为本发明包含第二观测时刻第三辐射值替代的实施例流程图。

具体实施方式

为使本申请的目的、技术方案和优点更加清楚，下面将结合本申请具体实施例及相应的附图对本申请技术方案进行清楚、完整地描述。显然，所描述的实施例仅是本申请一部分实施例，而不是全部的实施例。基于本申请中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本申请保护的范围。

以下结合附图，详细说明本申请各实施例提供的技术方案。

图1为本发明地表太阳短波辐射计算方法实施例流程图。

本申请实施例提供一种地表太阳短波辐射计算方法，包含以下步骤：

步骤101、提取卫星的观测像元，去掉有云的像元图像。

在步骤101中，提取卫星的观测像元点上全部像元图像，判断像元图像晴空还是有云，去掉有云的像元图像。所述观测像元点上的像元图像是通过静止卫星采集的观测区域的观测数据。

所述像元为卫星观测图像上的一个像素点。所述像元点为同一坐标上所有像元的集合。

目前地表太阳短波辐射测量仍是以地面观测为主，由于地基站点观测覆盖范围有限，并不适合大尺度测量。而采用静止卫星测量地表太阳短波辐射，则可获得大面积、连续高精度观测资料，提供区域尺度上地表能量较为详细的空间分布信息。

例如，所述静止卫星可以选择葵花-8，也可以选择风云-4，还可以选择其他静止卫星。

优选的，所述静止卫星为葵花-8，所述葵花-8是日本葵花系列静止气象卫星中的第8颗，隶属于日本气象局(jma)。葵花8搭载了ahi(advancedhimawariimager)多光谱传感器，拥有16个通道，10分钟的时间分辨率。由于其有多光谱、高空间分辨率优势，其对全球云遥感、气候变化、地表太阳辐射研究有重要意义。

需要说明的是，本申请实施例提供的方法只适用于晴空状态，若提取的观测像元为有云，则无法进行后续的计算。

例如，可以选用静止卫星附带的云识别产品来判断观测像元为晴空还是有云。优选地，通过葵花-8的二级云识别产品判断观测区域是晴空还是有云。

步骤102、建立第一查找表，利用查找表法计算出地表太阳短波辐射数据。

用第一查找表计算地表太阳短波辐射，作为第一辐射值；所述第一查找表，包含太阳天顶角、气溶胶光学厚度、气溶胶类型和对应的地表太阳短波辐射的数据。

所述第一辐射值数据是通过查找表法计算出的地表太阳短波辐射数据。通过将太阳天顶角、气溶胶光学厚度和气溶胶类型参数输入第一查找表，即可算出第一辐射值数据。

所述第一查找表参数包含：太阳天顶角a、气溶胶光学厚度b和气溶胶类型c和地表太阳短波辐射数据d。随着数据a、数据b、数据c的变化，会得到不同的数据d。知道数据d可以反过来查找一组数据a、数据b、数据c。

例如，所述太阳天顶角的参数范围为1-90度，间隔1度。所述气溶胶光学厚度的参数范围为0-3km，间隔为km。

所述查找表方法，基于大气辐射传输模式rstar构建查找表。所述大气辐射传输模式rstar是1986年由日本nakajima和tanaka提出的一套适用于平面平行大气辐射传输数值模型，该模型由6种常见的大气模式可选，设置高度可从海平面到120km分50层，模式中提供了，预设由9种基本成分组成的11种气溶胶类型(water,ice，dust_like，soot，volcanic_ash，75％h2so4，rural，seaspray，urban，tropo，yellowsand)。

步骤103、读取卫星的地表太阳短波辐射数据，为第二辐射值；与所述第一辐射值对比；若误差大于第一阈值，则该像元为校正像元。

所述第二辐射值，为卫星产品自带的地表太阳短波辐射数据。

需要说明的是，此处第一辐射值是计算值，第二辐射值是测试值，当第一辐射值和第二辐射值的差别较大时，认为误差较大。误差产生的原因可以是计算和或测试过程所引入，这里不做具体分析。所述第一阈值，例如可以为30w/m²。

例如，所述第二辐射值数据，为静止卫星产品自带的地表太阳短波辐射数据，作为测试值，由于精度很低，一般不作单独使用。此处，用均方根来分析测试值和计算值之间的一致性。

例如，定义的误差＞第一阈值，则说明此处的像元图像用查找表法计算的第一辐射值数据与测试值的差别过大，则使用参数化方法计算。

计算误差的公式举例为：

式中，rmse为均方根误差；

n：葵花8拍摄图像的像元总数，n＝5764801；

i：像元编号；

swri,rstar：编号为i的像元，采用本专利方法，由改进大气辐射传输模式rstar得到的第一查找表，计算得到的第一辐射值；

swri,h8_l3：编号为i的像元，来自葵花8卫星官方l3地表太阳短波辐射数据，即第二辐射值。

步骤104、对所述校正像元用参数化方法计算地表太阳短波辐射，为第三辐射值。

例如，所述参数化方法采用era-interim数据，参数化方法的计算参数包含：地表反射率、大气压、总臭氧含量、可降水量、气溶胶厚度和气溶胶类型。

所述era-interim是由ecmwf(欧洲中期天气预报中心)开发的一款全球大气再分析资料。时间序列为1979年至今，空间分辨率0.703°×0.702°，时间分辨率6小时和12小时，产品通常每月更新一次。

参数化方法的计算公式为：

e＝e0t(2)

其中e为地表太阳短波辐射；e0为大气顶层入射短波辐射；t为整层大气总透过率，包含臭氧吸收、rayleigh散射、水汽吸收、混合气体吸收、气溶胶吸收和散射。

图2为本发明包含气溶胶信息计算步骤的方法实施例流程图。

作为本申请进一步优化的实施例，包含步骤101、101a、102、103、104，在构建第一查找表(102)之前，还包含步骤101a。具体说明如下：

步骤101、提取卫星的观测像元，去掉有云的像元图像；

步骤101a、计算像元的地表反射率，用第二查找表计算出气溶胶厚度和气溶胶类型；所述第二查找表，包含表观反射率均值、地表反射率和对应的气溶胶厚度、气溶胶类型。

要计算出气溶胶厚度和气溶胶类型，需要计算出卫星的观测像元的地表反射率。现有技术中，通过一个卫星观测像元的全部卫星的观测通道参数，就可以计算出这个像元图像的地表反射率。

例如，葵花-8的观测通道包含ir、ir1、ir2、ir3、ir4、band1～16，经过大量实验证明，计算地表反射率只需要读取晴空像元图像的band1、band3、band6三个通道数据。

首先要读取晴空的像元图像的band1、band3、band6三个通道数据，得到像元的表观反射率ρ0.46、ρ0.64、ρ2.3。

对ρ0.46、ρ0.64、ρ2.3进行水汽、o3和co2的吸收订正，若满足0.01≤ρ2.3≤0.25，则被标识为暗像元，将像元图像内所有挑选出的暗像元按照红波段的反射率ρ0.64排序，并将其中的最亮的50％和最暗的20％暗像元剔除，这一处理目的是减少过亮或过暗的异常像元干扰。剩余的暗像元取其表观反射率的均值，作为该子区域的表观反射率。

地表反射率计算流程为基于上述暗像元在0.46μm、0.64μm与2.3μm地表反射率之间的经验关系，利用ρ2.13求取地表反射率。方程是：

在步骤101a中，第二查找表为气溶胶查找表，第二查找表包含表观反射率均值、地表反射率和对应的气溶胶厚度、气溶胶类型，反射率为需要输入4个参数，分别是表观反射率均值ρ0.46、ρ0.64；地表反射率通过将地表反射率等参数对照第二查找表，即可计算出气溶胶光学厚度和气溶胶类型。

需要说明的是，所述气溶胶类型，分为城市型气溶胶、乡村型气溶胶和海洋性气溶胶。不同的气溶胶类型有不同的参数集。

步骤102、用第一查找表计算地表太阳短波辐射，作为第一辐射值；所述第一查找表，包含太阳天顶角、气溶胶光学厚度、气溶胶类型和对应的地表太阳短波辐射的数据；

步骤103、读取卫星的地表太阳短波辐射数据，为第二辐射值；与所述第一辐射值对比；若误差大于第一阈值，则该像元为校正像元；

步骤104、对所述校正像元用参数化方法计算地表太阳短波辐射，为第三辐射值。

图3为本发明包含查表法校准步骤的实施例流程图。

作为本申请进一步优化的实施例，包含步骤101、102、103、104，且在对校正像元用参数化方法计算后，还包含步骤105，实现查表法校准。具体说明如下：

步骤101、提取卫星的观测像元，在第一观测时刻的图像中，去掉有云的像元图像；

步骤102、用第一查找表计算第一观测时刻的地表太阳短波辐射，作为第一辐射值；所述第一查找表，包含太阳天顶角、气溶胶光学厚度、气溶胶类型和对应的地表太阳短波辐射的数据；

步骤103、读取卫星第一观测时刻的地表太阳短波辐射数据，为第二辐射值；与所述第一辐射值对比；若误差大于第一阈值，则该像元为校正像元；若误差小于第一阈值，则该像元为一般像元；

步骤104、对所述校正像元用参数化方法计算第一观测时刻的地表太阳短波辐射，为第三辐射值；

步骤105、提取校正像元的第三辐射值和第一辐射值进行对比得到调整值，用所述调整值对第一辐射值数据进行校准。

在步骤105中，将所有经过参数化方法计算所得的地表太阳短波辐射数据和查找表法计算出的地表太阳短波辐射数据进行对比，得出的调整值，再用调整值对那些采用查找表法计算的像元(即包括一般像元)进行校正。

例如，一个卫星图像有10000个像元，其中有100个像元既用查找表方法计算、也用参数化方法计算，这100个像元为校正像元，其他9900个像元为一般像元。将校正像元的第三辐射值与第一辐射值进行对比，得到调整值，用于校正9900个一般像元的第一辐射值。校正像元所在位置的调整值是局部参数，而一般像元的位置范围更大，当多个校正像元分布在图像区域中时，本领域技术人员能够理解，通过对校正像元的误差分析获得的校正像元调整值，可以通过插值法、范数法、拟合法等方式扩展至全局，作用于一般像元的误差。

在本申请的实施例中，对于感兴趣的重点观测区域采用参数化方法，则可以获得较为精确的反演结果，对于一般观测区域采用查找表方法可以快速得到地表入射短波太阳辐射结果。进一步扩展本实施例的运用，所述重点区域，可以是所述校正像元所在的区域，也可以是指其他具有研究价值的区域，可以用其他感兴趣的指标(例如，权重值)代替本申请实施例所述误差值。

图4为本发明包含第二观测时刻气溶胶参数校正的实施例流程图。

作为本申请优化的实施例，包含步骤101～104、106。具体说明如下，

步骤101、提取卫星的观测像元，在第一观测时刻的图像中，去掉有云的像元图像；

步骤102、用第一查找表计算第一观测时刻的地表太阳短波辐射，作为第一辐射值；所述第一查找表，包含太阳天顶角、气溶胶光学厚度、气溶胶类型和对应的地表太阳短波辐射的数据；

步骤103、读取卫星第一观测时刻的地表太阳短波辐射数据，为第二辐射值；与所述第一辐射值对比；若误差大于第一阈值，则该像元为校正像元；

步骤104、对所述校正像元用参数化方法计算第一观测时刻的地表太阳短波辐射，为第三辐射值；

步骤106、对第一观测时刻的校正像元，计算第二观测时刻的第一辐射值，如果误差仍大于第一阈值，则以第一观测时刻所述校正像元的第三辐射值为基础，根据第一查找表，得出更新的太阳天顶角、气溶胶厚度和气溶胶类型；第一观测时刻和第二观测时刻之间的间隔小于第二阈值。

回归步骤101。将第二观测时刻作为新的第一观测时刻，迭代执行以上步骤。所述第二观测时刻卫星的观测像元的图像，是指与相同位置的第一观测时刻的采集时间相邻的像元图像。

在本实施例中，例如，两个像元的观测时间间隔小于第二阈值时，在用参数化方法求出一个像元在第一观测时刻的第三辐射值后，可以简化这个像元所在像元点其第二观测时刻像元计算过程。由于卫星观测图像是按周期采集的，例如为五分钟，地表太阳短波辐射数据的变化可以忽略不计。

图5为本发明包含第二观测时刻第三辐射值替代的实施例流程图。

作为本申请优化的实施例，包含步骤101～104、107。具体说明如下，

步骤101、提取卫星的观测像元，在第一观测时刻的图像中，去掉有云的像元图像；

步骤102、用第一查找表计算第一观测时刻的地表太阳短波辐射，作为第一辐射值；所述第一查找表，包含太阳天顶角、气溶胶光学厚度、气溶胶类型和对应的地表太阳短波辐射的数据；

步骤103、读取卫星第一观测时刻的地表太阳短波辐射数据，为第二辐射值；与所述第一辐射值对比；若误差大于第一阈值，则该像元为校正像元；

步骤104、对所述校正像元用参数化方法计算第一观测时刻的地表太阳短波辐射，为第三辐射值；

步骤107、对第一观测时刻的校正像元，计算第二观测时刻的第一辐射值，如果误差仍大于第一阈值，则以第一观测时刻所述校正像元的第三辐射值作为第二观测时刻所述校正像元的地表太阳短波辐射；第一观测时刻和第二观测时刻之间的间隔小于第二阈值。

回归步骤101。将第二观测时刻作为新的第一观测时刻，迭代执行以上步骤。所述第二观测时刻卫星的观测像元的图像，是指与相同位置的第一观测时刻的采集时间相邻的像元图像。

例如，卫星观测图像上一个像元点a，a点的第一观测时刻采集的第一辐射值误差较大，用参数化方法求出的第三辐射值。在计算a点的第二观测时刻的第一辐射值时，首先判断它是否需要用参数化方法，若需要用参数化方法，则直接采用第一观测时刻的第三辐射值数据，以此类推。

需要说明的是，在本申请的各项实施例中所述步骤的组合，例如涉及步骤105，也可以用在各实施例中。如一个实施例包含步骤101、102、103、104、105、106；或另一个实施例包含步骤101、102、103、104、105、107。

再例如涉及步骤101a，用于其他实施例中，一个实施例包含步骤101、101a、102、103、104、105、106；或另一个实施例包含步骤101、101a、102、103、104、105、107。

还需要说明的是，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、商品或者设备不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、商品或者设备所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括所述要素的过程、方法、商品或者设备中还存在另外的相同要素。

以上所述仅为本申请的实施例而已，并不用于限制本申请。对于本领域技术人员来说，本申请可以有各种更改和变化。凡在本申请的精神和原理之内所作的任何修改、等同替换、改进等，均应包含在本申请的权利要求范围之内。