CN112906284A - 一种浑浊水域气溶胶光学厚度与浑浊度的反演算法 - Google Patents
一种浑浊水域气溶胶光学厚度与浑浊度的反演算法 Download PDFInfo
- Publication number
- CN112906284A CN112906284A CN202110258636.6A CN202110258636A CN112906284A CN 112906284 A CN112906284 A CN 112906284A CN 202110258636 A CN202110258636 A CN 202110258636A CN 112906284 A CN112906284 A CN 112906284A
- Authority
- CN
- China
- Prior art keywords
- aerosol
- turbidity
- reflectivity
- aod
- zenith angle
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Pending
Links
- 239000000443 aerosol Substances 0.000 title claims abstract description 94
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 title claims abstract description 61
- 230000003287 optical effect Effects 0.000 title claims abstract description 34
- 230000005855 radiation Effects 0.000 claims abstract description 47
- 239000013049 sediment Substances 0.000 claims abstract description 20
- 238000000034 method Methods 0.000 claims abstract description 18
- 238000011160 research Methods 0.000 claims abstract description 10
- 238000012937 correction Methods 0.000 claims abstract description 6
- 230000002457 bidirectional effect Effects 0.000 claims abstract description 5
- 238000005315 distribution function Methods 0.000 claims abstract description 3
- 238000002310 reflectometry Methods 0.000 claims description 80
- 239000002245 particle Substances 0.000 claims description 24
- 230000005540 biological transmission Effects 0.000 claims description 19
- 238000004458 analytical method Methods 0.000 claims description 18
- 238000002834 transmittance Methods 0.000 claims description 12
- 230000008859 change Effects 0.000 claims description 10
- 239000000428 dust Substances 0.000 claims description 10
- 239000000049 pigment Substances 0.000 claims description 9
- 238000004364 calculation method Methods 0.000 claims description 8
- 239000013535 sea water Substances 0.000 claims description 8
- 239000004576 sand Substances 0.000 claims description 6
- 239000004071 soot Substances 0.000 claims description 6
- 238000010521 absorption reaction Methods 0.000 claims description 5
- 238000012417 linear regression Methods 0.000 claims description 5
- 235000002779 Morchella esculenta Nutrition 0.000 claims description 4
- 240000002769 Morchella esculenta Species 0.000 claims description 4
- 230000001419 dependent effect Effects 0.000 claims description 4
- 230000000694 effects Effects 0.000 claims description 4
- 230000003595 spectral effect Effects 0.000 claims description 3
- 230000002238 attenuated effect Effects 0.000 claims description 2
- 239000003245 coal Substances 0.000 claims description 2
- 230000000873 masking effect Effects 0.000 claims description 2
- 230000000630 rising effect Effects 0.000 claims description 2
- 239000000779 smoke Substances 0.000 claims description 2
- 238000007689 inspection Methods 0.000 claims 1
- 230000035945 sensitivity Effects 0.000 description 10
- 239000000725 suspension Substances 0.000 description 4
- 238000012360 testing method Methods 0.000 description 4
- RZVAJINKPMORJF-UHFFFAOYSA-N Acetaminophen Chemical compound CC(=O)NC1=CC=C(O)C=C1 RZVAJINKPMORJF-UHFFFAOYSA-N 0.000 description 3
- 239000003653 coastal water Substances 0.000 description 3
- 238000005259 measurement Methods 0.000 description 3
- 239000003643 water by type Substances 0.000 description 3
- 238000012952 Resampling Methods 0.000 description 2
- FAPWRFPIFSIZLT-UHFFFAOYSA-M Sodium chloride Chemical compound [Na+].[Cl-] FAPWRFPIFSIZLT-UHFFFAOYSA-M 0.000 description 2
- 239000005427 atmospheric aerosol Substances 0.000 description 2
- 239000003738 black carbon Substances 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000007613 environmental effect Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- RGNPBRKPHBKNKX-UHFFFAOYSA-N hexaflumuron Chemical compound C1=C(Cl)C(OC(F)(F)C(F)F)=C(Cl)C=C1NC(=O)NC(=O)C1=C(F)C=CC=C1F RGNPBRKPHBKNKX-UHFFFAOYSA-N 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 235000002639 sodium chloride Nutrition 0.000 description 2
- 239000011780 sodium chloride Substances 0.000 description 2
- 239000000126 substance Substances 0.000 description 2
- OKTJSMMVPCPJKN-UHFFFAOYSA-N Carbon Chemical compound [C] OKTJSMMVPCPJKN-UHFFFAOYSA-N 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000003915 air pollution Methods 0.000 description 1
- 230000008033 biological extinction Effects 0.000 description 1
- 229910052799 carbon Inorganic materials 0.000 description 1
- 239000011362 coarse particle Substances 0.000 description 1
- 230000000052 comparative effect Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 238000001514 detection method Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000010586 diagram Methods 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 235000010755 mineral Nutrition 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 239000000203 mixture Substances 0.000 description 1
- 150000002823 nitrates Chemical class 0.000 description 1
- 239000013618 particulate matter Substances 0.000 description 1
- 230000010287 polarization Effects 0.000 description 1
- 230000002265 prevention Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 238000000611 regression analysis Methods 0.000 description 1
- 230000004044 response Effects 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 230000007480 spreading Effects 0.000 description 1
- 150000003467 sulfuric acid derivatives Chemical class 0.000 description 1
- 238000012546 transfer Methods 0.000 description 1
- 238000012795 verification Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/25—Design optimisation, verification or simulation using particle-based methods
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2113/00—Details relating to the application field
- G06F2113/08—Fluids
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- General Engineering & Computer Science (AREA)
- General Physics & Mathematics (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
Abstract
本发明属于卫星遥感技术领域，具体为一种浑浊水域气溶胶光学厚度与浑浊度的反演算法。本发明方法利用在浑浊水域近红外波段(0.86μm)的离水辐射(WLR)不再近似为0，而短波近红外波段(2.1μm)处的WLR可近似为0这一特点，选取2.1μm通道来反演气溶胶光学厚度，通过海洋双向反射分布函数(BRDF)模式来定义浑浊度，从而代替传统的近红外波段大气校正的方法来研究泥沙量，使用0.86μm通道来反演出浑浊度，从而建立浑浊度与泥沙量之间的关系。本发明算法可以填补以及改进沿海区域卫星遥感的AOD产品，同时可以反映出沿海区域的泥沙量浓度，具有重要的实际应用意义。
Description
技术领域
本发明属于卫星遥感技术领域，具体涉及浑浊水域气溶胶光学厚度与浑浊度的反演算法。
背景技术
气溶胶光学厚度(AOD)是气溶胶最重要的参数之一，代表着介质的消光系数在垂直方向上的积分，它可以表征大气的环境质量和混浊程度。当前AOD的监测主要有地基观测和卫星遥感两种手段。地基观测方法主要使用太阳光度计进行测量。然而地面监测往往只能得到站点周围的信息，没有大范围的整体观测数据，并且地面站点的数量及分布也受到限制。但是气溶胶卫星遥感反演就不存在这些问题，并且空间分辨率高、探测范围广、精度更高、获取数据快(Wei et al.,2019)。近年来，随着遥感卫星技术的发展，对于大气气溶胶的反演研究不断深入；定量反演气溶胶光学厚度也取得了很大的进展，同时气溶胶的精确反演也对于研究气候变化与空气污染防治做出了重大的贡献。
现有气溶胶卫星反演算法都是针对不同地表类型和气溶胶类型组成的不同，以不同的原理来对气溶胶进行反演。目前，AOD反演算法主要有：暗像元法、深蓝算法、SARA算法、结构函数法、多角度遥感法、双星协同法、偏振特性遥感等。其中暗目标算法又称为浓密植被算法，其核心原理是植被覆盖地区的红光波段的地表反射率和蓝光波段的地表反射率与短波近红外波段的地表反射率存在一定的线性关系(Remer et al.,2008)。深蓝算法主要选取蓝光波段作为工作波段(Hsu et al.,2006)。深蓝算法是利用在蓝波段大气反射较强，地表反射较弱的特点，并假设同期的地表反射率不变，将晴天的地表反射率代入公式反演气溶胶。SARA算法在设定地表为朗伯体、近似单次散射和气溶胶单次散射率(SingleScattering Albedo，SSA)及不对称因子(Asymmetry Factor，AF)不存在空间变化前提条件下，引入大气后向散射率、气溶胶散射相函数的近似表达函数(Bilal,Nichol,&Chan,2014)。结构函数方法反演AOD的实质就是先找到该地区的地表反射率的分布规律。假设分布规律保持不变，则表观反射率的空间变化认为是大气气溶胶光学厚度产生的贡献，从而反演出气溶胶光学厚度(Tanre et al.,1988)。
为了方便对水色遥感的研究，根据水体中悬浮和溶解物对水体光学属性影响不同，Morel和Prieur(1977年)把海洋水体分为一类和二类水体，一类水体主要指光学属性变化基本来自于浮游植物影响的水体，通常呈深蓝色，例如大洋。而二类水体不仅受浮游植物和相关的颗粒影响，还受到悬浮颗粒和沉积物的影响，对于一些近海浑浊水域，目前利用卫星遥感技术反演AOD卫星是一个很大的挑战。浑浊水域主要指近岸、河口等受陆源物质排放影响较为严重的地方。作为最主要的气溶胶卫星遥感数据源之一的MODIS产品反演陆地与海洋Ⅰ类水体AOD算法较多而且比较成熟，但是现有的算法都无法反演Ⅱ类水体的气溶胶参数。
在一些海浑浊水域反演AOD十分困难，是因为气溶胶浓度关系，且气溶胶的组成成分也十分复杂且多样，包括沙漠尘埃、海盐物质等，这为AOD的精准遥感反演增加了难度，并且一些沿海由于大量泥沙的堆积导致水体比较复杂。MODIS反演海洋上的AOD目前采用的是MODIS DT-ocean算法，DT-ocean算法(Levy et al.,2013)使用了不同精细气溶胶粒子和粗气溶胶粒子的加权组合在6个波段(0.55、0.65、0.86、1.24、1.63和2.11μm)的光谱反射率，通过建立查找表反演出AOD与FMF精细模式分数。此外，查找表假设除0.55微米(使用0.005的固定离水反射率)之外，所有波段的water-leaving radiance为零。而传统大气校正中采用的近红外波段(0.86μm)water-leaving radiance近似为零的假设却只适用于公海，对于混浊的沿海水域(二类水体)就不再适用。因为在浑浊水域，浅水海底的反射(尤其是在0.66μm波段)和水中悬浮或溶解的颗粒物质的反射(尤其是在0.55、0.86μm波段)会贡献water-leaving radiance。因此，在DT-ocean算法中浑浊水域像素将会被标记，并且不被反演。
目前，一些沿海的AOD卫星产品的精度和适用性还很低。NASA的MODIS气溶胶产品在Ⅱ类水体上空也通常为空白值,或者与实际测量值存在较大的误差而无法使用。Wang(2017年)分析了在全球和区域范围内无云条件下的MODIS AOD反演不可用性，以揭示AOD不能被反演仅仅是因为水的混浊度(而不是其他因素，如云层)，通过分析发现在所有公海上，这种数据可用性几乎是完全的，并且朝着海岸线急剧下降(下降90-100％)，在全球平均水平上，AOD在沿海水域的不可用性约为20％。气溶胶是研究这些区域大气环境质量的重要因子之一和遥感定标的重要参数，并且约60％的人口生活在沿海地区，扩大这些地区上卫星遥感数据集至关重要，因此对Ⅱ类水体开展气溶胶性质反演及分布的研究具有重要意义。
沿海水域的特征通常是由于海床再次悬浮或河流排放的颗粒而产生的高浓度悬浮有机物和无机物，悬浮物浓度直接影响水体的浑浊度和颜色。长江每年平均向河口排放约9×1011m3的水和4×108吨的泥沙,传送的泥沙量主要受枯水期和汛期的影响，汛期(6月至9月)时输送的泥沙量约占全年输送泥沙总量的87％(chenet al.2003年)。夏季，长江汇入东海扩散之后形成了浑浊海域，向东延伸了104-105km2的区域(Zhang et al.2007)。长江约50％的泥沙堆积在河口处(Shen et al.2001)，SSC值跨越两个数量级，范围从20mg/l到2500mg/l甚至更多。
对于沿海悬浮物浓度的卫星遥感分析目前也已经有了一些进展。如Miller和McKee(2004)分析了大气校正后的MODIS第1波段反射率与总悬浮物(TSM)浓度之间的回归关系，发现MODIS Terra 250m波段1(620～670nm)数据与实地测量的TSM数据之间可以建立起良好的线性关系。Doxaranet al.(2002,2009)一种从卫星数据的可见光和近红外(NIR)波长来确定悬浮物浓度的方法，通过这些大量的现场测量，在可见光到近红外波段将遥感反射率和悬浮物浓度之间建立了经验关系，发现近红外波段XS3(790–890nm)的遥感反射率Rrs和可见光波段XS1(500–590nm)与XS2(790–890nm)遥感反射率的比(Rrs(XS3)/Rrs(XS1)和Rrs(XS3)/Rrs(XS2))与SPM之间有很好的相关性。这些算法可以适用于高悬浮物浓度的区域，但是用于低悬浮物浓度区域时会产生较大的误差。
针对以上问题，本发明提出了一个改进的浑浊水域气溶胶光学厚度反演算法，主要利用在浑浊水域近红外波段(0.86μm)的water-leaving radiance不再近似为0，而短波近红外波段(2.1μm)处的water-leaving radiance可近似为0这一特点，选取了2.1μm来反演气溶胶光学厚度，并且不采用传统的近红外波段大气校正的方法来研究泥沙量，而是通过海洋BRDF模式来定义浑浊度，并且反演出0.86μm的浑浊度，从而建立浑浊度与泥沙量之间的关系。
发明内容
本发明的目的在于提供一种能够正确反映大气的环境质量的浑浊水域气溶胶光学厚度与浑浊度的反演算法。
本发明提供的浑浊水域气溶胶光学厚度与浑浊度的反演算法，内容主要包括四部分：
一、通过四种气溶胶粒子(沙尘型，水溶性，海洋性，煤烟型)的不断迭代来确定研究区域的气溶胶类型，利用MODIS数据的2.1μm波段，通过6S辐射传输模型[1]生成查找表，来反演出AOD；
二、使用海洋BRDF模型[2]来模拟表面反射的情况，并且将海洋BRDF模型做出改进，定义一个新的参数：浑浊度；
三、再次利用6S辐射传输模式来反演出0.86μm处的浑浊度；
四、结合站点的实测数据验证该算法反演AOD的适用性，并进一步结合泥沙量数据，与反演得到的地表浑浊度进行回归分析，研究他们之间的响应关系。
具体步骤为：
(1)确定研究区域的气溶胶类型；将气溶胶类型看做沙尘型、水溶性、海洋性和煤烟型四种气溶胶粒子的组合；假设四种气溶胶粒子的体积浓度百分比范围后，让四种粒子在各自区间进行迭代，利用6S辐射传输模型计算出理论表观反射率，将研究区域每个格点的理论表观反射率与观测表观反射率最匹配时所对应的四种粒子百分比的组合记录下来，所有格点记录的四种粒子百分比的平均值就代表该地区的气溶胶类型；
(2)选择2.1μm通道进行反演；因为2.1μm处的离水辐射率可近似为0；利用6S辐射传输模型，输入相关参数计算出表观反射率理论值，得到AOD、地表反射率、太阳天顶角、卫星天顶角、相对方位角变化的查找表；根据每个格点上的地表反射率(由MOD09GA得到，这里，MOD09GA提供500m分辨率的MODIS波段1-7的每日表面反射率，并提供1km的观测和地理位置统计数据)、太阳天顶角、卫星天顶角、相对方位角(由MOD02得到，这里MOD02提供1km分辨率的MODIS各波段的表观反射率、地理坐标信息以及几何参数)，对查找表进行插值，找到表观反射率观测值(由MOD02得到)与表观反射率理论值最匹配时所对应的AOD；
(3)使用海洋BRDF模式，定义浑浊度；双向反射率表示遥感器观测到的地球表面反射率不仅取决于目标本身，还取决于太阳和传感器相对于目标的位置；双向反射分布函数(BRDF)模型是双向反射率的数学表达式，用于计算在不同的太阳、传感器几何参数下观察到的反射率；将海洋BRDF模型做出改进，定义一个新的参数：浑浊度；
(4)使用0.86μm反演出浑浊度；依旧利用6S辐射传输模型生成浑浊度、AOD、太阳天顶角、卫星天顶角、相对方位角变化的查找表；根据每个格点上的AOD(由前一步的反演得到)、太阳天顶角、卫星天顶角、相对方位角(由MOD02得到)，对查找表进行插值，找到表观反射率观测值(由MOD02得到)与表观反射率理论值最匹配时所对应的浑浊度；
(5)AOD反演结果的精度检验：对地面站点AOD资料进行时空匹配；将气溶胶反演结果与地基数据观测的AOD进行线性回归分析，分析反演结果与实测值的相关性；
(6)建立浑浊度与泥沙量之间的关系：对泥沙量资料进行时空匹配；通过经验模型分析法、理论模型分析法、半分析模型分析法等多种方法，找到最相匹配的浑浊度与泥沙量之间的经验公式；将经验公式用于其他个例以检验反演浑浊度的精度。
本发明的反演算法中，涉及有关算法原理有：
1、AOD反演原理
当通过卫星传感器进行测量时，大气和地表的贡献取决于大气和地面光学参数以及相对方位角、太阳天顶角以及卫星天顶角；卫星传感器所接收到的地面目标的总辐射亮度，是经大气衰减后的地表面辐射亮度和大气本身的路径辐射之和(Tanré et al.,1996)；大气程辐射是指太阳辐射在大气传输过程中各组分及气溶胶微粒散射后直接到达传感器的辐射。
对于无云且大气水平均一的地-气系统，朗伯地表(各向同性反射)上空卫星观测的TOA方程可以表达为(Kaufman et al.,1988)：
式中，θv，θs，
分别为观测天顶角、太阳天顶角、相对方位角，F(θs)是向下总透过率，T(θv)是向上总透过率，S代表大气半球反照率；ρa为大气反射率，又称路径辐射，它是由气溶胶和大气其他组分散射共同造成的；对于单次散射，并且进一步假设光学厚度很小，此时路径辐射ρa与气溶胶光学厚度τa、气溶胶散射相函数Pa以及单次散射反照率ω0有关，可以用公式表示：
其中，ρm是分子散射造成的路径辐射，μv是卫星天顶角余弦，μs是太阳天顶角余弦；公式(1)中F(θs)、T(θv)和S取决于τa、Pa以及ω0。将(2)代入(1)中，整理得到：
根据方程(3)可知，表观反射率不仅是气溶胶光学厚度的函数，也是地表反射率的函数。其中，角度数据与表观反射率可由卫星数据得到，而单次散射反照率、散射相函数等参数可通过符合实际情况的大气模式及气溶胶模型确定。实际操作中，我们通过辐射传输模型输入相应参数来获得这些量。因此只要求得地表反射率和气溶胶模型，理论上便可以通过公式(3)反演气溶胶光学厚度。
2、海洋BRDF模型
通过考虑白帽，太阳耀斑和浑浊度的影响来计算海洋表面的BRDF；假定海洋表面的反射率ρos(λ)是取决于三个分量的总和(Koepke，1984)：
ρos(θv,θs,Φ,λ)＝ρwc(λ)+{1-W}·ρgl(θv,θs,Φ,λ)+{1-ρwc(λ)}·ρsw(θv,θs,Φ,λ) (4)。
其中，ρwc(λ)是海洋白帽反射，ρgl(λ)是海洋表面的镜面反射率，ρsw(λ)是海水产生的散射反射率，W是覆盖有白帽的的相对面积，ρsw(θs,θv,φ,λ)是刚好在海面(level0+)上观测到的反射率，该反射率与Rw(Rw表面正下方的上升流光谱辐照度Eu(λ)与下降流辐照度Ed(λ)之比)有关；如果我们假设海洋是朗伯反射，则ρsw(θs,θv,φ,λ)可以表示为:
而Rw取决于海水的浑浊度，本发明将下面的色素浓度C定义为浑浊度。其中td是下行辐射的透射比，tu是上升辐射的透射比，都可以使用空气-水界面的菲涅耳反射系数Ra-w(θs,θv,φ)计算[3]。
辐照度反射率Rw(λ)特别取决于海水的固有光学特性：总吸收系数a(λ)[m-1]和总后向散射系数bb(λ)[m-1]。例如，Morel&Prieur(1977)表明，在一个很好的近似值中(当a(λ)<<1时)可以表示为：
Morel将总后向散射系数分为2个部分：
其中bw(λ)是水的分子散射系数，
是色素的反向散射比，与色素浓度C(Chl a+Pheo a)[mg·m-3]和波长[μm]有关：
其中b是色素的散射系数，计算公式为：
b＝0.3C0.62 (9)
总吸收系数写为：
a(λ)＝u(λ)·Kd(λ) (10)
其中u(λ)是与波长有关的函数，计算如下：
Kd(λ)是向下辐射的总漫射衰减系数，由下式给出：
Kd(λ)＝Kw(λ)+χc(λ)Ce(λ) (12)
根据Morel的模型，反射率Rw(λ)的计算仅取决于浑浊度C，并且越浑浊的水域，C值相应也就越大。
3、浑浊度的反演以及与泥沙量建立关系
使用0.86μm反演出浑浊度，因为在干净的水域，0.86μm的water-leavingradiance可近似为0，然而在浑浊水域2.1μm波段的water-leaving radiance可近似为0，但是0.86μm的water-leaving radiance却不能忽略，这其中就包含了水域浑浊度的信息。依旧利用6s辐射传输模型生成浑浊度、AOD、太阳天顶角、卫星天顶角、相对方位角变化的查找表；根据每个格点上的AOD(由前一步的反演得到)、太阳天顶角、卫星天顶角、相对方位角(由MOD02得到)，对查找表进行插值，找到表观反射率观测值(由MOD02得到)与表观反射率理论值最匹配时所对应的浑浊度。接下来对泥沙量资料进行时空匹配，通过经验模型分析法、理论模型分析法、半分析模型分析法等多种方法找到最相匹配的浑浊度与泥沙量之间的经验公式，最后将经验公式用于其他个例以检验反演浑浊度的精度。
附图说明
图1为2020年8月15日MYD04_L2 AOD分布图。
图2为算法流程图。
图3为绿色区域为本文的研究区域，红点为地面观测站点。
图4为SONET观测站2020年8月16日上海站点不同气溶胶组分百分比图。其中，BC：黑碳气溶胶；BrC：有机碳气溶胶；CM：粗粒子气溶胶(矿物沙尘或海盐)；FS：精细模式散射成分(如硫酸盐、硝酸盐和轻质非吸收有机物)；AW：水溶型气溶胶。
图5自定义气溶胶类型下表观反射率理论值与表观反射率观测值的线性回归。
图6为2020年8月15日上海沿海区域的2.1μm通道反射率分布图。其中，(a)为表观反射率分布图，(b)为地表反射率分布图。
图7为2.1μm表观反射率(a)与地表反射率(b)的统计图。
图8为表观反射率对于AOD的敏感性实验。其中，(a)是0.65μm通道，(b)是0.86μm通道。
图9为表观反射率对于AOD的敏感性实验(2.1μm通道)。
图10为表观反射率对于地表反射率的敏感性实验。其中，(a)是2.1μm通道，(b)是0.65μm通道。
图11为MYD04_L2 AOD产品以及反演得到的AOD(红色边框区域为反演区域)。
具体实施方式
下面结合附图和实施例对本算法进一步说明。
实施例1
国内外研究者通过不同的卫星数据来实现对AOD的反演已经取得了一定的成果，但是目前的算法还存在一定的适用性的问题，例如现在常见的海岸气溶胶反演算法只能适用于一类水体，而对于一些二类水体而言，这些算法还需要进行改进。在图1中可以看到MODIS的AOD产品在上海沿海区域就是缺测值。具体技术路线如图2所示。
研究区域和时间的选择：
地点选择上海沿海区域(经度120.9-122.5，纬度30.4-32.3)，也就是长江口到杭州湾这一区域时间选择2020年7,8,9月的晴天。
辐射传输模式的选择：
选择了6S大气辐射传输模型，6S大气辐射传输模型是模拟太阳辐射在地气系统中的传输过程的模型。主要作用是模拟传感器在无云条件下，考虑了非朗伯体情况下，真实大气各组分对于太阳辐射的吸收和散射作用后，卫星传感器理论上应该接受到的辐亮度。改进后的6S模型采用了最新的散射计算方法，考虑了多种散射，包括：分子散射、多次散射、气溶胶散射，使得计算得到的太阳光谱散射值精度更高(Lee et al.,2004)。
卫星数据的选择与处理：
卫星数据选择了MOD021km与MOD09GA；辐射定标：将原始的数字信号转换为反射率和亮温等实际物理参量；几何校正：将数据准确定位到特定的地理坐标系；数据重采样：将不同分辨率下的MODIS波段重采样到相同分辨率；裁剪：将数据裁剪成我们需要的范围。
气溶胶模式的选择：
气溶胶类型可以看做沙尘型，水溶性，海洋性，煤烟型气溶胶粒子的组合。根据Fan等人在2015年对选取的几个中国地区aeronet站点进行了对气溶胶各种组分的百分比的统计。发现每个地区水溶性所占比分最高，均在40％以上，而沙尘性均在40％以下，海洋性和煤烟性最少，均低于10％。同时根据SONET上海站点对不同气溶胶粒子的观测可以发现，水吸收性气溶胶粒子的占比是最高的，达到了80％，而黑碳型气溶胶粒子时占比最少的。
因此可以随机假设四种气溶胶粒子的体积浓度百分比为：沙尘型x1，水溶性x2，海洋性x3，煤烟型x4。对于各组分的取值可以限定为0≤x1≤20，40≤x2≤80，0≤x4≤5，x3＝100-(x1+x2+x4)。将x1,x2,x3,x4在各自区间进行迭代，迭代步长均为1，计算MOD02得到的四个角度的平均值作为输入的角度参数。使用站点资料的AOD作为输入。每迭代一次就运行6S得到F(θs)·T(θv)总透射率、S大气半球反照率、ρa路径辐射。现在共有151×203个格点，对于每个格点，结合地表反射率、F(θs)·T(θv)总透射率、S大气半球反照率、ρa路径辐射计算出理论表观反射率ρ*，然后与观测表观反射率值ρ进行对比。计算出ε＝|ρ*-ρ|，找到ε最小值的情况，并将对应的x1,x2,x3,x4记录下来。所有格点记录的x1,x2,x3,x4的平均值就代表该地区的气溶胶类型。经过计算得到了表1中的自定义气溶胶模型：
表1自定义的气溶胶类型
| 沙尘型 | 水溶型 | 海洋型 | 煤烟型 |
| 19 | 61 | 15 | 5 |
为了验证自定义的气溶胶模型的精确性，于是将自定义气溶胶模型作为6S模式的输出参数，计算得到这个区域的表观反射率理论值，再和MODIS观测的表观反射率进行线性回归分析，如图5所示，发现相关系数达到了0.8698，表观反射率理论值和观测值拟合得较好，说明自定义的气溶胶模型很好地反映了这个地区的气溶胶类型。
敏感性实验：
在反演之前需要做一个敏感性实验，验证表观反射率是否对气溶胶光学厚度以及地表反射率敏感。
图6显示的是2.1μm表观反射率以及地表反射率分布图，可以发现在长江口以及上海几个港口的下游都呈现表观反射率以及地表反射率数值的大值区，但是总体来说，表观反射率以及地表反射率的数值都在0.15以下。
图7是2.1μm表观反射率和地表反射率的统计分布图，可以看到表观反射率主要集中在0.06到0.08之间，而地表反射率主要集中在0.065到0.08之间，比常用的0.65μm的地表反射率值小得多。于是利用6S辐射传输模型，模拟在0.65μm通道和2.1μm通道下，表观反射率对气溶胶光学厚度以及地表反射率变化的变化。
如图8所示，对于0.65μm波段，在比较清澈的水域时(地表反射率为0.05)，表观反射率随AOD是明显的上升趋势，敏感性较好，但是当水域比较浑浊时(地表反射率为0.2)，表观反射率先是随AOD下降后趋于平滑，敏感性不再显著。对于0.86μm波段，无论是在清澈的水域还是浑浊水域，敏感性都不够显著。如图9所示，对于2.1μm波段，无论是在清澈的水域还是浑浊的水域，表观反射率都随AOD有明显的下降的趋势，呈现出较好的敏感性。因此在浑浊水域，使用2.1μm波段来反演比使用0.65μm波段更为准确。从图10可以看出无论是在0.65μm波段还是在2.1μm波段表观反射率都对地表反射率表现出较好的敏感性，且都呈现出表观反射率随地表反射率增大而增大的趋势。
建立AOD查找表：
本算法选择2.1微米通道进行反演，因为浑浊水域2.1μm处的water-leavingradiance可近似为0。确定好气溶胶模型之后，利用6S辐射传输模型输入相关参数计算出表观反射率理论值，得到AOD、地表反射率、太阳天顶角、卫星天顶角、相对方位角变化的查找表。其中建立的查找表中各变量的范围为：AOD:0.1-0.5(步长为0.01)，地表反射率：0-0.2(步长为0.02)，太阳天顶角：22-25(步长为1)，卫星天顶角：2-20(步长为3)，相对方位角：211-226(步长为3)。
对MODIS表观反射率数据进行气体校正、云掩膜后，根据每个格点上的地表反射率(由MOD09GA得到)、太阳天顶角、卫星天顶角、相对方位角(由MOD02得到)，对查找表进行插值，找到表观反射率观测值(由MOD02得到)与表观反射率理论值最匹配时所对应的AOD。
图11是反演得到的AOD(红色边框区域)以及MODIS AOD产品分布图，趋势大致吻合，根据反演的AOD分布，可以发现在长江口以及上海几个港口的下游地区的AOD都呈现大值，达到了0.45以上，并且这两个地方的地表反射率也是大值区(水域较浑浊)，说明这两个地方受到了沿海工业以及船舶的污染影响。
AOD的对比验证：
对地面站点AOD资料进行时空匹配。将气溶胶反演结果与地基数据观测的AOD进行线性回归分析，分析反演结果与实测值的相关性。
使用BRDF模式定义浑浊度：
通过考虑白帽，太阳耀斑和浑浊度的影响来计算海洋表面的BRDF。假定海洋表面的反射率ρos(λ)是取决于三个分量的总和(Koepke，1984)：
ρos(θv,θs,Φ,λ)＝ρwc(λ)+{1-W}·ρgl(θv,θs,Φ,λ)+{1-ρwc(λ)}·ρsw(θv,θs,Φ,λ)
如果我们假设海洋是朗伯反射，则ρsw(θs,θv,φ,λ)可以表示为:
如公式(4)-(12)所示，Rw取决于海水的浑浊度，本发明将该色素浓度C定义为浑浊度。
建立浑浊度查找表：
使用0.86μm反演出浑浊度，因为在干净的水域，0.86μm的water-leavingradiance可近似为0，然而在浑浊水域2.1μm波段的water-leaving radiance可近似为0，但是0.86μm的water-leaving radiance却不能忽略，这其中就包含了水域浑浊度的信息。依旧利用6s辐射传输模型生成浑浊度、AOD、太阳天顶角、卫星天顶角、相对方位角变化的查找表；根据每个格点上的AOD(由前一步的反演得到)、太阳天顶角、卫星天顶角、相对方位角(由MOD02得到)，对查找表进行插值，找到表观反射率观测值(由MOD02得到)与表观反射率理论值最匹配时所对应的浑浊度。
同泥沙量数据建立经验关系：
对泥沙量资料进行时空匹配，通过经验模型分析法、理论模型分析法、半分析模型分析法等多种方法找到最相匹配的浑浊度与泥沙量之间的经验公式，最后将经验公式用于其他个例以检验反演浑浊度的精度。
Claims (5)
1.一种浑浊水域气溶胶光学厚度与浑浊度的反演算法，其特征在于，具体步骤为：
(1)确定研究区域的气溶胶类型；
将气溶胶类型看做沙尘型、水溶性、海洋性和煤烟型四种气溶胶粒子的组合；假设四种气溶胶粒子的体积浓度百分比范围后，让四种粒子在各自区间进行迭代，利用6S辐射传输模型计算出理论表观反射率，将研究区域每个格点的理论表观反射率与观测表观反射率最匹配时所对应的四种粒子百分比的组合记录下来，所有格点记录的四种粒子百分比的平均值就代表该地区的气溶胶类型；
(2)选择2.1μm通道进行反演AOD；
由于2.1μm处的water-leaving radiance近似为0；利用6S辐射传输模型，输入相关参数计算出表观反射率理论值，得到AOD、地表反射率、太阳天顶角、卫星天顶角、相对方位角变化的查找表；根据每个格点上的地表反射率、太阳天顶角、卫星天顶角、相对方位角，对查找表进行插值，找到表观反射率观测值与表观反射率理论值最匹配时所对应的AOD；
(3)使用海洋BRDF模式，定义浑浊度；
双向反射率表示遥感器观测到的地球表面反射率不仅取决于目标本身，还取决于太阳和传感器相对于目标的位置；双向反射分布函数(BRDF)模型是双向反射率的数学表达式，用于计算在不同的太阳、传感器几何参数下观察到的反射率；根据海洋BRDF模型，定义一个新的参数：浑浊度；
(4)使用0.86μm反演出浑浊度；
依旧利用6S辐射传输模型生成浑浊度、AOD、太阳天顶角、卫星天顶角、相对方位角变化的查找表；根据每个格点上的AOD、太阳天顶角、卫星天顶角、相对方位角，对查找表进行插值，找到表观反射率观测值与表观反射率理论值最匹配时所对应的浑浊度；
(5)AOD反演结果的精度检验；
对地面站点AOD资料进行时空匹配；将气溶胶反演结果与地基数据观测的AOD进行线性回归分析，分析反演结果与实测值的相关性；
(6)建立浑浊度与泥沙量之间的关系；
对泥沙量资料进行时空匹配；通过经验模型分析法、理论模型分析法、半分析模型分析法等多种方法，找到最相匹配的浑浊度与泥沙量之间的经验公式；将经验公式用于其他个例以检验反演浑浊度的精度。
2.根据权利要求1所述的浑浊水域气溶胶光学厚度与浑浊度的反演算法，其特征在于，步骤(1)的具体流程为：随机假设四种气溶胶粒子的体积浓度百分比为：沙尘型x1，水溶性x2，海洋性x3，煤烟型x4；各组分的取值限定为0≤x1≤20，40≤x2≤80，0≤x4≤5，x3＝100-(x1+x2+x4)；将x1,x2,x3,x4在各自区间进行迭代，迭代步长均为1，计算MOD02得到的四个角度的平均值作为输入的角度参数；使用站点资料的AOD作为输入；每迭代一次就运行6S得到F(θs)·T(θv)总透射率、S大气半球反照率、ρa路径辐射；对于每个格点，结合地表反射率、F(θs)·T(θv)总透射率、S大气半球反照率、ρa路径辐射计算出理论表观反射率ρ*，然后与观测表观反射率值ρ进行对比；计算出ε＝|ρ*-ρ|，找到ε最小值的情况，并将对应的x1,x2,x3,x4记录下来；所有格点记录的x1,x2,x3,x4的平均值就代表该地区的气溶胶类型。
3.根据权利要求1所述的浑浊水域气溶胶光学厚度与浑浊度的反演算法，其特征在于，步骤(2)中所述利用6S辐射传输模型输入相关参数计算出表观反射率理论值，得到AOD、地表反射率、太阳天顶角、卫星天顶角、相对方位角变化的查找表；查找表中各变量的范围为：AOD:0.1-0.5，步长为0.01；地表反射率：0-0.2，步长为0.02；太阳天顶角：22-25，步长为1；卫星天顶角：2-20，步长为3；相对方位角：211-226，步长为3；
对MODIS表观反射率数据进行气体校正、云掩膜后，根据每个格点上的地表反射率、太阳天顶角、卫星天顶角、相对方位角，对查找表进行插值，找到表观反射率观测值与表观反射率理论值最匹配时所对应的AOD。
4.根据权利要求1所述的浑浊水域气溶胶光学厚度与浑浊度的反演算法，其特征在于，所述AOD反演，具体流程如下：
当通过卫星传感器进行测量时，大气和地表的贡献取决于大气和地面光学参数以及相对方位角、太阳天顶角以及卫星天顶角；卫星传感器所接收到的地面目标的总辐射亮度，是经大气衰减后的地表面辐射亮度和大气本身的路径辐射之和；大气程辐射是指太阳辐射在大气传输过程中各组分及气溶胶微粒散射后直接到达传感器的辐射；
对于无云且大气水平均一的地-气系统，朗伯地表上空卫星观测的TOA方程为：
式中，θv，θs，
分别为观测天顶角、太阳天顶角、相对方位角，F(θs)是向下总透过率，T(θv)是向上总透过率，S代表大气半球反照率；ρa为大气反射率，又称路径辐射，它是由气溶胶和大气其他组分散射共同造成的；对于单次散射，并且进一步假设光学厚度很小，此时路径辐射ρa与气溶胶光学厚度τa、气溶胶散射相函数Pa以及单次散射反照率ω0有关，公式表示为：
其中，ρm是分子散射造成的路径辐射，μv是卫星天顶角余弦，μs是太阳天顶角余弦；公式(1)中F(θs)、T(θv)和S取决于τa、Pa以及ω0；将(2)式代入(1)式中，整理得到：
根据方程(3)可知，表观反射率不仅是气溶胶光学厚度的函数，也是地表反射率的函数；其中，角度数据与表观反射率由卫星数据得到，而单次散射反照率、散射相函数等参数通过符合实际情况的大气模式及气溶胶模型确定；通过公式(3)反演气溶胶光学厚度。
5.根据权利要求1所述的浑浊水域气溶胶光学厚度与浑浊度的反演算法，其特征在于，步骤(3)中所述使用海洋BRDF模式定义浑浊度，具体内容如下：
考虑白帽，太阳耀斑和浑浊度的影响来计算海洋表面的BRDF，假定海洋表面的反射率ρos(λ)是取决于三个分量的总和：
ρos(θv,θs,Φ,λ)＝ρwc(λ)+{1-W}·ρgl(θv,θs,Φ,λ)+{1-ρwc(λ)}·ρsw(θv,θs,Φ,λ) (4)
其中，ρwc(λ)是海洋白帽反射，ρgl(λ)是海洋表面的镜面反射率，ρsw(λ)是海水产生的散射反射率，W是覆盖有白帽的相对面积，ρsw(θs,θv,φ,λ)是刚好在海面(level0+)上观测到的反射率，该反射率与辐照度反射率Rw有关，Rw是表面正下方的上升流光谱辐照度Eu(λ)与下降流辐照度Ed(λ)之比；假设海洋是朗伯反射，则ρsw(θs,θv,φ,λ)表示为:
而Rw取决于海水的浑浊度，即色素浓度C，定义C为浑浊度；其中，td是下行辐射的透射比，tu是上升辐射的透射比，都可以使用空气-水界面的菲涅耳反射系数Ra-w(θs,θv,φ)计算；
辐照度反射率Rw(λ)特别取决于海水的固有光学特性：总吸收系数a(λ)和总后向散射系数bb(λ)；在一个近似值中，当a(λ)<<1时，Rw(λ)表示为：
其中，a定义如下：
为了最小化计算，取为常数a＝0.485；
Morel将总后向散射系数bb(λ)分为2个部分：
其中,bw(λ)是水的分子散射系数，
是色素的反向散射比，与色素浓度C和波长有关：
其中,b是色素的散射系数，计算公式为：
b＝0.3C0.62, (10)
总吸收系数写为：
a(λ)＝u(λ)·Kd(λ), (11)
其中,u(λ)是与波长有关的函数，计算如下：
Kd(λ)是向下辐射的总漫射衰减系数，由下式给出：
Kd(λ)＝Kw(λ)+χc(λ)Ce(λ) (13)
根据Morel的模型，反射率Rw(λ)的计算仅取决于浑浊度C，并且越浑浊的水域，C值相应也就越大。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110258636.6A CN112906284A (zh) | 2021-03-09 | 2021-03-09 | 一种浑浊水域气溶胶光学厚度与浑浊度的反演算法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110258636.6A CN112906284A (zh) | 2021-03-09 | 2021-03-09 | 一种浑浊水域气溶胶光学厚度与浑浊度的反演算法 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| CN112906284A true CN112906284A (zh) | 2021-06-04 |
Family
ID=76107077
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110258636.6A Pending CN112906284A (zh) | 2021-03-09 | 2021-03-09 | 一种浑浊水域气溶胶光学厚度与浑浊度的反演算法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112906284A (zh) |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20080067402A (ko) * | 2007-01-16 | 2008-07-21 | 연세대학교 산학협력단 | 단일 가시광 채널을 갖는 정지궤도 인공위성을 이용한에어러솔 광학깊이 산출방법 |
| CN101329173A (zh) * | 2008-07-07 | 2008-12-24 | 武汉大学 | 一种浑浊水体大气校正方法 |
| CN102073792A (zh) * | 2011-01-19 | 2011-05-25 | 安徽师范大学 | 一种利用modis图像反演海岸带气溶胶光学特性的方法 |
| CN104019753A (zh) * | 2014-06-17 | 2014-09-03 | 杭州电子科技大学 | 一种基于modis数据反演城市大气气溶胶光学厚度的方法 |
| US20170294011A1 (en) * | 2016-04-08 | 2017-10-12 | University Of Electronic Science And Technology Of China | Method for Retrieving Atmospheric Aerosol Based on Statistical Segmentation |
| CN110186823A (zh) * | 2019-06-26 | 2019-08-30 | 中国科学院遥感与数字地球研究所 | 一种气溶胶光学厚度反演方法 |
-
2021
- 2021-03-09 CN CN202110258636.6A patent/CN112906284A/zh active Pending
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| KR20080067402A (ko) * | 2007-01-16 | 2008-07-21 | 연세대학교 산학협력단 | 단일 가시광 채널을 갖는 정지궤도 인공위성을 이용한에어러솔 광학깊이 산출방법 |
| CN101329173A (zh) * | 2008-07-07 | 2008-12-24 | 武汉大学 | 一种浑浊水体大气校正方法 |
| CN102073792A (zh) * | 2011-01-19 | 2011-05-25 | 安徽师范大学 | 一种利用modis图像反演海岸带气溶胶光学特性的方法 |
| CN104019753A (zh) * | 2014-06-17 | 2014-09-03 | 杭州电子科技大学 | 一种基于modis数据反演城市大气气溶胶光学厚度的方法 |
| US20170294011A1 (en) * | 2016-04-08 | 2017-10-12 | University Of Electronic Science And Technology Of China | Method for Retrieving Atmospheric Aerosol Based on Statistical Segmentation |
| CN110186823A (zh) * | 2019-06-26 | 2019-08-30 | 中国科学院遥感与数字地球研究所 | 一种气溶胶光学厚度反演方法 |
Non-Patent Citations (3)
| Title |
|---|
| YI WANG 等: "MODIS Retrieval of Aerosol Optical Depth over Turbid Coastal Water" * |
| 徐梦溪;许宝华;郑胜男;刘翔龙;石爱业;: "基于MODIS卫星遥感数据的大气气溶胶光学厚度优选反演方法" * |
| 王毅;石汉青;黄思训;: "中国东南地区及近海海域气溶胶反演遥感研究" * |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Zhang et al. | A Landsat 8 OLI-based, semianalytical model for estimating the total suspended matter concentration in the slightly turbid Xin’anjiang Reservoir (China) | |
| Zhang et al. | Retrieval of total suspended matter concentration in the Yellow and East China Seas from MODIS imagery | |
| Lafon et al. | SPOT shallow water bathymetry of a moderately turbid tidal inlet based on field measurements | |
| Mao et al. | A regional remote sensing algorithm for total suspended matter in the East China Sea | |
| Doxaran et al. | Spectral signature of highly turbid waters: Application with SPOT data to quantify suspended particulate matter concentrations | |
| Zhu et al. | Inversion of chromophoric dissolved organic matter from EO-1 Hyperion imagery for turbid estuarine and coastal waters | |
| Hakvoort et al. | Towards airborne remote sensing of water quality in The Netherlands—validation and error analysis | |
| Forbrich et al. | Cross-evaluation of measurements of peatland methane emissions on microform and ecosystem scales using high-resolution landcover classification and source weight modelling | |
| Vos et al. | Multiplatform optical monitoring of eutrophication in temporally and spatially variable lakes | |
| CN104406686B (zh) | 复杂地形条件下太阳短波入射辐射估算方法 | |
| Kearney et al. | The effects of tidal inundation on the reflectance characteristics of coastal marsh vegetation | |
| Wong et al. | Modeling of suspended solids and sea surface salinity in Hong Kong using Aqua/MODIS satellite images | |
| Jiao et al. | Estimation of chlorophyll‐a concentration in Lake Tai, China using in situ hyperspectral data | |
| CN113553907A (zh) | 一种基于遥感技术的森林生态环境状况评价方法 | |
| CN113639716A (zh) | 一种基于深度残差收缩网络的水深遥感反演方法 | |
| Li et al. | Evaluation of the Quasi-Analytical Algorithm (QAA) for estimating total absorption coefficient of turbid inland Waters in Northeast China | |
| CN114297938B (zh) | 一种基于神经网络的光学浅水水底深度的反演方法 | |
| CN116087115A (zh) | 一种复杂环境下的河湖水质快速检测方法 | |
| Sterckx et al. | Retrieval of suspended sediment from advanced hyperspectral sensor data in the Scheldt estuary at different stages in the tidal cycle | |
| Froidefond et al. | Method for the quantification of suspended sediments from AVHRR NOAA-11 satellite data | |
| Li et al. | Analyzing the distribution and variation of Suspended Particulate Matter (SPM) in the Yellow River Estuary (YRE) using Landsat 8 OLI | |
| CN111650128B (zh) | 一种基于地表反射率库的高分辨率大气气溶胶反演方法 | |
| CN112906284A (zh) | 一种浑浊水域气溶胶光学厚度与浑浊度的反演算法 | |
| Tester et al. | Airborne detection of ecosystem responses to an extreme event: Phytoplankton displacement and abundance after hurricane induced flooding in the Pamlico-Albemarle Sound system, North Carolina | |
| De Cauwer et al. | Optical remote sensing in support of eutrophication monitoring in the southern North Sea. |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| WD01 | Invention patent application deemed withdrawn after publication | ||
| WD01 | Invention patent application deemed withdrawn after publication |
Application publication date: 20210604 |