CN106372434B - 一种被动微波遥感瞬时地表发射率估算方法及装置 - Google Patents

一种被动微波遥感瞬时地表发射率估算方法及装置 Download PDF

Info

Publication number
CN106372434B
CN106372434B CN201610800047.5A CN201610800047A CN106372434B CN 106372434 B CN106372434 B CN 106372434B CN 201610800047 A CN201610800047 A CN 201610800047A CN 106372434 B CN106372434 B CN 106372434B
Authority
CN
China
Prior art keywords
emissivity
mpdi
polarization
under
instantaneous
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.)
Active
Application number
CN201610800047.5A
Other languages
English (en)
Other versions
CN106372434A (zh
Inventor
邱玉宝
郭华东
施建成
石利娟
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Remote Sensing and Digital Earth of CAS
Original Assignee
Institute of Remote Sensing and Digital Earth of CAS
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Institute of Remote Sensing and Digital Earth of CAS filed Critical Institute of Remote Sensing and Digital Earth of CAS
Priority to CN201610800047.5A priority Critical patent/CN106372434B/zh
Publication of CN106372434A publication Critical patent/CN106372434A/zh
Application granted granted Critical
Publication of CN106372434B publication Critical patent/CN106372434B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16Z—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS, NOT OTHERWISE PROVIDED FOR
    • G16Z99/00—Subject matter not provided for in other main groups of this subclass

Landscapes

  • Radiation Pyrometers (AREA)

Abstract

本发明提供了一种被动微波遥感瞬时地表发射率估算方法及装置。该方法或装置包括：用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据，计算晴空条件下瞬时地表发射率的步骤或装置；获取不同频率水平和垂直极化瞬时地表发射率的步骤或装置；根据同频率水平和垂直极化的亮温确定微波极化差指数和极化率的步骤或装置；根据瞬时地表发射率的水平极化发射率、以及微波极化差指数和极化率建立二次多项式和线性回归方程，计算云覆盖条件下的地表发射率的水平发射率的步骤或装置。采用本发明技术方案中的方法或装置对全球进行了分区的局部分块测试，实现全球云层覆盖条件下地表发射率的估算，为地表温度反演和大气校正计算提供重要的数据基础。

Description

一种被动微波遥感瞬时地表发射率估算方法及装置
技术领域
本发明涉及遥感数据地学参数反演技术，具体涉及一种被动微波遥感瞬时地表发射率估算方法及装置。
背景技术
在对地的观测中，对于不同频率、不同极化和不同时间的地表向外发射的微波发射率的准确模拟、估算和认识理解是进一步开展微波遥感应用的基础，也是应用微波数据进行生态系统模型、辐射传输模型、气候模型、陆面模型和从卫星被动微波观测资料进行地面上大气反演的重要输入参数，特别是在大气温湿廓线反演、地表温度反演、植被变化监测、积雪特征、全球陆表降雨估算以及地表土壤水分等方面，地表发射率的特征性信息具有重要的作用。
在全球地表发射率的估算研究中，根据局部地区的大气降水、云中水汽、陆地表温度和土壤水分等研究的需要，研究局部地区地表发射率的情况较为普遍，且为多月或半月内的平均地表发射率，在大气影响考虑方面往往多采用地面探空数据和月平均大气状况等。 Jones,A.S.对平均地表发射率进行了局部地区的估算；Choudhury B.J. 对非洲的沙漠和两个植被覆盖区研究了19GHz和37GHz的地表发射率月平均计算；Prigent,C.使用SMM/I进行了全球地表发射率的月平均的估算；在Prigent,C.的方法基础之上Karbou使用AMSU微波测量进行了半月平均地表发射率的估算；Norouzi,H.在采用地面探空数据的基础之上，形成了月平均的全球陆表发射率的估算。在地表瞬时发射率估算方面，由于同步观测的大气水汽探空数据获取和瞬时大气状况的获取比较困难，个别方法集中在较小区域，如Xiang X.等提出了SMM/I 局部瞬时地表发射率的反演方法，主要是基于地表发射率的昼夜变化不显著的假设，大气的辅助数据采用了探空气球的数据。
由于大气参数中大气水汽、温度廓线以及云的识别、夜间的露水、土壤水分的变化等均会对地表发射率的估算直接产生影响，而最为直接影响地表发射率产生变化的是由于大气云层覆盖等，使得对于其瞬间地表发射率估算存在着一定的难度，在云层覆盖条件下，下垫面的信息在传输过程中受到影响，从而使得陆面参数化反演(如土壤水分、地表温度等)出现误差。
发明内容
针对现有技术中的缺陷，本发明提供一种被动微波遥感瞬时地表发射率估算方法及装置，采用星载被动微波传感器AMSR-E和中分辨率成像光谱仪MODIS的两种仪器实现同步观测地表数据，计算获得晴空条件下的瞬时地表发射率，根据瞬时地表发射率的水平极化发射率、微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程实现对云层覆盖条件下的地表瞬时发射率的估算从而形成对整个全球瞬时地表发射率的估算能力。
为实现上述目的，本发明提供以下技术方案：
一种被动微波遥感瞬时地表发射率估算方法，包括如下步骤：
使用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据，根据所述地表数据计算晴空条件下的瞬时地表发射率；
获取所述晴空条件下的不同频率水平和垂直极化瞬时地表发射率；
根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR；
根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程；
采用所述二次多项式和线性回归方程计算云层覆盖条件下的地表发射率的水平发射率。
根据所述地表数据计算晴空条件下的瞬时地表发射率的步骤包括根据如下公式计算获取晴空条件下的瞬时地表发射率的步骤：
Figure BDA0001108842260000031
其中，Tbp为卫星观测亮温，v为频率，θ为辐射计观测角度，Ts为地表温度，Γ(v,θ)为大气的吸收因子，Tatm↓(v,θ)和Tatm↑(v,θ)为向下或向上的大气辐射，TCB为宇宙背景辐射。
所述根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR的步骤包括根据如下公式计算微波极化差指数MPDI和极化率PR的步骤：
PR＝Tbh/Tbv
MPDI＝(Tbv-Tbh)/(Tbv+Tbh)
其中，Tbh和Tbv为相同频率下的水平极化和垂直极化的亮温。
根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程的步骤包括如下步骤：
eh＝B0+B1·PR/MPDI+B2·PR/MPDI2
eh1＝B0+B1·PR/MPDI6.9+B2·PR/MPDI10.7+B3·PR/MPDI18.7
其中，eh为晴空条件下的地表发射率的水平极化发射率，eh1云层覆盖条件下的地表发射率水平发射率，MPDI为微波极化差指数， MPDI6.9,MPDI10.7,MPDI18.7分别为在频率为C波段，X波段，Ku波段下的微波极化差指数，B0,B1,B2,B3分别为晴空条件下的局部地区地表参数，其中晴空条件下的局部地区地表参数与云层覆盖条件下的该局部地区地表参数相同。
所述亮温Tbp采用AMSR-E L2A双极化升/降轨道的亮度温度数据，所述亮度温度数据采用算术平均方法转换到同一投影格式；
大气参数采用MODIS/Aqua 5-Min Level2的20层大气廓线数据，所述大气廓线数据采用算术平均方法转换到同一投影格式；
地表温度数据采用MODIS/Aqua Land Surface Temperature/Emissivity 5-MinL2Swath-1km测量的数据，所述测量的数据采用算术平均方法转换到同一投影格式，测量数据中的质量控制数据QA作为晴空条件的判别原则。
另一方面，本发明提供了一种被动微波遥感瞬时地表发射率估算的装置，所述装置包括：
用于使用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据、根据所述地表数据计算晴空条件下的瞬时地表发射率的装置；
用于获取所述晴空条件下的不同频率水平和垂直极化瞬时地表发射率的装置；
用于根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR的装置；
用于根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程的装置；
用于所述采用二次多项式和线性回归方程计算云层覆盖条件下的地表发射率的水平发射率的装置。
用于使用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据、根据所述地表数据计算晴空条件下的瞬时地表发射率的装置包括根据如下公式计算获取晴空条件下的瞬时地表发射率的装置：
Figure BDA0001108842260000041
其中，Tbp为卫星观测亮温，v为频率，θ为辐射计观测角度，Ts为地表温度，Γ(v,θ)为大气的吸收因子，Tatm↓(v,θ)和Tatm↑(v,θ)为向下或向上的大气辐射，TCB为宇宙背景辐射。
用于根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR的装置包括根据如下公式计算微波极化差指数MPDI和极化率PR的装置：
PR＝Tbh/Tbv
MPDI＝(Tbv-Tbh)/(Tbv+Tbh)
其中，Tbh和Tbv为相同频率下的水平极化和垂直极化的亮温。
用于根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程的装置包括根据如下公式计算云层覆盖条件下的地表发射率水平发射率方程的装置：
eh＝B0+B1·PR/MPDI+B2·PR/MPDI2
eh1＝B0+B1·PR/MPDI6.9+B2·PR/MPDI10.7+B3·PR/MPDI18.7
其中，eh为晴空条件下的地表发射率的水平极化发射率，eh1云层覆盖条件下的地表发射率水平发射率，MPDI为微波极化差指数， MPDI6.9,MPDI10.7,MPDI18.7分别为在频率为C波段，X波段，Ku波段下的微波极化差指数，B0,B1,B2,B3分别为晴空条件下的局部地区地表参数，其中晴空条件下的局部地区地表参数与云层覆盖条件下的该局部地区地表参数相同。
所述亮温Tbp采用装置AMSR-E L2A双极化升/降轨道的亮度温度数据；
大气参数采用装置MODIS/Aqua 5-Min Level2的20层大气廓线数据；
地表温度数据采用装置MODIS/Aqua Land Surface Temperature/Emissivity 5-Min L2Swath-1km测量的数据。
由上述技术方案可知，本发明所述的一种被动微波遥感瞬时地表发射率估算方法及装置，采用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据；采集数据速度快、精准以及图像分辨率高，更加适用于计算晴空条件下的瞬时地表发射率，根据晴空条件下的瞬时地表发射率的水平极化发射率、微波极化差指数MPDI和极化率PR 建立二次多项式和线性回归方程来实现全球云层覆盖覆盖条件下地表发射率的估算。为地表温度反演和大气校正计算提供重要的数据基础，也为微波数据进行生态系统模型、辐射传输模型、气候模型、陆面模型和从卫星被动微波观测资料进行地面上大气反演提供重要精准的输入参数。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1是本发明实施例提供的被动微波遥感瞬时地表发射率估算方法的流程图；
图2是本发明实施例中步骤S101的具体流程图；
图3是本发明实施例中PR和MPDI大气敏感性模拟分析曲线图；
图4是本发明实施例中步骤S103的具体流程图；
图5是本发明实施例中计算发射率和通过方程估算发射率的散点图；
图6是本本发明实施例中瞬时发射率统计直方图；
图7是本本发明实施例中半月平均发射率统计直方图；
图8是本本发明实施例中预测/晴空瞬时发射率和14天平均水平发射率差的统计直方图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在对地的观测中，地表发射率的特征性信息在大气温湿廓线反演、地表温度反演、植被变化监测、积雪特征、全球陆表降雨估算以及地表土壤水分等方面具有重要的作用。由于大气参数中大气水汽、温度廓线以及云的识别、夜间的露水、土壤水分的变化等均会对地表发射率的估算直接产生影响，而最为直接影响地表发射率产生变化的是由于大气云层覆盖等，使得对于其瞬间地表发射率估算存在着一定的难度，在云层覆盖条件下，下垫面的信息在传输过程中受到影响，从而使得陆面参数化反演出现误差。为了解决该问题，本发明下述实施例提出了一种被动微波遥感瞬时地表发射率估算方法及装置。
图1示出了本发明实施例提供的被动微波遥感瞬时地表发射率估算方法的流程图，参见图1，本发明实施例提供的被动微波遥感瞬时地表发射率估算方法包括如下步骤：
S101：使用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据，根据所述地表数据计算晴空条件下的瞬时地表发射率；
如图2所示，在本步骤中，由星载被动微波传感器接收的晴空条件下地表温度数据、晴空条件下大气参数数据和亮温数据并将三种数据经过投影变换；
亮温采用原始的AMSR-E L2A(Advanced Microwave Scanning Radiometer-EarthObserving System Level-2A product)双极化升/降轨道的亮度温度数据(Ashcroft，2006)。AMSR-E L2A全球数据每天有 14或者15轨道(升/降轨道)，在处理计算的时候由于不同频率的分辨率不同，因此我们选用分辨率最高的数据(Res.1 6.925GHz,Res.2 for10.65GHz，18.7GHz，Res.3 23.8GHz，36.5GHz和Res.4 89.0GHz)，所有的轨道产品数据使用算术平均方法转换到同一投影格式，如 EASE-GRID(Equal-Area Scalable Earth Grid)；
大气参数采用MODIS/Aqua 5-Min Level2的20层大气廓线数据 (MYD07L2)，地表温度数据采用MYD11_L2(MODIS/Aqua Land Surface Temperature/Emissivity 5-MinL2Swath-1km)产品(Wan Z, 2002)，同时该产品的质量控制数据(QA)作为晴空条件的判别原则。获得的MODIS数据为每5分钟的地表温度数据和大气分层廓线数据也被投影到相同的同一投影格式，如EASE-GRID格点中。
由于大气的厚度以及地面海拔高度和位势高度的差异，在辐射传输过程中，使用GOTOP30DEM(Digital Elevation Model)数据对地表高程进行处理。然后，剔除全球的水体(海洋和较大的内陆水体)，采用了NSIDC(National Snow and Ice Data Center)AMSR-E土壤水分算法中的水体辅助数据(阀值采用>80％)地表发射率被看为是有效的地表发射率，这主要式因为AMSR-E的视场角是由不同的地表覆盖物组成，同时忽略了微波穿透深度的不同，假设地面温度和陆地表面的辐射温度相等。
根据如下公式计算取晴空条件下的瞬时地表发射率：
Figure BDA0001108842260000081
Tbp为卫星观测亮温，频率v，极化p，辐射计观测角度θ，Ts为地表温度，Γ(v,θ)大气的吸收因子，Tatm↓(v,θ)和Ta t m↑(v,θ)为向下或向上的大气辐射，TCB宇宙背景辐射(2.7K)
吸收系数Γ(v,θ)＝e-τ(v,0,∞)·secθ为光学厚度τ(v,0,∞)的负指数关系，向上和向下的大气贡献如下
Figure BDA0001108842260000082
Figure BDA0001108842260000083
这里，T(z')在z'km高度的大气温度，γabs(z',v)在高度z'km单位长度和频率的大气吸收影响因子，光学厚度τ(0,∞,v,θ)则可以表达为
Figure BDA0001108842260000084
利用大气吸收模型计算γabs(z',v)，τ(v,0,∞)和Γ(v,θ)，利用Edington近似计算和，结合输入的亮温，地表温度等数据，可计算发射率es,p(v,θ)。
S102：获取所述晴空条件下的不同频率水平和垂直极化瞬时地表发射率；
在本步骤中，采用正交分解法计算瞬时地表发射率的水平极化瞬时地表发射率和垂直极化瞬时地表发射率；
S103：根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR；
在本步骤中，低频情况下水平极化云层覆盖条件下的地表发射率可以通过10.7GHz/18.7GHz的水平极化发射率和PR/MPDI的线性或者二次多项式的关系进行预测。
其极化率PR和微波极化差指数MPDI的表达式如下两式：
PR＝Tbh/Tbv (4.5)
MPDI＝(Tbv-Tbh)/(Tbv+Tbh) (4.6)
其中，Tbh和Tbv为同一频率下的水平和垂直极化的亮温。
基于公式(4.1)、(4.5)和(4.6)，地表温度和大气的光学厚度对于PR/MPDI的影响分析，可表达为
PR＝(Ts·eh·Γ+T↑+T↓·(1-eh)·Γ)/(Ts·eν·Γ+T↑+T↓·(1-ev)·Γ) (4.8)
MPDI＝(eν-eh)/(ev+eh+g) (4.9)
其中eν和eh为地表的V极化和H极化的发射率，
g＝2(T↑+T↓·Γ)/(Ts-T↓·Γ) (4.10)
假设大气的有效温度和地表温度相等，也即Ts＝Ta(Ta为大气的有效温度)，同时可以认为上行和下行辐射相等，并且T↑＝T↓＝Ta·(1-Γ)。公式(4.8)和(4.9)可以重写为
PR＝(1-Γ2+eh·Γ)/(1-Γ2+eν·Γ) (4.11)
MPDI＝(eν-eh)/(eν+eh+1/Γ-1) (4.12)
当大气的透过率Γ趋近于1的时候(如6.9GHz和10.7GHz情况下的大气透过率Γ趋近于1的值)，然后，可以得到
PR≈eh/ev (4.13)
MPDI≈(ev-eh)/(ev+eh) (4.14)
为分析MPDI/PR对大气影响的敏感性，把大气的透过率作为变量，使其从不透明的Γ＝0.1开始，到完全透明Γ＝1.0，如图3中的a和b。据不同的大气水汽和云水含量的发射率进行模拟，如图3中的c和d。图 3中，在给定的两种地表发射率情况下(eh＝0.85，ev＝0.95)和 (eh＝0.80，eh＝0.95)对大气透过率的模拟使用US62标准大气情况下，增加不同的大气水汽和云中液态水(云)，从标准大气(US62)发射率的模拟可以看出，对于6.9GHz和10.7GHz,针对不同的综合水汽含量，从6.5g/cm2到2.5g/cm2,其透过率的变化较小，分别为Γ6.9GHz的 0.980-0.990和Γ10.7GHz的0.960-0.977。使PR/MPDI在0.006和0.001间变化
而当云水含量从0.22g/cm2到0.07g/cm2变化时，Γ6.9GHz、Γ10.7GHz和Γ18.7GHz变化分别为0.950-0.985,0.905-0.975和0.735-0.940，从而引起 PR/MPDI的变化分别为0.0112/0.0018(6.9GHz)，0.0216/0.0036 (10.7GHz)和0.0554/0.0116(18.7GHz)。
对于PR/MPDI的这两种公式(4.13)和(4.14)是基于 10.7GHz/18.7GHz在大气穿透率相对高且大气和非散射情况下层覆盖的假设而成立。
从上面分析可以看出PR和MPDI和地表温度无关，并且对于大气影响不敏感等特点，初步分析，其包含的信息和地表发射率可以说具有一定相关性，进一步从近似公式(4.12)和(4.13)可以看出，PR 和MPDI均是由地表发射率(土壤水分、粗糙度和植被等参数)决定的。通过敏感性分析认为MPDI对于大气的影响来说具有比较好的稳定性。
S104：根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程；
如图4所示，在本步骤中，通过在晴空条件下，使用基于MODIS (ModerateResolution Imaging Spectroradiometer)的地表发射率估算的水平极化10.7GHz和18.7GHz发射率的值来建立预测方程来估算 6.9GHz、10.7GHz和18.7GHz云层覆盖条件下的水平极化发射率的值。然后结合发射率之间的相互关系确定其他频率和极化的地表发射率(6.9GHz、10.7GHz、18.7GHz和36.5GHz)。
根据分析，使用MPDI或者PR提出的可能的估算方程为式(4.15) 的二次多项式和(4.16)的多变量线性回归方程(或单变量线性回归)，
eh＝B0+B1·PR/MPDI+B2·PR/MPDI2 (4.15)
eh1＝B0+B1·PR/MPDI6.9+B2·PR/MPDI10.7+B3·PR/MPDI18.7 (4.16)
其中，eh为晴空条件下的地表发射率的水平极化发射率，eh1云层覆盖条件下的地表发射率水平发射率，MPDI为微波极化差指数， MPDI6.9,MPDI10.7,MPDI18.7分别为在频率为C波段，X波段，Ku波段下的微波极化差指数，B0,B1,B2,B3分别为晴空条件下的局部地区地表参数，由于C波段和X波段的强穿透性，其中晴空条件下的局部地区地表参数与云层覆盖条件下的该局部地区地表参数相同；
所述C波段的频率优选为6.9GHz，X波段的频率优选为10.7GHz， Ku波段的频率优选为18.7GHz。
S105：采用所述步骤S104中的二次多项式和线性回归方程计算云层覆盖条件下的地表发射率的水平发射率。
另一方面，本发明提供了一种被动微波遥感瞬时地表发射率估算的装置，包括：
用于使用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据、根据所述地表数据计算晴空条件下的瞬时地表发射率的装置；
用于获取所述晴空条件下的不同频率水平和垂直极化的瞬时地表发射率的装置；
用于根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR的装置；
用于根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程的装置；
用于采用所述二次多项式和线性回归方程计算云层覆盖条件下的地表发射率的水平发射率的装置。
用于使用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据、根据所述地表数据计算晴空条件下的瞬时地表发射率的装置包括根据如下公式计算获取晴空条件下的瞬时地表发射率的装置：
Figure BDA0001108842260000121
其中，Tbp为卫星观测亮温，v为频率，θ为辐射计观测角度，Ts为地表温度，Γ(v,θ)为大气的吸收因子，Tatm↓(v,θ)和Tatm↑(v,θ)为向下或向上的大气辐射，TCB为宇宙背景辐射。
用于根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR的装置包括根据如下公式计算微波极化差指数MPDI和极化率PR的装置：
PR＝Tbh/Tbv
MPDI＝(Tbv-Tbh)/(Tbv+Tbh)
其中，Tbh和Tbv为相同频率下的水平极化和垂直极化的亮温。
用于根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程的装置包括根据如下公式计算云层覆盖条件下的地表发射率水平发射率方程的装置：
eh＝B0+B1·PR/MPDI+B2·PR/MPDI2
eh1＝B0+B1·PR/MPDI6.9+B2·PR/MPDI10.7+B3·PR/MPDI18.7
其中，eh为晴空条件下的地表发射率的水平极化发射率，eh1云层覆盖条件下的地表发射率水平发射率，MPDI为微波极化差指数， MPDI6.9,MPDI10.7,MPDI18.7分别为在频率为C波段，X波段，Ku波段下的微波极化差指数，B0,B1,B2,B3分别为晴空条件下的局部地区地表参数，其中晴空条件下的局部地区地表参数与云层覆盖条件下的该局部地区地表参数相同。
所述亮温Tbp采用装置AMSR-E L2A双极化升/降轨道的亮度温度数据；大气参数采用装置MODIS/Aqua 5-Min Level2的20层大气廓线数据；地表温度数据采用装置MODIS/Aqua Land Surface Temperature/Emissivity 5-Min L2 Swath-1km测量的数据。
采用本发明技术方案中的装置对全球进行了分区的局部分块测试，实现全球云层覆盖条件下地表发射率的估算。为地表温度反演和大气校正计算提供重要的数据基础。
实施例一，本发明实施例选择在中国的两块地方进行10.7GHz和 18.7GHz的发射率估算性测试，测试区中地表覆盖主要是草地(青藏高原东南地区)和作物(中国中部)，采用测试的瞬时地表发射率分布图 (时间为12-08-2006升轨)，使用这两块区域来测试PR和MPDI的发射率的预测能力。数据格式均为EASE-GRID(Equal-Area Scalable EarthGrid)网格。
在2008年8月12日升轨的有效发射率采样点分别为505和157点。同时使用PR和MPDI两种极化参量建立起公式4.15(二次多项式)和4.16 (多/单变量线性回归)的发射率预测方程。方程所获得发射率和其水平极化发射率的均方根误差(RMSE)如表1所示。从表1中可以看出，对于PR和MPDI在发射率的预测能力上具有相当的水平，根据RMSE来看，二次多项式具有更高的精度，二次多项式更加能描述到在地表上具有一层稀薄植被层的非线性影响。
表1测试区中四个低频率中不同预测方程均拟合的均方根误差
Figure BDA0001108842260000131
Figure BDA0001108842260000141
选择亚洲其他森林进行进一步的研究分析，我们发现在落叶阔叶林，发射率的PR线性方程估算的RMSE在10.7GHz的时候为0.0077，而使用PR的二次多项式进行估算，其RMSE为0.00535，而且在其他的森林和植被覆盖区上使用PR的二次多项式拟合均比线性方程的拟合具有更小的RMSE。可以发现，低频(6.9GHz/10.7GHz)的发射率估算具有比较差的RMSE，这主要是由于低频对土壤水分具有比较强的反映，这导致亮温的变化性比较大。
结合大气的敏感性分析，36.5GHz以及更高的频率上，由于大气影响会比低频的要大，同时，PR/MPDI二次多项式大气透过率所导致的发射率的误差估算为0.02/0.003(6.9GHz)，0.04/0.007(10.7GHz)和 0.1/0.01(18.7GHz)。6.9GHz、10.7GHz和18.7GHz的水平极化发射率 MPDI的二次多项式拟合方法具有很好的精度。
通过对上述方法的权衡，又由于6.9GHz对于地表土壤水分的敏感性很强，因此选择最小的RMSE的10.7GHz和18.7GHz的MPDI二次多项式的发射率估算(在18.7Hz，使用MPDI的二次多项式水平极化的发射率估算效果比较好)，估算的发射率和观测计算的发射率如图5所示。通过RMSE和相关系数，发现在两个测试区域的水平发射率的预测在 18.7GHz上最好。
从图5可以看到，中国中部地区作物的MPDI显示了其很好的预测性，基本上其RMSE均小于0.009，这导致亮温的变化不会超过2.7K(假设地表为300K)。而在青藏高原草地上，其RMSE超过了0.01，这对于亮温在草地上的模拟计算会带来较大的误差，但依然在0.02的精度范围。
分别对全球的14种主要地物的发射率的估算能力进行了分析，除了雪冰、裸地和半沙漠地区的水平发射率预测误差非常之大，其余的均在允许误差范围之列。10.7GHz全球范围内不同地表发射率的预测能力在森林覆盖区域很可靠，其RMSE<0.01，而在稀疏灌丛，稀树大草原、作物和城市和建成区。RMSE显示其全球的发射率估算方程并不是很好。
实施例二，通过采用上述方法的计算全球的地表发射率在36.5GHz 垂直和极化的发射率的结果。该结果显示了预想的空间结构，在V极化大于H极化，这符合模型的描述，冰雪覆盖地区的极化差大，如南极大陆、格陵兰岛冰雪表面，这是由于在冰雪表面的散射作用的影响。在沙漠和半干旱地区也显示了同样的现象。青藏高原在两种极化都出现了低发射率的情况，同样的现象也出现在沙漠和半干旱地带，如撒哈拉沙漠、澳大利亚中部等沙漠地区。
基于上述分析，可以看出瞬时和半月平均的基于MODIS的地表发射率都没有很好的全球覆盖，这些不能取值的点均为云或者雨的影响，可以看出这些区域在长达半个月的时间长期被云层覆盖，因此这就是需要进行大气校正工作的直接缘由之一。
从表2可以看出，相比于平均发射率，瞬时发射率值的动态范围相对来说比较小，在每个波段的最小瞬时发射率都大于最小的平均发射率，而最大的瞬时发射率却小于最大的平均发射率。半个月的平均发射率比12-08-2006的瞬时发射率要大0.0015左右，这个变化对于发射率的总体来说，其变化性比较小。
表2地表发射率统计值
Figure BDA0001108842260000151
Figure BDA0001108842260000161
更加精确的统计结果如图6和图7所示，从统计图看以看出，大部分的陆面发射率都在0.9左右，垂直极化大于水平极化的值，很小一部分值出现了异常现象(>1.0的情况)。随着频率的增加，异常值的比例，在V极化从2.6％降到了0.39％，而H极化从0.14％降到了0.09％。相比于发射率的情况，平均发射率的结果也出现了相同的趋势，对于V 极化，从3.82％降到了0.73％，而H极化从0.13％降到了0.09％。这种情形可以解释为具有长波长的的低频具备更深的穿透能力，并且对于地表土壤水分具有更强的敏感性。
使用上述测试的MPDI的二次多项式发射率估算方法，在全球范围内按照不同的分区(按照地理位置)进行全球瞬时地表发射率云层覆盖条件下的预测估算。首先我们把全球按照地理位置分为亚洲(南部)、欧亚大陆(北部)、北美、南美和非洲，然后再通过晴空条件确定MPDI二次多项式来分别确定方程。
对于瞬时发射率的有效性进行验证性统计分析，分别对预测区和晴空区和14天平均值进行统计对比，其结果如图8所示，其显示了 10.7GHz和18.7GHz情况下，预测发射率和瞬时发射率相对于平均发射率差的统计分布图。
如图8的统计分布直方图，可以看出，相比于瞬时发射率和平均发射率差值统计直方图来说，二次多项式MPDI方法的预测结果能提供合理的结果，但是也存在着一定的差别，MPDI的二次多项式估算发射率略微要大于平均发射率，而基于MODIS的地表发射率相比情况下，比其他发射率来说略小于平均值。
使用MODIS的地表发射率作为参考，当对比瞬时发射率和平均发射率的时候，发现95％的离差都小于0.05，该值是一天地表发射率的变化值，MPDI二次多项式估算的发射率显示了更大的离差，在 10.7GHz的时候约为0.075，而在18.7GHz的时候为0.05，而MODIS 晴空条件下的发射率和平均发射率的差值在10.7GHz的时候为0.05，而在18.7GHz的时候小于0.04。
对18.7GHz来说，由于其对地表土壤水分的变化敏感性不如 10.7GHz，因此在计算MPDI二次多项式以及MODIS的瞬时发射率相对于平均发射率的时候，发现其相差的范围会较10.7GHz小。获得的地表发射率还需要通过其他方法进行验证，比如地表温度、降水或云液态水的反演，间接的验证该发射率估算的有效性。
以上实施例仅用于说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (8)

1.一种被动微波遥感瞬时地表发射率估算方法，其特征在于，所述方法包括如下步骤：
使用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据，根据所述地表数据计算晴空条件下的瞬时地表发射率；
获取所述晴空条件下的不同频率水平和垂直极化瞬时地表发射率；
根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR；
根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程；
采用所述二次多项式和线性回归方程计算云层覆盖条件下的地表发射率的水平发射率；
其中，所述根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程的步骤包括如下步骤：
eh＝B0+B1·PR/MPDI+B2·PR/MPDI2
eh1＝B0+B1·PR/MPDI6.9+B2·PR/MPDI10.7+B3·PR/MPDI18.7
其中，eh为晴空条件下的地表发射率的水平极化发射率，eh1云层覆盖条件下的地表发射率水平发射率，MPDI为微波极化差指数，MPDI6.9,MPDI10.7,MPDI18.7分别为在频率为C波段，X波段，Ku波段下的微波极化差指数，B0,B1,B2,B3分别为晴空条件下的局部地区地表参数，其中晴空条件下的局部地区地表参数与云层覆盖条件下的该局部地区地表参数相同。
2.根据权利要求1所述的方法，其特征在于，根据所述地表数据计算晴空条件下的瞬时地表发射率的步骤包括根据如下公式计算获取晴空条件下的瞬时地表发射率的步骤：
Figure FDA0002081884690000021
其中，Tbp为卫星观测亮温，v为频率，θ为辐射计观测角度，Ts为地表温度，Γ(v,θ)为大气的吸收因子，Tatm↓(v,θ)和Tatm↑(v,θ)为向下或向上的大气辐射，TCB为宇宙背景辐射。
3.根据权利要求2所述的方法，其特征在于，所述根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR的步骤包括根据如下公式计算微波极化差指数MPDI和极化率PR的步骤：
PR＝Tbh/Tbv
MPDI＝(Tbv-Tbh)/(Tbv+Tbh)
其中，Tbh和Tbv为相同频率下的水平极化和垂直极化的亮温。
4.根据权利要求1～3中任一项所述的方法，其特征在于，所述亮温Tbp采用AMSR-E L2A双极化升/降轨道的亮度温度数据，所述亮度温度数据采用算术平均方法转换到同一投影格式；
大气参数采用MODIS/Aqua 5-Min Level2的20层大气廓线数据，所述大气廓线数据采用算术平均方法转换到同一投影格式；
地表温度数据采用MODIS/Aqua Land Surface Temperature/Emissivity 5-Min L2Swath-1km测量的数据，所述测量的数据采用算术平均方法转换到同一投影格式，测量数据中的质量控制数据QA作为晴空条件的判别原则。
5.一种被动微波遥感瞬时地表发射率估算的装置，其特征在于，所述装置包括：
用于使用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据、根据所述地表数据计算晴空条件下的瞬时地表发射率的装置；
用于获取所述晴空条件下的不同频率水平和垂直极化瞬时地表发射率装置；
用于根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR的装置；
用于根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程的装置；
用于采用二次多项式和线性回归方程计算云层覆盖条件下的地表发射率的水平发射率的装置；
其中，用于根据晴空条件下的瞬时地表发射率的水平极化发射率、以及所述微波极化差指数MPDI和极化率PR建立二次多项式和线性回归方程的装置包括根据如下公式计算云层覆盖条件下的地表发射率水平发射率方程的装置：
eh＝B0+B1·PR/MPDI+B2·PR/MPDI2
eh1＝B0+B1·PR/MPDI6.9+B2·PR/MPDI10.7+B3·PR/MPDI18.7
其中，eh为晴空条件下的地表发射率的水平极化发射率，eh1云层覆盖条件下的地表发射率水平发射率，MPDI为微波极化差指数，MPDI6.9,MPDI10.7,MPDI18.7分别为在频率为C波段，X波段，Ku波段下的微波极化差指数，B0,B1,B2,B3分别为晴空条件下的局部地区地表参数，其中晴空条件下的局部地区地表参数与云层覆盖条件下的该局部地区地表参数相同。
6.根据权利要求5所述的装置，其特征在于，用于使用星载被动微波传感器和中分辨率成像光谱仪同步观测地表数据、根据所述地表数据计算晴空条件下的瞬时地表发射率的装置包括根据如下公式计算获取晴空条件下的瞬时地表发射率的装置：
Figure FDA0002081884690000031
其中，Tbp为卫星观测亮温，v为频率，θ为辐射计观测角度，Ts为地表温度，Γ(v,θ)为大气的吸收因子，Tatm↓(v,θ)和Tatm↑(v,θ)为向下或向上的大气辐射，TCB为宇宙背景辐射。
7.根据权利要求6所述的装置，其特征在于，用于根据相同频率下水平极化的亮温Tbh和垂直极化的亮温Tbv确定微波极化差指数MPDI和极化率PR的装置包括根据如下公式计算微波极化差指数MPDI和极化率PR的装置：
PR＝Tbh/Tbv
MPDI＝(Tbv-Tbh)/(Tbv+Tbh)
其中，Tbh和Tbv为相同频率下的水平极化和垂直极化的亮温。
8.根据权利要求5～7中任一项所述的装置，其特征在于，所述亮温Tbp采用装置AMSR-EL2A双极化升/降轨道的亮度温度数据；
大气参数采用装置MODIS/Aqua 5-Min Level2的20层大气廓线数据；
地表温度数据采用装置MODIS/Aqua Land Surface Temperature/Emissivity 5-MinL2 Swath-1km测量的数据。
CN201610800047.5A 2016-08-31 2016-08-31 一种被动微波遥感瞬时地表发射率估算方法及装置 Active CN106372434B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610800047.5A CN106372434B (zh) 2016-08-31 2016-08-31 一种被动微波遥感瞬时地表发射率估算方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610800047.5A CN106372434B (zh) 2016-08-31 2016-08-31 一种被动微波遥感瞬时地表发射率估算方法及装置

Publications (2)

Publication Number Publication Date
CN106372434A CN106372434A (zh) 2017-02-01
CN106372434B true CN106372434B (zh) 2020-03-06

Family

ID=57900269

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610800047.5A Active CN106372434B (zh) 2016-08-31 2016-08-31 一种被动微波遥感瞬时地表发射率估算方法及装置

Country Status (1)

Country Link
CN (1) CN106372434B (zh)

Families Citing this family (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107451348B (zh) * 2017-07-25 2020-06-09 中国科学院东北地理与农业生态研究所 一种l波段森林微波辐射传输特性的优化估算方法
CN108827878B (zh) * 2018-04-08 2021-05-18 中国科学院国家空间科学中心 一种地表气压的被动微波遥感探测方法
CN108982548B (zh) * 2018-07-20 2020-01-17 浙江大学 一种基于被动微波遥感数据的地表土壤水分反演方法
US10996179B2 (en) 2019-03-11 2021-05-04 Skaha Remote Sensing Ltd. System and method to detect ground moisture
CN110865391B (zh) * 2019-11-14 2021-09-21 清华大学 用于目标增强的毫米波太赫兹多极化成像方法及成像系统
CN110879923B (zh) * 2019-12-04 2023-08-11 北京中科技达科技有限公司 阴天条件下长波下行辐射估算方法、存储介质及电子设备
CN111337531B (zh) * 2020-02-26 2021-04-27 中国科学院遥感与数字地球研究所 土壤冻融状态确定方法、装置和电子设备
CN111983333B (zh) * 2020-08-20 2021-05-25 自然资源部国土卫星遥感应用中心 一种分离被动微波遥感混合像元组分土壤辐射信号的方法
CN112395989B (zh) * 2020-11-18 2023-07-14 北京师范大学 一种用于多卫星传感器的积雪覆盖度混合像元分解方法
CN113063806A (zh) * 2021-03-22 2021-07-02 自然资源部国土卫星遥感应用中心 一种多通道微波辐射数据反演裸露土壤水分剖面的方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102411143A (zh) * 2010-09-20 2012-04-11 中国科学院遥感应用研究所 分离多角度卫星观测信号中的植被信息和地表发射信号的方法
CN103149220A (zh) * 2013-01-30 2013-06-12 中国科学院对地观测与数字地球科学中心 一种单频率微波辐射计土壤水分的反演方法
CN103424405A (zh) * 2013-08-27 2013-12-04 华南农业大学 一种基于hj-1a/1b ccd数据的干旱监测方法
CN103969268A (zh) * 2014-04-04 2014-08-06 中国科学院华南植物园 地表土壤物理参数的被动微波遥感反演方法
CN105466957A (zh) * 2016-01-07 2016-04-06 中国农业科学院农业资源与农业区划研究所 一种fy-3c被动微波数据估算土壤湿度的方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102411143A (zh) * 2010-09-20 2012-04-11 中国科学院遥感应用研究所 分离多角度卫星观测信号中的植被信息和地表发射信号的方法
CN103149220A (zh) * 2013-01-30 2013-06-12 中国科学院对地观测与数字地球科学中心 一种单频率微波辐射计土壤水分的反演方法
CN103424405A (zh) * 2013-08-27 2013-12-04 华南农业大学 一种基于hj-1a/1b ccd数据的干旱监测方法
CN103969268A (zh) * 2014-04-04 2014-08-06 中国科学院华南植物园 地表土壤物理参数的被动微波遥感反演方法
CN105466957A (zh) * 2016-01-07 2016-04-06 中国农业科学院农业资源与农业区划研究所 一种fy-3c被动微波数据估算土壤湿度的方法

Non-Patent Citations (9)

* Cited by examiner, † Cited by third party
Title
A Practical Method for Retrieving Land Surface Temperature From AMSR-E Over the Amazon Forest;Huilin Gao等;《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》;20080131;第46卷(第1期);正文第II部分 *
An examination of methods for estimating land surface microwave emissivity;Yudong Tian等;《Journal of Geophysical Research: Atmospheres》;20151110;第120卷(第21期);正文第3.2节 *
THE AMSR-E INSTANTANEOUS EMISSIVITY ESTIMATION AND ITS CORRELATION, FREQUENCY DEPENDENCY ANALYSIS OVER DIFFERENT LAND COVERS;Yubao Qiu等;《Geoscience and Remote Sensing Symposium》;20080711;第749-752页 *
Using a Microwave Emission Model to Estimate Soil Moisture from ESTAR Observations during SGP99;HUILIN GAO等;《JOURNAL OF HYDROMETEOROLOGY》;20040229;第5卷(第1期);第49-63页 *
Using the Special Sensor Microwave/Imager to Monitor Land Surface Temperatures, Wetness, and Snow Cover;ALAN BASIST等;《JOURNAL OF APPLIED METEOROLOGY》;19980930;第37卷(第9期);第888-911页 *
基于AMSR-E 的微波波段地表发射率反演——以青藏高原为例;王永前等;《中国科学:地球科学》;20130220;第43卷(第2期);第271-279页 *
基于AMSR-E的全球陆表被动微波发射率数据集;邱玉宝等;《遥感技术与应用》;20160820;第31卷(第4期);正文第1-2节、第4.4节 *
基于被动微波的寒旱区地表温度反演;彭丹青等;《冰川冻土》;20090415;第31卷(第2期);第233-238页 *
被动微波遥感反演地表发射率研究进展;吴莹等;《国土资源遥感》;20121215(第4期);第1-7页 *

Also Published As

Publication number Publication date
CN106372434A (zh) 2017-02-01

Similar Documents

Publication Publication Date Title
CN106372434B (zh) 一种被动微波遥感瞬时地表发射率估算方法及装置
Ukhurebor et al. Review of methodology to obtain parameters for radio wave propagation at low altitudes from meteorological data: new results for Auchi area in Edo State, Nigeria
He et al. Validation of MODIS derived aerosol optical depth over the Yangtze River Delta in China
Liu et al. A multi-sensor study of water vapour from radiosonde, MODIS and AERONET: a case study of Hong Kong.
Xian et al. Fengyun meteorological satellite products for earth system science applications
Sayer et al. Global and regional evaluation of over-land spectral aerosol optical depth retrievals from SeaWiFS
Sorek-Hamer et al. Assessment of PM2. 5 concentrations over bright surfaces using MODIS satellite observations
Seo et al. Estimation of PM 10 concentrations over Seoul using multiple empirical models with AERONET and MODIS data collected during the DRAGON-Asia campaign
CN102539336B (zh) 基于环境一号卫星的可吸入颗粒物估算方法及系统
Adediji et al. Distribution of radio refractivity gradient and effective earth radius factor (k-factor) over Akure, South Western Nigeria
Mengistu Tsidu et al. Observations of precipitable water vapour over complex topography of Ethiopia from ground-based GPS, FTIR, radiosonde and ERA-Interim reanalysis
Yu et al. Supplement of the radiance-based method to validate satellite-derived land surface temperature products over heterogeneous land surfaces
Huang et al. High resolution surface radiation products for studies of regional energy, hydrologic and ecological processes over Heihe river basin, northwest China
Feng et al. Merging ground-based sunshine duration observations with satellite cloud and aerosol retrievals to produce high-resolution long-term surface solar radiation over China
Guzinski et al. Remotely sensed land-surface energy fluxes at sub-field scale in heterogeneous agricultural landscape and coniferous plantation
Merrikhpour et al. Improving the algorithm of extracting regional total precipitable water vapor over land from MODIS images
Zibordi et al. An evaluation of marine regions relevant for ocean color system vicarious calibration
De Abreu et al. Comparison of in situ and AVHRR-derived broadband albedo over Arctic sea ice
Wang et al. Recent declines in global water vapor from MODIS products: Artifact or real trend?
Ruiz-Arias et al. Bias induced by the AOD representation time scale in long-term solar radiation calculations. Part 1: Sensitivity of the AOD distribution to the representation time scale
He et al. Direct estimation of land surface albedo from simultaneous MISR data
Wang et al. Analysis of photosynthetically active radiation in Northwest China from observation and estimation
Damiani et al. Satellite-derived UV irradiance for a region with complex morphology and meteorology: comparison against ground measurements in Santiago de Chile
Borbas et al. MODIS atmospheric profile retrieval algorithm theoretical basis document
CN110658325A (zh) 一种确定sm监测中lst/fvc空间的方法

Legal Events

Date Code Title Description
C06 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