CN113642191B - 一种基于短波红外的遥感蒸散模型构建方法 - Google Patents
一种基于短波红外的遥感蒸散模型构建方法 Download PDFInfo
- Publication number
- CN113642191B CN113642191B CN202110982239.3A CN202110982239A CN113642191B CN 113642191 B CN113642191 B CN 113642191B CN 202110982239 A CN202110982239 A CN 202110982239A CN 113642191 B CN113642191 B CN 113642191B
- Authority
- CN
- China
- Prior art keywords
- vegetation
- remote sensing
- index
- evapotranspiration
- soil
- 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.)
- Expired - Fee Related
Links
- 238000010276 construction Methods 0.000 title claims description 9
- 239000002689 soil Substances 0.000 claims abstract description 74
- 238000000034 method Methods 0.000 claims abstract description 69
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 claims abstract description 49
- 230000003287 optical effect Effects 0.000 claims abstract description 22
- 238000011160 research Methods 0.000 claims description 33
- 238000004364 calculation method Methods 0.000 claims description 28
- 230000005855 radiation Effects 0.000 claims description 28
- 238000001704 evaporation Methods 0.000 claims description 26
- 230000008020 evaporation Effects 0.000 claims description 26
- 238000002310 reflectometry Methods 0.000 claims description 12
- 230000005068 transpiration Effects 0.000 claims description 12
- 230000001932 seasonal effect Effects 0.000 claims description 7
- 238000012937 correction Methods 0.000 claims description 6
- 238000001556 precipitation Methods 0.000 claims description 6
- 239000000470 constituent Substances 0.000 claims description 5
- 238000012545 processing Methods 0.000 claims description 5
- 230000004907 flux Effects 0.000 claims description 4
- 230000014759 maintenance of location Effects 0.000 claims description 4
- 230000008033 biological extinction Effects 0.000 claims description 3
- 238000002329 infrared spectrum Methods 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 3
- 238000013480 data collection Methods 0.000 claims description 2
- LMDZBCPBFSXMTL-UHFFFAOYSA-N 1-ethyl-3-(3-dimethylaminopropyl)carbodiimide Chemical compound CCN=C=NCCCN(C)C LMDZBCPBFSXMTL-UHFFFAOYSA-N 0.000 claims 3
- 229920006395 saturated elastomer Polymers 0.000 claims 1
- 230000008859 change Effects 0.000 abstract description 7
- 238000011156 evaluation Methods 0.000 abstract description 7
- 230000000694 effects Effects 0.000 abstract description 3
- 230000002262 irrigation Effects 0.000 abstract 2
- 238000003973 irrigation Methods 0.000 abstract 2
- 230000007547 defect Effects 0.000 description 6
- 230000006872 improvement Effects 0.000 description 4
- 239000000919 ceramic Substances 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 238000004422 calculation algorithm Methods 0.000 description 2
- 238000012512 characterization method Methods 0.000 description 2
- 238000005094 computer simulation Methods 0.000 description 2
- 239000010410 layer Substances 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 239000002356 single layer Substances 0.000 description 2
- 238000000859 sublimation Methods 0.000 description 2
- 230000008022 sublimation Effects 0.000 description 2
- 238000009834 vaporization Methods 0.000 description 2
- 230000008016 vaporization Effects 0.000 description 2
- 238000010521 absorption reaction Methods 0.000 description 1
- 235000013405 beer Nutrition 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 238000006243 chemical reaction Methods 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 238000005194 fractionation Methods 0.000 description 1
- 230000031877 prophase Effects 0.000 description 1
- 238000004088 simulation Methods 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 238000012360 testing method Methods 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
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
本发明提出一种基于短波红外的遥感蒸散模型构建方法，特别是结合遥感短波红外波段信息构建了可刻画下垫面水分变化特征的光学指数，基于提出的光学指数构建了区域蒸散发的参数化模型。本发明提出的遥感蒸散模型构建方法可改进传统蒸散发参数化模型对下垫面土壤水分状况考虑不足的问题，可提高干旱区蒸散发模型在干旱区的估算精度。本发明可服务于农业灌区管理、农业用水效率评估、区域生态环境修复、农业节水、区域人类取用水活动影响评估、灌区现代化管理等方面的应用领域。
Description
技术领域
本发明属于定量遥感和水资源基础研究技术领域，具体为一种基于短波红外遥感蒸散模型的构建方法，特别是结合遥感短波红外波段信息构建了可刻画下垫面水分变化特征的光学指数，基于提出的光学指数构建了区域蒸散发的参数化模型。本发明可服务于农业灌区水资源管理、农业用水效率评估、区域生态环境修复、农业节水、区域人类取用水活动影响评估、灌区现代化管理等方面的应用领域。
背景技术
卫星遥感是获取非均匀下垫面辐射和能量平衡信息的重要手段，被广泛用于区域蒸散发的反演和模拟研究之中。随着卫星遥感技术的进步和定量遥感模型的发展，基于地表能量平衡的遥感蒸散模型和基于地表导度参数化的蒸散模型被广泛应用于区域蒸散的估算研究之中。(1)基于地表能量平衡的遥感蒸散模型主要包括单层和双层模型。剩余阻抗参数化的可靠性和适用性仍然是单层模型应用的最大障碍，目前仍然没有适用于不同条件的剩余阻抗参数化方法。双层模型方法涉及复杂的阻抗网络计算，阻抗网络不合理的参数化常会导致蒸散估算存在较大误差。此外，由于基于地表能量平衡模型方法需结合卫星过境瞬时的地表温度，估算结果为瞬时蒸散发值。因此，基于地表能量平衡的遥感蒸散模型还需要结合时间尺度扩展方法进行日、月尺度的扩展。基于地表能量平衡的遥感蒸散模型的应用主要受制于地表温度的可获得性。由于受到云、雾等天气因素的影响，地表温度的获取往往存在困难，在一定程度限制了此类模型方法的应用。(2)基于地表导度参数化的蒸散模型得到了广泛的应用，全球和区域蒸散产品被相关单位生产和发布。但相关研究表明地表导度参数化模型对蒸散发的估算存在较大偏差，特别是土壤蒸发的估算误差较大，这与模型方法对土壤水分控制因子和土壤蒸发的刻画不足密切相关。土壤蒸发的估算与土壤水分动态密切相关。目前，基于地表导度参数化的蒸散模型土壤水分控制因子的参数化通常使用结合空气相对湿度的经验性参数化方法、前期累计降水量法和土壤水分动态模拟法。但这些方法在刻画干旱区的土壤水分动态存在明显不足。虽然基于地表导度参数化的蒸散模型取得了显著进展，但在实现极端干旱区蒸散估算方面仍然有待改进。由于不能有效的刻画下垫面土壤水分的动态变化，特别是干旱区季节性河流的土壤水分动态变化特征，地表导度参数化蒸散模型在估算干旱区的蒸散发上存在明显不足。
综上，目前，基于地表能量平衡的遥感蒸散模型受到云、雾等天气因素的影响，地表温度的获取往往存在困难，在一定程度限制了此类模型方法的应用。基于地表导度参数化的蒸散模型土壤水分控制因子的参数化方法在刻画干旱区的土壤水分动态存在明显不足。
发明内容
针对上述技术缺陷，特别是目前基于遥感的区域蒸散发估算存在的严重不足的问题，本发明提出一种基于短波红外的遥感蒸散模型构建方法，结合遥感短波红外波段信息构建了可刻画下垫面水分变化特征的光学指数，基于提出的光学指数构建了区域蒸散发的参数化模型。本发明的目的是通过以下技术方案实现的：
一种基于短波红外的遥感蒸散模型构建方法，包括以下步骤：
步骤1、研究区基础资料、气象驱动及土地利用数据收集：包括研究区的地理位置及边界、研究区的主要农作物种植情况及主要农作物的物候期、研究区内的气象站网观测资料或中国气象局陆面数据同化系统驱动数据集；
步骤2、研究区卫星遥感数据收集与处理：根据研究需要收集研究区光学卫星遥感数据，所述光学卫星遥感数据包含短波红外波段，短波红外波谱范围为1.3-3.0μm；光学遥感数据预处理包括辐射定标和大气校正；
步骤3、研究区遥感地表参数的计算：包括地表反照率、植被指数、全球植被水分指数、叶面积指数、植被覆盖度的计算；
步骤4、遥感地表参数时空重建：采用线性内插法进行时间序列重建，使用时空重建方法可对地表反照率、植被指数、全球植被水分指数、叶面积指数、植被覆盖度进行重建，生成时空连续的区域地表参数；
步骤5、气象要素场生成：基于步骤1收集的研究区气象站网观测资料或中国气象局陆面数据同化系统驱动数据集，进一步基于最近邻方法插值法生成和遥感数据分辨率一致的气象驱动数据；
步骤6、植被高度计算：植被高度计算结合重建后的NDVI进行计算：
h＝hmin+(hmax－hmin)/(NDVIveg－NDVIsoil)*(NDVI－NDVI_soil)
其中，h为植被高度，m；hmin和hmax分别为植被生长周期内的高度的最小值和最大值，NDVIveg和NDVIsoil分别为植被全覆盖条件下和裸土条件下的NDVI值，NDVI为归一化植被指数；植被高度用于后续步骤8中植被组分潜在蒸腾量计算中冠层空气动力学导度的计算；
步骤7、基于短波红外的下垫面水分供给指数计算：采用全球植被水分指数构建刻画下垫面土壤水分供给状况的短波红外指数WSCI：
gvmimax＝max(gvmimax,seasonal,0.625)
其中，gvmi为全球植被水分指数，gvmimax为全球植被水分指数的最大值；gvmimax,seasonal为像元点生长季内gvmi的最大值；
步骤8、蒸散模型参数化计算：包括净辐射Rn计算、能量分配计算、土壤组分潜在蒸发量Esp计算、植被组分潜在蒸腾量Ecp计算及植被截留蒸发量Ei计算；
步骤9、区域地表蒸散发计算：当下垫面为非冰雪和非水域覆盖时，土壤组分的蒸发量Esoil(mm/d)计算为：
Esoil＝WSCI×Esp
植被组分的蒸腾量Eveg(mm/d)计算为：
Eveg＝WSCI×(Ecp-Ei)
地表蒸散发ET(mm/d)计算为：
ET＝Ei+Eveg+Esoil。
进一步的，步骤1中的气象站网观测资料或中国气象局陆面数据同化系统驱动数据集包括的气象要素为温度、压强、湿度、风速和辐射。
进一步的，步骤2中光学卫星遥感数据可选择MODIS、FY-3、Landsat8、Sentinel2的卫星数据。
进一步的，步骤3中，基于Sentinel 2的地表反照率计算如下式：
α＝0.2266*B2+0.1236*B3+0.1573*B4+0.3417*B8+0.1170*B11+0.0338*B12
其中，B2，B3，B4，B8，B11，B12分别为哨兵波段2、3、4、8、11和12的波段反射率；
基于Sentinel 2的归一化植被指数计算如下式：
其中，B4，B8为哨兵波段4和8的波段反射率；
基于Sentinel 2的全球植被水分指数计算如下式：
其中，gvmi为全球植被水分指数；NIR为近红外波段的反射率，对应于Sentinel2的第12波段；SWIR为短波红外的反射率，对应于Sentinel2的第8波段；
基于Sentinel 2的叶面积指数计算如下式：
其中，LAI为叶面积指数，NDVI为归一化植被指数；
植被覆盖度计算如下式：
其中，fc为植被覆盖度；NDVImax和NDVImin为研究区内NDVI的最大值和最小值。
进一步的，步骤4中，采用线性内插法进行时间序列重建的具体方法为：对于每个像元的遥感地表参数构成的一维时间序列，根据数据序列中需要重建点在固定窗口内的数据线性插值进行数值估计：
Xi＝X0+k(ti-t0)
其中，Xi为需要重建的遥感地表参数，X0和X1为选定时间窗口内t0和t1时刻的遥感地表参数，k为时间序列窗口的斜率。
进一步的，步骤8中，具体的计算方案如下：
净辐射Rn的计算：
Rn＝(1-α)Rs+(εa-εs)σ(273.15+T)4
εa＝1-0.26exp(-7.77×10-4T2)
其中，α为地表反照率，Rs为短波下行辐射(W/m2)，εs为地表比辐射率，εa为大气比辐射率，T为空气温度(℃)，σ为斯特凡波尔兹曼常数；
能量分配计算：
Rns＝Rnexp(-kaLAI)
Rnc＝Rn-Rns
其中，Rn为净辐射(W/m2)，Rnc为冠层组分的净辐射(W/m2)，Rns为土壤组分的净辐射(W/m2)，ka为消光系数，LAI为叶面积指数；
土壤组分潜在蒸发量Esp(mm/d)计算：
其中，Δ为平均气温时饱和水汽压与温度相关曲线的斜率(kPa/℃)，γ为湿度计常数(kPa/℃)，Rns为土壤组分的净辐射(W/m2)。G为土壤热通量(W/m2)；
植被组分潜在蒸腾量Ecp(mm/d)计算：
其中，ρ为空气密度，ρ为空气密度(kg/m3)，Cp为空气的比热容(J/(kg·K))，es为饱和水汽压(kPa)，ea为实际水汽压(kPa)，Ga为冠层空气动力学导度(m/s)，Gc为冠层水汽导度(m/s)；
其中，zm为风速观测高度(m)，zh为湿度观测高度(m)，d为零平面取代高度(m)，zom为动量粗糙长度(m)，zoh为热力和水汽粗糙长度(m)，k为冯卡门常数，uz为风速度(m/s)；
d＝0.67*h
zom＝0.123*h
zoh＝0.1*zom
其中，h为植被高度(m)；
Gc＝gs,maxRhLAI
其中，LAI为叶面积指数，Rh为相对湿度(％)，gs,max为最大气孔导度(mm s-1)；
植被截留蒸发量计算：
Ei＝min(Di,d,Pd)
Di,d＝min(Smax,Ep)
Smax＝0.935+0.498LAI-0.00575LAI2
其中，Ei为截留蒸发量(mm/d)，Pd为日降水量(mm/d)，Di,d为截留阈值(mm/d)，Smax为最大截留量(mm/d)，Ep为日潜在蒸散量(mm/d)。
本发明的有益效果：
本发明提出的蒸散模型方法使用短波红外，由于短波热红外往往具有比热红外波段更高的空间分辨率，且短波热红外的获取相比热红外遥感更为便捷。本发明提出的遥感蒸散模型方法可改进传统蒸散发参数化模型对下垫面土壤水分状况考虑不足的问题，可提高干旱区蒸散发模型在干旱区的估算精度。因此本发明提出的模型方法更具可操作性。
附图说明
图1是本发明的方法流程图；
图2为本发明提出的蒸散模型方法与华北平原馆陶站观测数据的对比；
图3为本发明方法与PML_V2蒸散模型(PML_V2)在馆陶站的对比；
图4为本发明方法应用于和田河流域蒸散发的估算应用。
具体实施方式
本发明的基于短波红外的遥感蒸散模型构建方法，特别是结合遥感短波红外波段信息构建了可刻画下垫面水分变化特征的光学指数，基于提出的光学指数构建了区域蒸散发的参数化模型。本发明提出的遥感蒸散模型构建方法可改进传统蒸散发参数化模型对下垫面土壤水分状况考虑不足的问题，可提高干旱区蒸散发模型在干旱区的估算精度。本发明提出蒸散发模型利用对植被和土壤水分含量非常敏感的短波红外波段刻画下垫面水分供给信息的变化特征，通过结合不同时空分辨率的卫星遥感数据和高精度气象驱动场，可满足蒸散估算的不同需要，具有应用的灵活性。本发明可服务于农业灌区水资源管理、农业用水效率评估、区域生态环境修复、农业节水、区域人类取用水活动影响评估、灌区现代化管理等方面的应用领域。此外，本发明提出的蒸散模型使用短波红外，由于短波热红外往往具有比热红外波段更高的空间分辨率，且短波热红外的获取相比热红外遥感更为便捷。因此本发明提出的模型方法更具可操作性。
下面结合附图对本发明的具体实施方式以及工作原理作进一步详细说明。此处所描述的具体实施例仅用以解释本发明，并不用于限定本发明。
实施例1
下面结合附图对本发明的技术方案做进一步的详细说明，本实施例的流程如：
(一)研究区基础资料、气象驱动及土地利用数据收集
根据研究所需确定研究区的地理位置及边界。收集研究区的主要农作物种植情况及主要农作物的物候期等资料。收集研究区内的气象站网观测资料，包括温度、压强、湿度、风速和辐射等气象要素。如研究区内气象资料收集困难，可以使用CLDAS气象驱动场代替。中国气象局陆面数据同化系统(CLDAS)0.0625×0.0625°分辨率的大气驱动数据。气象驱动产品包括逐小时的2米气温、2米比湿、10米风速、地面气压、降水、短波辐射六个要素，可满足本实施例的模型驱动的需求。
(二)研究区卫星遥感数据收集与处理
根据研究需要收集研究区光学卫星遥感数据。由于本发明使用对土壤和植被含水量非常敏感的短波红外信息做支撑，因此光学卫星数据需包含短波红外波段。短波红外波谱范围为1.3-3.0μm。光学卫星数据可选择MODIS、FY-3、Landsat8、Sentinel 2等。其中，MODIS、FY-3卫星数据可满足500m至1km分辨率的区域蒸散估算需求，Landsat8、Sentinel 2可满足区域20-30米分辨率蒸散估算需求。本实施例后续论述将以Sentinel2为例。
光学遥感数据预处理包括辐射定标和大气校正。辐射定标是将传感器接收到的遥感数据，转换为辐射亮度或天顶反射率，根据所获取数据的传感器类型和相应的定标系数，对光学遥感数据进行辐射定标计算。大气校正处理旨在于消除大气的吸收、散射等效应对地表反射率的影响，消除由大气影响所造成的辐射误差，计算地物的反射率。可采用ENVI/FLAASH对收集的遥感数据进行大气校正处理。依据ENVI/FLAASH输入要求和卫星遥感数据获取的头文件信息对遥感数据进行大气校正处理，获得地表反射率数据。
由于美国地质调查局和欧空局已经开始提供Landsat8和Sentinel 2的地表反射率数据产品，可根据研究所需直接下载使用。本实施例采用Sentinel2遥感数据，可选择ESA提供的共享网站获得相应的地表反射率数据。
(三)研究区地表参数计算
本发明的技术方法所需的遥感地表参数包括：地表反照率、植被指数、叶面积指数、覆盖度、全球植被水分指数。
基于Sentinel 2的地表反照率参考[S.Bonafoni and A.Sekertekin,"AlbedoRetrieval From Sentinel-2by New Narrow-to-Broadband Conversion Coefficients,"in IEEE Geoscience and Remote Sensing Letters,vol.17,no.9,pp.1618-1622,Sept.2020,doi:10.1109/LGRS.2020.2967085.]计算如下：
α＝0.2266*B2+0.1236*B3+0.1573*B4+0.3417*B8+0.1170*B11+0.0338*B12 (1)
其中，B2，B3，B4，B8，B11，B12为哨兵波段2、3、4、8、11和12的波段反射率。
基于Sentinel 2的归一化植被指数计算如下：
其中，B4，B8为哨兵波段4和8的波段反射率。
全球植被水分指数(gvmi)由Ceccato等(2002)，基于Sentinel 2的全球植被水分指数计算如下：
其中，NIR为近红外波段的反射率，对应于Sentinel2的第12波段；SWIR为短波红外的反射率，对应于Sentinel2的第8波段。
基于Sentinel 2的叶面积指数计算如下：
其中，NDVI为归一化植被指数。
植被覆盖度为植物垂直投影面积与地域面积之比，可使用NDVI计算：
其中，NDVImax和NDVImin为研究区内NDVI的最大值和最小值。
(四)遥感地表参数时空重建
光学卫星容易受到云、雾、雨、雪等天气因素的影响，在直接获取区域高时空分辨率地表观测方面存在缺陷。遥感地表参数时空重建方法是使用最近邻插值、时序重建等方法获取时空连续的地表观测的重要手段。本发明采用线性内插法进行时间序列重建，该算法的优点是无需其他辅助数据，利用遥感地表参数时间序列时间窗口内的数据进行线性插值来填补天气因素造成的数据缺失。对于每个像元的遥感地表参数构成的一维时间序列，根据数据序列中需要重建点在固定窗口内的数据线性插值进行数值估计：
Xi＝X0+k(ti-t0) (6)
其中，X代表遥感地表参数时间序列，Xi为需要重建的遥感地表参数数据值，X0和X1为选定时间窗口内t0和t1时刻的遥感地表参数，k为时间序列窗口的斜率。使用时空重建方法可对反照率、植被指数、叶面积指数、覆盖度等遥感地表参数进行重建，生成时空连续的区域地表参数。
(五)气象要素场生成
基于步骤(一)收集的研究区气象站网资料或中国气象局陆面数据同化系统驱动数据集，进一步基于最近邻方法插值方法生成和遥感数据分辨率一致的气象驱动数据，以满足遥感蒸散模型方法的计算。本实施例采用Python SciPy的框架下使用其interpolate进行插值实现。
(六)植被高度计算
植被高度的计算结合重建后的NDVI进行计算：
h＝hmin+(hmax－hmin)/(NDVI_veg－NDVI_soil)*(NDVI－NDVI_soil) (8)
其中，hmin和hmax分别为植被生长周期内的最小值和最大值(m)，NDVI_veg和NDVI_soil分别为植被全覆盖条件下和裸土条件下的NDVI值。植被高度用于冠层空气动力学导度的计算。
(七)基于短波红外的下垫面水分供给指数计算
基于地表导度参数化的蒸散模型通过建立遥感植被指数或叶面积指数和地表导度的关系，使用Penman-Monteith，Shuttleworth-Wallace或Priestley公式来实现区域蒸散的估计。基于地表导度参数化的蒸散模型得到了广泛的应用，全球和区域蒸散产品被相关单位生产和发布。但研究表明地表导度参数化模型对蒸散发的估算存在较大偏差，特别是土壤蒸发的估算误差较大，这与模型方法对土壤水分控制因子和土壤蒸发的刻画不足密切相关。土壤蒸发的估算与土壤水分动态密切相关。目前，土壤水分控制因子的参数化通常使用结合空气相对湿度的经验性参数化方法、前期累计降水量法和土壤水分动态模拟法。但这些方法在刻画干旱区的土壤水分动态存在明显不足，特别是对干旱区的季节性河流区域，以上三种方法都难以真实准确刻画土壤水分控制。短波红外(波谱范围为1.3-3.0μm)对土壤水分和植被含水量非常敏感，基于短波红外波段信息可构建土壤水分控制因子。在这方面的初步探索工作，Yao Yunjun等(2018)和Marshall等(2020)采用MODIS短波红外波段构建了土壤水分控制因子改进了PT-JPL模型的反演精度，但改进方法仅适用于植被覆盖区域，在植被稀疏地区和裸土区域并不适用。
对此，本发明提出一种新型短波红外下垫面土壤水分供给指数的计算方法。本发明采用Ceccato等(2002)提出全球植被水分指数(global vegetation moisture index，gvmi)构建刻画下垫面土壤水分供给状况的短波红外指数(WSCI，water supply conditionindex)：
gvmimax＝max(gvmimax,seasonal,0.625) (10)
其中，gvmi为全球植被水分指数，gvmimax为全球植被水分指数的最大值；
gvmimax,seasonal为像元点生长季内gvmi的最大值。
(八)蒸散模型参数化计算
本发明的蒸散发模型参数化计算包括净辐射计算、能量分配计算、土壤组分潜在蒸发计算和植被组分潜在蒸腾计算等，具体的计算方案如下：
净辐射Rn采用Cleugh等(2007)[Cleugh H A,Leuning R,Mu Q,et al.Regionalevaporation estimates from flux tower and MODIS satellite data[J].RemoteSensing of Environment,2007,106(3):285-304.]的参数化方案：
Rn＝(1-α)Rs+(εa-εs)σ(273.15+T)4 (11)
εa＝1-0.26exp(-7.77×10-4T2) (12)
其中，α为地表反照率，Rs为短波下行辐射(W/m2)，εs为地表比辐射率，εa为大气比辐射率，T为空气温度(℃)，σ为斯特凡波尔兹曼常数(5.67×10-8W/m2/K4)。
能量分配是将净辐射分配到土壤和植被组分的过程，本发明采用比尔定律进行能量分配计算：
Rn＝Rnc+Rns (13)
Rns＝Rnexp(-kaLAI) (14)
Rnc＝Rn-Rns (15)
其中，Rn为净辐射(W/m2)，Rnc为冠层组分的净辐射(W/m2)，Rns为土壤组分的净辐射(W/m2)，ka为消光系数，LAI为叶面积指数。
土壤组分潜在蒸发Esp(mm/d)参数化：
其中，λ为气化潜热(MJ/kg),Δ为平均气温时饱和水汽压与温度相关曲线的斜率(kPa/℃)，γ为湿度计常数(kPa/℃)，Rns为土壤组分的净辐射(W/m2)。G为土壤热通量(W/m2)，参数化方案采用Mu等(2011)的土壤热通量参数化方案。[Mu Q,Zhao M,Running SW.Improvements to a MODIS global terrestrial evapotranspiration algorithm[J].Remote Sensing of Environment,2011,115(8):1781-1800.]
植被组分潜在蒸腾Ecp(mm/d)参数化：
其中，λ气化潜热(MJ/kg),Rnc为冠层组分的净辐射(W/m2)，ρ为空气密度(kg/m3)，Cp为空气的比热容(J/(kg·K))，es为饱和水汽压(kPa)，ea为实际水汽压(kPa)，Ga为冠层空气动力学导度(m/s)，Gc为冠层水汽导度(m/s)。
冠层空气动力学导度计算如下：
其中，zm为风速观测高度(m)，zh为湿度观测高度(m)，d为零平面取代高度(m)，zom为动量粗糙长度(m)，zoh为热力和水汽粗糙长度(m)，k为冯卡门常数，uz为风速度(m/s)。
d＝0.67*h (19)
zom＝0.123*h (20)
zoh＝0.1*zom (21)
其中，h为植被高度(m)。
冠层水汽导度Gc采用Yan等(2012)提出的简化方法参数化[A,H.Yan,et al."Global estimation of evapotranspiration using a leaf area index-based surfaceenergy and water balance model."Remote Sensing of Environment 124.1(2012):581-595.]：
Gc＝gs,maxRhLAI (22)
其中，LAI为叶面积指数，Rh为相对湿度(％)，gs,max为最大气孔导度(mm s-1)。Yan等(2012)认为gs,max可设定为12mm s-1。但在实际检验中，发现对不同作物类型固定gs,max会造成一定误差。因此，本发明中，采用Zhang等(2017)拟合的针对不同植被类型的最大气孔导度作为参考。具体可参考[Zhang,Y.,F.H.S.Chiew,J.
F.Sun,H.Li,andR.Leuning(2017),Global variation of transpiration and soil evaporation andthe role of their major climate drivers,J.Geophys.Res.Atmos.,122,6868–6881,doi:10.1002/2017JD027025.]。
植被截留蒸发使用Mianabadi等(2019)改进的Gerrits模型进行参数化[[1]Mianabadi A,Coenders-Gerrits M,Shirazi P,et al.Aglobal Budyko model topartition evaporation into interception and transpiration[J].Hydrology andEarth System Sciences,2019,23(12):4983-5000.]：
Ei＝min(Di,d,Pd) (23)
Di,d＝min(Smax,Ep) (24)
Smax＝0.935+0.498LAI-0.00575LAI2 (25)
其中,Ei为截留蒸发量(mm/d)，Pd为日降水量(mm/d),Di,d为截留阈值(mm/d)，Smax为最大截留量(mm/d)，Ep为日潜在蒸散量(mm/d)。
(九)区域地表蒸散发计算
当下垫面覆盖为水域时，地表蒸散发计算采用彭曼公式进行计算。
当下垫面覆盖为冰雪时，冰雪升华的参数化使用Kuzmin等提出计算方法：
Eice/snow＝(0.18+0.098u10)(es-ea) (26)
其中，Eice/snow为冰雪升华(mm/d)，u10为10米风速(m/s)，es为饱和水汽压(kPa)，ea为实际水汽压(kPa)。
当下垫面为非冰雪和非水域覆盖时，土壤组分的蒸发量计算(mm/d)：
植被组分的蒸腾Eveg(mm/d)计算：
Eveg＝WSCI×(Ecp-Ei) (28)
地表蒸散发ET(mm/d)计算：
ET＝Ei+Eveg+Esoil (29)
(十)区域实例验证
为对本发明方法的地表蒸散进行验证，以华北平原和和田河流域为研究区，使用遥感数据基于本发明方法估算了华北平原馆陶站和和田河的蒸散发。图2为本发明提出的蒸散模型方法与华北平原馆陶站观测数据的对比，图3为本发明方法与PML_V2蒸散模型(PML_V2)在馆陶站的对比，从图中可以看出相比PML_V2，本发明提出的模型估算结果与观测值更为吻合，说明本发明提出的模型方法精度较高，结果合理可信。图4为本发明方法应用于和田河流域蒸散发的估算应用，相比PML_V2和MOD16蒸散发模型，本发明提出的考虑下垫面供水状况的蒸散发参数化模型可以较好揭示和田河季节性河流地表蒸散的空间格局，特别是很好的揭示了下游河道及附近的蒸散分布状况，进一步说明本发明方法可以揭示干旱区季节河流的蒸散发格局。
以上实施例仅用于说明本发明，而并非对本发明的限制。应当指出，本领域技术人员在不脱离本发明的精神和范围的情况下还可以做出各种变化和变型，因此所有不脱离本发明范畴所做出的改进和修改也在本发明的保护范围之内。
Claims (6)
1.一种基于短波红外的遥感蒸散模型构建方法，其特征在于，包括以下步骤：
步骤1、研究区数据收集：包括1)研究区的地理位置及边界、2)研究区的主要农作物种植情况及主要农作物的物候期、3)研究区内的气象站网观测资料或中国气象局陆面数据同化系统驱动数据集；
步骤2、研究区卫星遥感数据收集与处理：根据研究需要收集研究区光学卫星遥感数据，所述光学卫星遥感数据包含短波红外波段，短波红外波谱范围为1.3-3.0μm；光学遥感数据预处理包括辐射定标和大气校正；
步骤3、研究区遥感地表参数的计算：包括地表反照率、植被指数、全球植被水分指数、叶面积指数、植被覆盖度的计算；
步骤4、遥感地表参数时空重建：采用线性内插法进行时间序列重建，使用时空重建方法对地表反照率、植被指数、全球植被水分指数、叶面积指数、植被覆盖度进行重建，生成时空连续的区域地表参数；
步骤5、气象要素场生成：基于步骤1收集的研究区气象站网观测资料或中国气象局陆面数据同化系统驱动数据集，进一步基于最近邻方法插值法生成和遥感数据分辨率一致的气象驱动数据；
步骤6、植被高度计算：
h＝hmin+(hmax-hmin)/(NDVIveg-NDVIsoil)*(NDVI-NDVIsoil)
其中，h为植被高度；hmin和hmax分别为植被生长周期内的高度的最小值和最大值，NDVIveg和NDVIsoil分别为植被全覆盖条件下和裸土条件下的NDVI值，NDVI为归一化植被指数；植被高度用于步骤8中蒸散模型参数化计算；
步骤7、基于短波红外的下垫面水分供给指数计算：
gvmimax＝max(gvmimax，seasonal，0.625)
其中，WSCI为短波红外指数；gvmi为全球植被水分指数，gvmimax为全球植被水分指数的最大值；gvmimax，seasonal为像元点生长季内gvmi的最大值；
步骤8、蒸散模型参数化计算：包括净辐射Rn计算、能量分配计算、土壤组分潜在蒸发量Esp计算、植被组分潜在蒸腾量Ecp计算及植被截留蒸发量Ei计算；
步骤9、区域地表蒸散发计算：当下垫面为非冰雪和非水域覆盖时，
土壤组分的蒸发量Esoil计算：
Esoil＝WSCI×Esp
植被组分的蒸腾量Eveg计算：
Eveg＝WSCI×(Ecp-Ei)
地表蒸散发量ET计算：
ET＝Ei+Eveg+Esoil
Esoil、Eveg、ET的单位均为mm/d。
2.根据权利要求1所述的基于短波红外的遥感蒸散模型构建方法，其特征在于，步骤1中的气象站网观测资料或中国气象局陆面数据同化系统驱动数据集包括的气象要素为温度、压强、湿度、风速和辐射。
3.根据权利要求1所述的基于短波红外的遥感蒸散模型构建方法，其特征在于，步骤2中光学卫星遥感数据可选择MODIS、FY-3、Landsat8、Sentinel 2的卫星数据。
4.根据权利要求1所述的基于短波红外的遥感蒸散模型构建方法，其特征在于，步骤3中，基于Sentinel 2的地表反照率计算如下式：
α＝0.2266*B2+0.1236*B3+0.1573*B4+0.3417*B8+0.1170*B11+0.0338*B12
其中，B2，B3，B4，B8，B11，B12分别为哨兵波段2、3、4、8、11和12的波段反射率；
基于Sentinel 2的归一化植被指数计算如下式：
其中，B4，B8为哨兵波段4和8的波段反射率；
基于Sentinel 2的全球植被水分指数计算如下式：
其中，NIR为近红外波段的反射率，对应于Sentinel2的第12波段；SWIR为短波红外的反射率，对应于Sentinel2的第8波段；
基于Sentinel 2的叶面积指数计算如下式：
其中，NDVI为归一化植被指数；
植被覆盖度计算如下式：
其中，NDVImax和NDVImin为研究区内NDVI的最大值和最小值。
5.根据权利要求1所述的基于短波红外的遥感蒸散模型构建方法，其特征在于，步骤4中，采用线性内插法进行时间序列重建的具体方法为：对于每个像元的遥感地表参数构成的一维时间序列，根据数据序列中需要重建点在固定窗口内的数据线性插值进行数值估计：
Xi＝X0+kx(ti-t0)
其中，Xi为需要重建的遥感地表参数，X0和X1为选定时间窗口内t0和t1时刻的遥感地表参数，kx为时间序列窗口的斜率。
6.根据权利要求1所述的基于短波红外的遥感蒸散模型构建方法，其特征在于，步骤8中，具体的计算方案如下：
Rn＝(1-α)Rs+(εa-εs)σ(273.15+T)4
εa＝1-0.26exp(-7.77×10-4T2)
其中，α为地表反照率；Rs为短波下行辐射，W/m2；εs为地表比辐射率；εa为大气比辐射率；T为空气温度，℃；σ为斯特凡波尔兹曼常数：5.67×10-8W/m2/K4；Rn为净辐射，W/m2；
能量分配计算：
Rns＝Rnexp(-kaLAI)
Rnc＝Rn-Rns
其中，Rn为净辐射，W/m2；Rnc为冠层组分的净辐射，W/m2；Rns为土壤组分的净辐射，W/m2；ka为消光系数，LAI为叶面积指数；
土壤组分潜在蒸发量Esp计算：
其中，λ为水的蒸发潜热，MI/kg；△为平均气温时饱和水汽压与温度相关曲线的斜率，kPa/℃；γ为湿度计常数，kPa/℃；Rns为土壤组分的净辐射，W/m2；G为土壤热通量，W/m2；Esp单位为mm/d；
植被组分潜在蒸腾量Ecp计算：
其中，λ为蒸发潜热，MJ/kg；ρ为空气密度，kg/m3；Cp为空气的比热容，J/kg-1·K-1；es为饱和水汽压，kPa；ea为实际水汽压，kPa；Ga为冠层空气动力学导度，m/s；Gc为冠层水汽导度，m/s；Ecp单位为mm/d；
其中，zm为风速观测高度，m；zh为湿度观测高度，m；d为零平面取代高度，m；zom为动量粗糙长度，m；zoh为热力和水汽粗糙长度，m；k为冯卡门常数；uz为风速度，m/s；
d＝0.67*h
zom＝0.123*h
zoh＝0.1*zom
其中，h为植被高度；
Gc＝gs，maxRhLAI
其中，LAI为叶面积指数；Rh为相对湿度，％；gs，max为最大气孔导度，mm s-1；植被截留蒸发量计算：
Ei＝min(Di，d，Pd)
Di，d＝min(Smax，Ep)
Smax＝0.935+0.498LAI-0.00575LAI2
其中，Ei为截留蒸发量，mm/d；Pd为日降水量，mm/d；Di，d为截留阈值，mm/d；Smax为最大截留量，mm/d；Ep为日潜在蒸散量，mm/d。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110982239.3A CN113642191B (zh) | 2021-08-25 | 2021-08-25 | 一种基于短波红外的遥感蒸散模型构建方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110982239.3A CN113642191B (zh) | 2021-08-25 | 2021-08-25 | 一种基于短波红外的遥感蒸散模型构建方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113642191A CN113642191A (zh) | 2021-11-12 |
| CN113642191B true CN113642191B (zh) | 2022-03-22 |
Family
ID=78423915
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110982239.3A Expired - Fee Related CN113642191B (zh) | 2021-08-25 | 2021-08-25 | 一种基于短波红外的遥感蒸散模型构建方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113642191B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN115797797B (zh) * | 2023-02-09 | 2023-04-21 | 水利部交通运输部国家能源局南京水利科学研究院 | 一种蒸散发塔基遥感监测方法系统及存储介质 |
Citations (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7242803B2 (en) * | 2003-01-21 | 2007-07-10 | The United States Of America As Represented By The Secretary Of The Navy | System and method for significant dust detection and enhancement of dust images over land and ocean |
| US7491944B1 (en) * | 2005-07-14 | 2009-02-17 | Sandia Corporation | Method to analyze remotely sensed spectral data |
| CN103196862A (zh) * | 2013-02-25 | 2013-07-10 | 北京师范大学 | 基于ASAR和Hyperion数据反演植被覆盖下土壤水分的方法及系统 |
| CN104200080A (zh) * | 2014-08-22 | 2014-12-10 | 中国科学院遥感与数字地球研究所 | 一种农业区水资源环境诊断方法及其系统 |
| CN105608307A (zh) * | 2015-12-09 | 2016-05-25 | 中国水利水电科学研究院 | 一种水资源干旱评价系统 |
| CN105844000A (zh) * | 2016-03-18 | 2016-08-10 | 江苏铨铨信息科技有限公司 | 一种mcc 表面海流反演方法 |
| CN106295221A (zh) * | 2016-08-19 | 2017-01-04 | 中国科学院遥感与数字地球研究所 | 一种城市水体类别提取方法和系统 |
| CN107065036A (zh) * | 2017-04-19 | 2017-08-18 | 中国农业科学院农业资源与农业区划研究所 | 一种联合遥感和气象数据获取全天候蒸散发的方法 |
| CN109187360A (zh) * | 2018-09-17 | 2019-01-11 | 中国科学院地理科学与资源研究所 | 基于反射率-植被覆盖度二维空间的蒸散发遥感反演方法和应用 |
| CN109993237A (zh) * | 2019-04-15 | 2019-07-09 | 中国水利水电科学研究院 | 基于高分卫星光学遥感数据的水体快速提取方法及系统 |
| CN111611965A (zh) * | 2020-05-29 | 2020-09-01 | 中国水利水电科学研究院 | 一种基于Sentinel-2影像的陆表水体提取方法 |
| CN111999251A (zh) * | 2020-08-14 | 2020-11-27 | 中国水利水电科学研究院 | 基于热红外遥感的区域植被蒸腾和土壤蒸发反演的遥感模型方法 |
| CN113128401A (zh) * | 2021-04-19 | 2021-07-16 | 中国水利水电科学研究院 | 基于光学和雷达遥感数据的区域实际灌溉面积监测方法 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20060081639A1 (en) * | 2004-10-14 | 2006-04-20 | Lifetime Hoan Corporation | Thermochromic cookware |
| CN102176002B (zh) * | 2010-12-30 | 2013-07-24 | 中国科学院地理科学与资源研究所 | 基于地表水热通量遥感反演的干旱监测方法及系统 |
| CN103363962B (zh) * | 2013-07-23 | 2015-04-08 | 中国科学院遥感与数字地球研究所 | 一种基于多光谱影像的湖泊水储量遥感估算方法 |
| CN111401336A (zh) * | 2020-05-09 | 2020-07-10 | 中国地质科学院岩溶地质研究所 | 基于反射率-植被覆盖度二维空间的蒸散发遥感反演应用 |
| CN112991247B (zh) * | 2021-03-04 | 2022-06-03 | 河南省气象科学研究所 | 冬小麦蒸散遥感反演及与作物模型同化方法 |
-
2021
- 2021-08-25 CN CN202110982239.3A patent/CN113642191B/zh not_active Expired - Fee Related
Patent Citations (13)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US7242803B2 (en) * | 2003-01-21 | 2007-07-10 | The United States Of America As Represented By The Secretary Of The Navy | System and method for significant dust detection and enhancement of dust images over land and ocean |
| US7491944B1 (en) * | 2005-07-14 | 2009-02-17 | Sandia Corporation | Method to analyze remotely sensed spectral data |
| CN103196862A (zh) * | 2013-02-25 | 2013-07-10 | 北京师范大学 | 基于ASAR和Hyperion数据反演植被覆盖下土壤水分的方法及系统 |
| CN104200080A (zh) * | 2014-08-22 | 2014-12-10 | 中国科学院遥感与数字地球研究所 | 一种农业区水资源环境诊断方法及其系统 |
| CN105608307A (zh) * | 2015-12-09 | 2016-05-25 | 中国水利水电科学研究院 | 一种水资源干旱评价系统 |
| CN105844000A (zh) * | 2016-03-18 | 2016-08-10 | 江苏铨铨信息科技有限公司 | 一种mcc 表面海流反演方法 |
| CN106295221A (zh) * | 2016-08-19 | 2017-01-04 | 中国科学院遥感与数字地球研究所 | 一种城市水体类别提取方法和系统 |
| CN107065036A (zh) * | 2017-04-19 | 2017-08-18 | 中国农业科学院农业资源与农业区划研究所 | 一种联合遥感和气象数据获取全天候蒸散发的方法 |
| CN109187360A (zh) * | 2018-09-17 | 2019-01-11 | 中国科学院地理科学与资源研究所 | 基于反射率-植被覆盖度二维空间的蒸散发遥感反演方法和应用 |
| CN109993237A (zh) * | 2019-04-15 | 2019-07-09 | 中国水利水电科学研究院 | 基于高分卫星光学遥感数据的水体快速提取方法及系统 |
| CN111611965A (zh) * | 2020-05-29 | 2020-09-01 | 中国水利水电科学研究院 | 一种基于Sentinel-2影像的陆表水体提取方法 |
| CN111999251A (zh) * | 2020-08-14 | 2020-11-27 | 中国水利水电科学研究院 | 基于热红外遥感的区域植被蒸腾和土壤蒸发反演的遥感模型方法 |
| CN113128401A (zh) * | 2021-04-19 | 2021-07-16 | 中国水利水电科学研究院 | 基于光学和雷达遥感数据的区域实际灌溉面积监测方法 |
Non-Patent Citations (4)
| Title |
|---|
| Contribution of understorey vegetation to evapotranspiration partitioning in apple orchards under Mediterranean climatic conditions in South Africa;Z. Ntshidi 等;《Agricultural Water Management》;20201010;第1-11页 * |
| Evolution of evapotranspiration models using thermal and shortwave remote sensing data;Jing M. Chen 等;《Remote Sensing of Environment》;20191231;第1-20页 * |
| 利用短波红外进行地基CO2反演的影响因子分析;霍彦峰 等;《光谱学与光谱分析》;20150930;第2500-2504页 * |
| 区域蒸散发遥感估算方法及验证综述;张荣华 等;《地球科学进展》;20121231;第1296-1307页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN113642191A (zh) | 2021-11-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Minacapilli et al. | A time domain triangle method approach to estimate actual evapotranspiration: Application in a Mediterranean region using MODIS and MSG-SEVIRI products | |
| Vinukollu et al. | Global estimates of evapotranspiration for climate studies using multi-sensor remote sensing data: Evaluation of three process-based approaches | |
| Boutin et al. | Sea surface freshening inferred from SMOS and ARGO salinity: Impact of rain | |
| Suzuki et al. | Interannual covariability between actual evapotranspiration and PAL and GIMMS NDVIs of northern Asia | |
| Shu et al. | Estimation of regional evapotranspiration over the North China Plain using geostationary satellite data | |
| Schlenz et al. | Uncertainty assessment of the SMOS validation in the upper Danube catchment | |
| Chandrapala et al. | Satellite measurements supplemented with meteorological data to operationally estimate evaporation in Sri Lanka | |
| Peng et al. | Estimation of evapotranspiration from MODIS TOA radiances in the Poyang Lake basin, China | |
| Song et al. | A revised surface resistance parameterisation for estimating latent heat flux from remotely sensed data | |
| Yeh et al. | Determining the precipitable water vapor thresholds under different rainfall strengths in Taiwan | |
| El-Shirbeny et al. | Assessment of the mutual impact between climate and vegetation cover using NOAA-AVHRR and Landsat data in Egypt | |
| Wu et al. | An improved method for deriving daily evapotranspiration estimates from satellite estimates on cloud-free days | |
| Mallick et al. | Latent heat flux estimation in clear sky days over Indian agroecosystems using noontime satellite remote sensing data | |
| Zheng et al. | Global evapotranspiration derived by ETMonitor model based on earth observations | |
| Tan et al. | A method for downscaling daily evapotranspiration based on 30-m surface resistance | |
| CN113642191B (zh) | 一种基于短波红外的遥感蒸散模型构建方法 | |
| Khesali et al. | A method in near-surface estimation of air temperature (NEAT) in times following the satellite passing time using MODIS images | |
| CN103678884A (zh) | 一种基于hj卫星的地表实际蒸散发动态监测方法 | |
| Ferreira et al. | Assessment of water demands for irrigation using energy balance and satellite data fusion models in cloud computing: A study in the Brazilian semiarid region | |
| Liu et al. | Actual daily evapotranspiration estimated from MERIS and AATSR data over the Chinese Loess Plateau | |
| Hoeltgebaum et al. | Evaluation of the storage and evapotranspiration terms of the water budget for an agricultural watershed using local and remote-sensing measurements | |
| Song et al. | Reconstruction of remotely sensed daily evapotranspiration data in cloudy-sky conditions | |
| Ramos et al. | Assessment of satellite precipitation products in the Philippine archipelago | |
| Zhou et al. | Estimation of Hourly Air Temperature in China Based on LightGBM and Himawari-8 | |
| Chávez et al. | Maize Evapotranspiration Estimates Using Planet Dove Mini-Satellites and Field-Level InfraRed Thermometers |
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 | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| CF01 | Termination of patent right due to non-payment of annual fee | ||
| CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20220322 |