一种陆地云光学厚度的遥感反演方法及系统

摘要

本发明提供了一种陆地云光学厚度的遥感反演方法及系统，该方法包括：获取陆地的可见光影像，生成所述可见光影像的地表反射率；根据预设参数利用辐射传输模型建立查找表；对所述可见光影像进行辐射定标，转换为可见光影像的星上辐照度；结合多光谱影像数据，将所述可见光影像生成云掩膜图像；根据所述可见光影像的地表反射率和星上辐照度，利用所述查找表，将所述云掩膜图像中的云覆盖像元的星上辐照度反演为云光学厚度。本发明提出的陆地云光学厚度的反演方法及系统，基于云掩膜影像的提取精度，可以快速准确的实现日间云光学厚度的提取，且独立于特定的传感器和卫星，在高亮地表，如冰、雪等下垫面情形下，取得较好的反演结果。

说明书

技术领域

本发明涉及大气遥感技术领域，尤其涉及一种陆地云光学厚度的遥感反演方法及系统。

背景技术

云性质及其在时空尺度上的变化对全球气候变化研究十分重要。云光学厚度是计算云辐射效应和光学作用的重要参量之一，是地气辐射系统的主要调节参量。遥感技术为大范围实时云光学参数反演提供了有力手段。在云光学厚度大的情况下，云的反射函数可写为符合渐进理论的表达形式，其中反射辐射在非吸收的大气波段是云光学厚度、地面反射率和非对称因子的函数。当云光学厚度大于9后，无论是非守恒大气还是守恒大气，反射函数的渐进理论表达式均可以适用。Ou等人在辐射传输模式和参数化的基础上，选取了AVHRR的3.7μm和10.9μm两个通道数据，提出了反演卷云光学厚度的算法，并将该算法应用于第一次国际卫星云气候计划区域试验卷云加密观测中，结果表明，反演得到的云参数与地基和机载探测结果较为一致。Rosenfeld等人利用NOAA-AVHRR的0.65通道数据反演了云顶附近的云光学厚度信息，并与雷达回波的结果进行了对比分析。第一个全球范围的光学厚度产品由ISCCP提供，该方法是针对有效粒子半径为10μm的水云，通过查找表的方式得到，其光学厚度值域为0.5～100。

20世纪80年代以来，我国利用飞机、雷达和气象卫星遥感资料，开展了一系列云辐射特性方面的研究工作，取得了众多研究成果。利用卫星可见光通道资料反演云光学厚度的基础理论已趋于成熟，即在可见光波段，云的反射函数主要依赖于云的光学厚度，但受多种因素的共同影响，云的辐射特性复杂，在空间和时间上具有高度的可变性。现有MODIS产品的业务算法也采用近红外和可见光两个通道同时实现光学厚度和有效粒子半径的反演，其理论基础依然是平面平行理论，其反演方法会高估向阳方向的云光学厚度，低估背光方向的云光学厚度。

在实现本发明过程中，发明人发现现有技术中至少存在如下问题：

目前而言，卫星遥感云光学厚度的精度还存在着一定的问题，其中最主要的问题是云的宏观和微观不均匀性对云光学厚度的影响的程度至今尚无准确结论。对冰云而言，球形粒子假定和Mie理论已不再适用，加上冰云的时空分布变化大和其水平取向的多样性，使得冰云的反演难度更大。Nakajima指出，对于光学厚度大于4，有效粒子半径大于6μm的层状水云来说，可以用0.75μm反演光学厚度，而对于光学薄云层来说，反演结果是不确定的，解不是唯一的，由0.75反演的云光学厚度代表整层云光学厚度的20％～40％。

发明内容

(一)要解决的技术问题

本发明提出一种陆地云光学厚度的遥感反演方法及系统，可以快速准确的实现日间云光学厚度的提取，对于高亮地表，如冰、雪等下垫面情形，得到较好的反演结果。

(二)技术方案

为解决上述技术问题，本发明提供了一种陆地云光学厚度的遥感反演方法，所述方法包括：

获取陆地的可见光影像，生成所述可见光影像的地表反射率；

根据预设参数利用辐射传输模型建立查找表；

对所述可见光影像进行辐射定标，转换为可见光影像的星上辐照度；

结合多光谱影像数据，将所述可见光影像生成云掩膜图像；

根据所述可见光影像的地表反射率和星上辐照度，利用所述查找表，将所述云掩膜图像中的云覆盖像元的星上辐照度反演为云光学厚度。

优选地，所述可见光影像为极轨卫星采集的可见光波段影像数据或静止卫星采集的可见光波段影像数据。

优选地，所述生成所述可见光影像的地表反射率数据具体包括：

在预设时间段内，选取逐日同一成像时刻的可见光影像；

对所述预设时间段内选取的所述可见光影像进行大气校正；

根据校正后的可见光影像，提取地表反射率的最小值作为地表反射率。

优选地，所述预设时间段的时间步长为一个月，在所述预设时间段内地表反射率特性变化最小，且该预设时间段内至少有一天是晴空状况。

优选地，所述根据预设参数利用辐射传输模型建立查找表，具体包括：

根据预设参数，所述预设参数包括多个太阳天顶角、多个观测天顶角、相对方位角、多个地表高程、多个气溶胶光学厚度和多个云光学厚度，采用大气辐射传输模型SBDART计算多组不同地表、大气和观测几何下的云光学厚度与星上辐照度的值，生成查找表。

优选地，根据云相态的不同，所述查找表包括：水云查找表和冰云查找表。

优选地，结合多光谱影像数据，生成云掩膜图像具体包括：

通过多光谱影像数据不同波段的波段反射率或亮温与背景值的差值进行比较，预设不同的差值阈值，利用光谱阈值法生成第一层云掩膜；

通过对多光谱影像划分不同的像元矩阵，在划分出的每块像元矩阵内设定不同的动态阈值，高于所述动态阈值的像元矩阵判定为有云像元，生成第二层云掩膜；

利用云的高反射率和低亮温特性，针对可见光波段和热红外波段设定不同的覆盖阈值，高于所述覆盖阈值的像元判定为云覆盖像元，生成第三层云掩膜；

将上述三层云掩膜进行图像融合，生成云掩膜图像。

优选地，根据所述可见光影像的地表反射率和星上辐照度，利用所述查找表，将所述云掩膜图像中的云覆盖像元的星上辐照度反演为云光学厚度，具体包括：

对所述可见光影像进行地形识别，识别出下垫面为陆地的区域；

获取所述区域的地表高程和气溶胶光学厚度；

计算当前卫星的观测天顶角与观测方位角以及太阳天顶角与太阳方位角，根据观测方位角和太阳方位角计算相对方位角；

根据3.7μm通道的亮温进行云相态识别，选择相应的查找表；

根据所述可见光影像的地表反射率、星上辐照度、地表高程、观测天顶角、太阳天顶角、相对方位角以及气溶胶光学厚度，在所述查找表中选取对应的光学厚度；

利用线性插值法获取所述云掩膜图像中的云覆盖像元的光学厚度。

优选地，所述根据3.7μm通道的亮温进行云相态识别，选择相应的查找表，具体包括：

如果所述3.7μm通道的亮温高于238K，则选用水云查找表进行云参数反演；如果所述3.7μm通道的亮温低于或等于238K，则选用冰云查找表进行云参数反演。

相应的，本发明还提出了一种陆地云光学厚度的遥感反演系统，所述系统包括：

地表反射率生成模块，用于获取陆地的可见光影像，生成所述可见光影像的地表反射率；

查找表建立模块，用于根据预设参数利用辐射传输模型建立查找表；

辐射定标模块，用于对所述可见光影像进行辐射定标，转换为可见光影像的星上辐照度；

云掩膜图像生成模块，用于结合多光谱影像数据，将所述可见光影像生成云掩膜图像；

反演模块，用于根据所述可见光影像的地表反射率和星上辐照度，利用所述查找表，将所述云掩膜图像中的云覆盖像元的星上辐照度反演为云光学厚度。

(三)有益效果

本发明提供了一种陆地云光学厚度的遥感反演方法及系统，基于云掩膜影像的提取，可以快速准确的实现日间云光学厚度的提取，且独立于特定的传感器和卫星，在高亮地表，如冰、雪等下垫面情形下，取得较好的反演结果。

附图说明

通过参考附图会更加清楚的理解本发明的特征和优点，附图是示意性的而不应理解为对本发明进行任何限制，在附图中：

图1为本发明提出的一种陆地云光学厚度的遥感反演方法流程图；

图2为本发明提出的一种陆地云光学厚度的遥感反演系统模块图。

具体实施方式

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

图1为本发明提出的一种陆地云光学厚度的遥感反演方法流程图，如图1所示，该方法包括：

S101，获取陆地的可见光影像，生成所述可见光影像的地表反射率；所述可见光影像为极轨卫星采集的可见光波段影像数据或静止卫星采集的可见光波段影像数据。

其中，生成所述可见光影像的地表反射率数据具体包括：在预设时间段内，选取逐日同一成像时刻的可见光影像；对所述预设时间段内选取的所述可见光影像进行大气校正；根据校正后的可见光影像，提取地表反射率的最小值作为地表反射率，所述预设时间段的时间步长为一个月，在所述预设时间段内地表反射率特性变化最小，且该预设时间段内至少有一天是晴空状况。

本发明实施例中，结合遥感影像可见光通道以及辐射传输模型，综合考虑地表反射、大气层的反射、日-地-星间的几何关系等参数，求算出可见光波段的地表反射率。以一个月的时间窗口作为步长，将连续一年每个月高于0.005的最低地表反射率作为晴空背景值，并录入数据库。其中，高反射地表(大于0.4)的像元标识为冰雪覆盖区。

S102，根据预设参数利用辐射传输模型建立查找表，具体包括：

根据预设参数，所述预设参数包括多个太阳天顶角、多个观测天顶角、相对方位角、多个地表高程、多个气溶胶光学厚度和多个云光学厚度，采用大气辐射传输模型SBDART计算多组不同地表、大气和观测几何下的云光学厚度与星上辐照度的值，生成查找表，根据云相态的不同，所述查找表包括：水云查找表和冰云查找表。

本发明实施例中，通过设置不同日-地-星间的几何参数(包括19个太阳天顶角、19个观测天顶角和19个相对方位角)，7个地表高程、11个地表反射率、7个气溶胶光学厚度、13个云光学厚度值，采用SBDART辐射传输模型计算多组不同地表、大气和观测几何下的云光学厚度与星上辐照度的值而成N×M维的查找表。根据云相态的不同，分别针对水云和冰云构建不同的查找表。

S103，对所述可见光影像进行辐射定标，转换为可见光影像的星上辐照度；

本发明实施例中，利用辅助文件中提供的定标系数，将影像DN值转换为星上辐照度R(W/m².sr)，具体如下：

R＝DN×gain+bias

其中，gain为增益量，bias为偏移量，定标系数根据静止或极轨卫星以及仪器灵敏度的不同而会不定期的进行更新。

S104，结合多光谱影像数据，将所述可见光影像生成云掩膜图像，具体包括：

通过多光谱影像数据不同波段的波段反射率或亮温与背景值的差值进行比较，预设不同的差值阈值，利用光谱阈值法生成第一层云掩膜；通过对多光谱影像划分不同的像元矩阵，在划分出的每块像元矩阵内设定不同的动态阈值，高于所述动态阈值的像元矩阵判定为有云像元，生成第二层云掩膜；利用云的高反射率和低亮温特性，针对可见光波段和热红外波段设定不同的覆盖阈值，高于所述覆盖阈值的像元判定为云覆盖像元，生成第三层云掩膜；将上述三层云掩膜进行图像融合，生成云掩膜图像。

本发明实施例中，为保证云掩膜图像的精度，结合可见光、热红外波段3.7μm和12.0μm波段影像进行云识别。首先根据不同波段，通过该波段反射率或亮温与背景值的差值进行比较，根据不同的阈值，利用光谱阈值法生成第一层云掩膜；通过对影像划分不同的像元矩阵，针对每块子区域内设定不同的动态阈值，高于该阈值的像元判定为有云像元，生成第二层云掩膜；利用云的高反射率和低亮温特性，针对可见光波段和热红外波段设定不同的阈值，高于该阈值的像元判定为云覆盖像元，生成第三层云掩膜。三层云掩膜图像融合即为区域的云掩膜影像。

S105，根据所述可见光影像的地表反射率和星上辐照度，利用所述查找表，将所述云掩膜图像中的云覆盖像元的星上辐照度反演为云光学厚度，具体包括：

对所述可见光影像进行地形识别，识别出下垫面为陆地的区域；

获取所述区域的地表高程和气溶胶光学厚度；可见光影像对应的像元区域的地表高程和气溶胶光学厚度为现有数据。

计算当前卫星的观测天顶角与观测方位角以及太阳天顶角与太阳方位角，根据观测方位角和太阳方位角计算相对方位角；

根据3.7μm通道的亮温并进行云相态识别，选择相应的查找表；

根据所述可见光影像的地表反射率、星上辐照度、地表高程、观测天顶角、太阳天顶角、相对方位角以及气溶胶光学厚度，在所述查找表中选取对应的光学厚度；

利用线性插值法获取所述云掩膜图像中的云覆盖像元的光学厚度。

其中，根据3.7μm通道的亮温进行云相态识别，选择相应的查找表，具体包括：如果所述3.7μm通道的亮温高于238K，则选用水云查找表进行云参数反演；如果所述3.7μm通道的亮温低于或等于238K，则选用冰云查找表进行云参数反演。

本发明实施例中，采用的星上辐照度转光学厚度(τ_c)的算法如Lambertain面上的光学厚度与反射函数的表达式：

$>R({\tau }_c;\mu ,{\mu }_0,\phi )=R_{\infty }(\mu ,{\mu }_0,\phi )-\frac{m[(1-A_g{A}^{*})l-A_g{\mathrm{mn}}^2]K\left(\mu \right)K\left({\mu }_0\right)e^{-2k{\tau }_c}}{[(1-A_g{A}^{*})(1-l^2{e}^{-2k{\tau }_c})+A_g{\mathrm{mn}}^{2}l{e}^{-2k{\tau }_c}]}$>

式中，K为逃逸函数，A_g为地面反射率，g为非对称因子，μ为观测天顶角余弦，μ₀为太阳天顶角余弦，φ为相对方位角，k为散射指数(特征向量)，A*为半无限大气的球面反照率，m、n、l为常数。上述k、A*、m、n以及l等五个参数与单次散射反照率(ω₀)紧密相关，与非对称因子(g)关系较小，并且可以由相似参数s表示，具体如下：

$>s=\sqrt{\frac{1-{\omega }_0}{1-{\omega }_{0}g}}$>

云层的反射特性主要由光学厚度指数和相似参数确定，同时还与地表反射率有关。而在波长小于1μm的可见光波段，相似参数s近乎为0，因此云光学厚度可以由这一波段的反射函数获得。

本发明实施例中，根据地形识别结果，判断下垫面情形，区分陆地、内陆水体和海洋，只对标记为陆地的区域进行像元计算，提取地表高程。

本发明实施例中，以静止卫星为例，详细介绍计算卫星天顶角与方位角，具体如下：

静止卫星的天顶角Zen_v可以表示为：

Zen_v＝π/2-tan^-1[(r_s/r_e-cosγ)/sinγ]-γ

cos(γ)＝cos(L_e)cos(L_s)cos(l_s-l_e)+sin(L_e)sin(L_s)

式中，r_s为卫星高度，r_e为地球半径，L_s为静止卫星所处的纬度，L_e为像元某点所处的纬度，l_s为静止卫星所处经度，l_e为像元某点所处经度。

静止卫星的方位角可以表示为：

$>\alpha ={\tan }^{-1}\left[\frac{\tan \left|(l_s-l_e)\right|}{\sin \left(L_e\right)}\right]$>

卫星方位角定义为从正北方向起，顺时针旋转，则各个象限的方位角为：

静止卫星的天顶角和方位角是固定的，因此只需要计算一次即可。

本发明实施例中，计算太阳天顶角与方位角，具体如下：

根据太阳视轨道方程，太阳高度角Zen_s可以表示为：

式中，为地理纬度，δ为太阳赤纬，ω为太阳时角。

太阳方位角Azi_s的计算为：

$>{\mathrm{Azi}}_s=\frac{\cos \delta \sin \omega }{\cos \left({\mathrm{Zen}}_s\right)}$>

由于有云天气下水汽、臭氧等参数对可见光波段的影响远低于云的影响，在本算法中都取默认值。根据读取的太阳天顶角、太阳方位角、观测天顶角、观测方位角、地表高程和气溶胶光学厚度，结合得到的云相态以及得到的地表反射率，在查找表中选取相应的数据，通过对查找表不断降维处理的方法，采用线性插值得到整幅图像云覆盖像元的光学厚度。

图2为本发明提出的一种陆地云光学厚度的遥感反演系统模块图，如图2所示，该系统包括：

地表反射率生成模块201，用于获取陆地的可见光影像，生成所述可见光影像的地表反射率；

查找表建立模块202，用于根据预设参数利用辐射传输模型建立查找表；

辐射定标模块203，用于对所述可见光影像进行辐射定标，转换为可见光影像的星上辐照度；

云掩膜图像生成模块204，用于结合多光谱影像数据，将所述可见光影像生成云掩膜图像；

反演模块205，用于根据所述可见光影像的地表反射率和星上辐照度，利用所述查找表，将所述云掩膜图像中的云覆盖像元的星上辐照度反演为云光学厚度。

采用本发明提供了陆地云光学厚度的遥感反演方法及系统，可以快速准确的实现日间云光学厚度的提取，且独立于特定的传感器和卫星，对于高亮地表，如冰、雪等下垫面情形下的反演结果较好。其反演精度主要依赖于云掩膜影像的提取精度，可以实现太阳天顶角在75°以下像元云光学厚度的自动反演提取。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的发明保护范围应由权利要求限定。