CN115524763B - 一种多时相高分辨率山地卫星影像地形辐射校正方法 - Google Patents

一种多时相高分辨率山地卫星影像地形辐射校正方法 Download PDF

Info

Publication number
CN115524763B
CN115524763B CN202211181020.4A CN202211181020A CN115524763B CN 115524763 B CN115524763 B CN 115524763B CN 202211181020 A CN202211181020 A CN 202211181020A CN 115524763 B CN115524763 B CN 115524763B
Authority
CN
China
Prior art keywords
angle
satellite
radiation
reflectivity
coefficient
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
CN202211181020.4A
Other languages
English (en)
Other versions
CN115524763A (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 Mountain Hazards and Environment IMHE of CAS
Original Assignee
Institute of Mountain Hazards and Environment IMHE 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 Mountain Hazards and Environment IMHE of CAS filed Critical Institute of Mountain Hazards and Environment IMHE of CAS
Priority to CN202211181020.4A priority Critical patent/CN115524763B/zh
Publication of CN115524763A publication Critical patent/CN115524763A/zh
Application granted granted Critical
Publication of CN115524763B publication Critical patent/CN115524763B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V13/00Manufacturing, calibrating, cleaning, or repairing instruments or devices covered by groups G01V1/00 – G01V11/00
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01NINVESTIGATING OR ANALYSING MATERIALS BY DETERMINING THEIR CHEMICAL OR PHYSICAL PROPERTIES
    • G01N21/00Investigating or analysing materials by the use of optical means, i.e. using sub-millimetre waves, infrared, visible or ultraviolet light
    • G01N21/17Systems in which incident light is modified in accordance with the properties of the material investigated
    • G01N21/25Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands
    • G01N21/27Colour; Spectral properties, i.e. comparison of effect of material on the light at two or more different wavelengths or wavelength bands using photo-electric detection ; circuits for computing concentration
    • G01N21/274Calibration, base line adjustment, drift correction
    • YGENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02ATECHNOLOGIES FOR ADAPTATION TO CLIMATE CHANGE
    • Y02A90/00Technologies having an indirect contribution to adaptation to climate change
    • Y02A90/10Information and communication technologies [ICT] supporting adaptation to climate change, e.g. for weather forecasting or climate simulation

Landscapes

  • Physics & Mathematics (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Physics & Mathematics (AREA)
  • Analytical Chemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Health & Medical Sciences (AREA)
  • Pathology (AREA)
  • Chemical & Material Sciences (AREA)
  • Immunology (AREA)
  • Biochemistry (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Engineering & Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Manufacturing & Machinery (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Image Processing (AREA)

Abstract

本发明提供了一种多时相高分辨率山地卫星影像地形辐射校正方法,所述方法包括,获取第一核驱动模型和多时相高分辨率山地卫星图像,所述第一核驱动模型用于表征二向反射率与观测角度之间的关系;基于所述多时相高分辨率山地卫星图像,对所述第一核驱动模型进行反演得到第二核驱动模型,所述第二核驱动模型用于表征倾斜地表任意观测角度下的二向反射率;采用所述第二核驱动模型对所述多时相高分辨率山地卫星图像进行地形辐射校正;具有降低多时相高分辨率山地卫星图像中地形影响因素的效果。

Description

一种多时相高分辨率山地卫星影像地形辐射校正方法
技术领域
本发明涉及多光谱遥感卫星影像的计算机处理领域,具体而言,涉及一种多时相高分辨率山地卫星影像地形辐射校正方法。
背景技术
对光学遥感而言,遥感成像过程中除了受传感器自身电子元器件老化造成的辐射性能衰减外,传感器记录的电磁波信号由于两次经过大气,因此,消除大气对光学遥感影像的影响以获取地表真实的反射率是遥感定量化应用的重要步骤。此外,除了受到大气影响,电磁波到达地面后还受到地形影响,山地光谱辐射畸变十分复杂。在山地光学遥感影像中,地形起伏导致不同朝向的坡面接收的太阳辐射不同,致使相同的地物在不同坡度、坡向条件有较大的光谱差异。这种地形效应将会显著影响地表参量提取精度以及定量遥感产品在山区的应用;地形对山地地表反射率的影响表现在改变目标像元的局地成像几何条件及入射辐射能量两个方面,导致目标像元的成像几何条件发生改变。而同一地表覆被类型条件下,不同坡度和坡向下像元的局地观测角度以及观测得到的反射率明显不同;地形不仅会改变成像几何,还会对目标像元入射能量产生显著影响。
对陆地遥感而言,传感器记录的电磁信号需要消除上述各种因素的影响,以恢复为真实的地表反射率。随着遥感在变化监测等领域应用的深入,以及不同传感器影像数据的对比分析需求,由太阳与观测角度变化所引起的地物反射率变化必须予以考虑并消除。传统高分辨率卫星影像由于回归周期长,且受云雾干扰,难以构成多角度观测,进而难以消除观测几何造成的二向反射差异。
发明内容
本发明的目的在于提供一种多时相高分辨率山地卫星影像地形辐射校正方法,其目的在于消除多时相密集高分辨率山地卫星影像中的地形影响因素。
本发明的实施通过以下技术方案实现:
一种多时相高分辨率山地卫星影像地形辐射校正方法,所述方法包括,
获取第一核驱动模型和多时相高分辨率山地卫星图像,所述第一核驱动模型用于表征二向反射率与观测角度之间的关系;
基于所述多时相高分辨率山地卫星图像,对所述第一核驱动模型进行反演得到第二核驱动模型,所述第二核驱动模型用于表征倾斜地表任意观测角度下的二向反射率;
采用所述第二核驱动模型对所述多时相高分辨率山地卫星图像进行地形辐射校正。
可选地,所述基于所述多时相高分辨率山地卫星图像,对所述第一核驱动模型进行反演得到第二核驱动模型,包括:
获取水平地表的卫星观测角和太阳观测角,所述卫星观测角包括卫星观测天顶角和卫星观测方位角,所述太阳观测角包括太阳入射天顶角和太阳入射方位角;
基于所述第一核驱动模型,将水平地表与倾斜地表的角度关系进行归一化;
获取所述第一核驱动模型中水平地表的二向反射率;
基于所述归一化的角度关系和所述水平地表的二向反射率,获取所述倾斜地表的二向反射率;
获取所述水平地表的大气辐射传输模型,所述大气辐射传输模型包括水平地表朗伯体假设情况的大气上界反射率、水平地表非朗伯体假设情况的大气上界反射率和水平地表非朗伯体假设情况的大气上界辐射亮度;
基于所述水平地表的大气辐射传输模型和所述倾斜地表的二向反射率,获取所述倾斜地表的地表反射率;
基于所述倾斜地表的地表反射率、所述第一核驱动模型和所述倾斜地表的二向反射率,得到所述第二核驱动模型中各向同性系数、体散射核系数和几何散射核系数的关系表达式;
基于所述第二核驱动模型中各向同性系数、体散射核系数和几何散射核系数的关系表达式和所述多时相高分辨率山地卫星图像,得到所述第二核驱动模型。
可选地,所述基于所述第一核驱动模型,将所述水平地表与倾斜地表的角度关系进行归一化,包括:
基于所述第一核驱动模型,获取水平地表条件下的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角的角度关系表达式;
获取倾斜地表的坡角数据和坡向数据;
基于所述水平地表条件下的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角的角度关系表达式及所述倾斜地表的坡角数据和坡向数据,得到所述倾斜地表条件下太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角的归一化角度关系表达式。
可选地,所述基于所述水平地表的大气辐射传输模型和所述倾斜地表的二向反射率,获取所述倾斜地表的地表反射率,包括:
基于所述水平地表朗伯体假设情况的大气上界反射率、水平地表非朗伯体假设情况的大气上界反射率和所述倾斜地表的二向反射率,获取所述倾斜地表的大气上界辐射表达式;
基于所述水平地表非朗伯体假设情况的大气上界辐射亮度,获取所述大气上界辐射亮度经倾斜地表重新分配后,所述倾斜地表的太阳辐射与所述水平地表的太阳辐射之间的相互关系式;所述倾斜地表的大气上界辐射亮度包括太阳直射辐射、天空散射辐射和临近地表反射辐射;
基于所述相互关系式,获取所述倾斜地表的大气上界辐射亮度中太阳直射辐射、天空散射辐射和临近地表反射辐射与所述倾斜地表受到的太阳总辐射之间的比例关系式;
基于所述倾斜地表的大气上界辐射表达式、相互关系式和所述比例关系式,获取所述倾斜地表的地表反射率。
可选地,所述基于所述相互关系式,获取所述倾斜地表的大气上界辐射亮度中太阳直射辐射、天空散射辐射和临近地表反射辐射与所述倾斜地表受到的太阳总辐射之间的比例关系式,包括:
基于所述归一化角度关系表达式,获取所述太阳辐射的地形可见因子系数;
获取所述倾斜地表的邻近地表反射率;
获取所述太阳辐射的大气透过率系数;
基于所述地形可见因子系数、邻近地表反射率、大气透过率系数和所述相互关系式,获取所述比例关系式。
可选地,所述基于所述倾斜地表的大气上界辐射表达式、相互关系式和所述比例关系式,获取所述倾斜地表的地表反射率,包括:
联立所述倾斜地表的大气上界辐射表达式、相互关系式和所述比例关系式,获取所述倾斜地表的地表反射率表达式。
可选地,所述基于所述倾斜地表的地表反射率、所述第一核驱动模型和所述倾斜地表的二向反射率,得到所述第二核驱动模型中各向同性系数、体散射核系数和几何散射核系数的关系表达式,包括:
将所述倾斜地表的二向反射率带入所述倾斜地表的地表反射率表达式中,并获取所述地表反射率的一阶近似;
将所述倾斜地表的坡角数据、坡向数据和地表反射率的一阶近似带入所述地表反射率表达式中,获取所述各向同性系数、体散射核系数和几何散射核系数的关系表达式。
可选地,所述基于所述第二核驱动模型中各向同性系数、体散射核系数和几何散射核系数的关系表达式和所述多时相高分辨率山地卫星图像,得到所述第二核驱动模型,包括:
获取所述多时相高分辨率山地卫星图像中的大气上界反射率、太阳入射天顶角、太阳入射方位角、卫星观测天顶角、卫星观测方位角、坡角和坡向;
将所述多时相高分辨率山地卫星图像中的大气上界反射率、太阳入射天顶角、太阳入射方位角、卫星观测天顶角、卫星观测方位角、坡角和坡向带入所述各向同性系数、体散射核系数和几何散射核系数的关系表达式,采用最小二乘法对所述各向同性系数、体散射核系数和几何散射核系数进行拟合,得到所述各向同性系数、体散射核系数和几何散射核系数;
将所述各向同性系数、体散射核系数和几何散射核系数带入所述第一核驱动模型,得到所述第二核驱动模型。
可选地,所述采用所述第二核驱动模型对所述多时相高分辨率山地卫星图像进行地形辐射校正,包括:
获取所述多时相高分辨率山地卫星图像中每个像元倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角;
基于所述第二核驱动模型和所述每个像元中倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角,获取每个像元所对应的地表二向反射率;
基于所述每个像元所对应的地表二向反射率,对所述多时相高分辨率山地卫星图像进行同一太阳观测角度和传感器观测角度的归一化协同校正。
可选地,所述获取所述多时相高分辨率山地卫星图像中每个像元倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角,包括:
逐个像元获取所述多时相高分辨率山地卫星图像中水平地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角;
获取所述多时相高分辨率山地卫星图像中每个所述像元的地理纬度、地理经度、赤纬、时角和卫星经度;
基于所述多时相高分辨率山地卫星图像中每个像元水平地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角及所述多时相高分辨率山地卫星图像中每个所述像元的地理纬度、地理经度、赤纬、时角和卫星经度,获取所述多时相高分辨率山地卫星图像中每个像元倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角。
本发明实施例的技术方案至少具有如下优点和有益效果:本发明基于核驱动模型基本理论,根据多时相山地卫星影像形成的多角度观测数据,可反演获得适于山地地表的第二核驱动模型的系数,该系数与地表的状态,如粗糙度、导电系数、植被覆盖等信息有关,与观测角度和地形无关,可根据第二核驱动模型可将多时相山地卫星影像观测角度归一化到统一观测几何条件;从而达到消除多时相密集高分辨率山地卫星影像中的地形影响因素的效果。
附图说明
图1为本发明其中一个实施例提供的一种多时相高分辨率山地卫星影像地形辐射校正方法的步骤流程示意图;
图2为本发明其中一个实施例中水平地表的角度示意图;
图3为本发明其中一个实施例中倾斜地表的角度示意图;
图4为本发明其中一个实施例对多时相高分辨率山地卫星图像校正前的示意图;
图5为本发明其中一个实施例对多时相高分辨率山地卫星图像校正后的示意图;
图6为本发明其中一个实施例提供的实验区地形条件及选择的HJ-1A/B影像;
图7为本发明其中一个实施例提供的实验区地形条件及选择的HJ-1A/B影像;
图8为本发明其中一个实施例提供的方法与MODIS的反演结果对比示意图;
图9为本发明其中一个实施例校正前二向反射率参数累积强度示意图;
图10为本发明其中一个实施例校正后二向反射率参数累积强度示意图。
具体实施方式
为使本发明实施例的目的、技术方案和优点更加清楚,下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例是本发明一部分实施例,而不是全部的实施例。通常在此处附图中描述和示出的本发明实施例的组件可以以各种不同的配置来布置和设计。
实施例1
本申请实施例提供一种多时相高分辨率山地卫星影像地形辐射校正方法,参照图1,所述方法包括以下步骤,
S1,获取第一核驱动模型和多时相高分辨率山地卫星图像,所述第一核驱动模型用于表征二向反射率与观测角度之间的关系;
多时相高分辨率山地卫星影像指通过遥感卫星在短时间内,高频率获取的高分辨率的影像。
在同一地表覆被类型条件下,不同坡度和坡向下像元的局地观测角度以及观测得到的反射率明显不同,地形不仅会改变成像几何,还会对目标像元入射能量产生显著影响。
本实施例中,采用核驱动模型对卫星图像中的反射率进行校正,减少卫星图像中反射率带来的影响,从而得到更加清晰的卫星图像。
S2,基于所述多时相高分辨率山地卫星图像,对所述第一核驱动模型进行反演得到第二核驱动模型,所述第二核驱动模型用于表征倾斜地表任意观测角度下的二向反射率;
第一核驱动模型是传统的模型,仅能够针对地形为水平地表的像元进行校正,若是采用第一核驱动模型校正多时相高分辨率山地卫星图像,其得到的结果会受到地形和观测角度的影响;通过对第一核驱动模型进行反演,得到的第二核驱动模型消除了地形的影响,得到的二向反射率能够更准确的对多时相高分辨率山地卫星图像进行校正,得到地表反射率更加一致的卫星图像。
其反演过程包括下述子步骤:
S11,获取水平地表的卫星观测角和太阳观测角,所述卫星观测角包括卫星观测天顶角和卫星观测方位角,所述太阳观测角包括太阳入射天顶角和太阳入射方位角;
如图2和图3所示,图2中示出了水平地表的角度关系,图3示出了山地倾斜地表的角度关系;在水平地表条件下,卫星观测天顶角、卫星观测方位角、太阳入射天顶角和太阳入射方位角可根据垂直于水平地表的法线进行确定;而在山地倾斜坡面条件下,地形起伏会导致局地太阳和卫星传感器角度关系发生变化,因为此时法线垂直于坡面而非水平面。
S12,基于所述第一核驱动模型,将所述水平地表与倾斜地表的角度关系进行归一化;
在本实施例中,归一化包括以下步骤:
基于所述第一核驱动模型,获取水平地表条件下的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角的角度关系表达式;
获取倾斜地表的坡角数据和坡向数据;
倾斜地表的坡角数据和坡向数据采用DEM衍生;数字高程模型(DigitalElevation Model),简称DEM,是通过有限的地形高程数据实现对地面地形的数字化模拟,它是用一组有序数值阵列形式表示地面高程的一种实体地面模型,是数字地形模型(Digital Terrain Model,简称DTM)的一个分支,其它各种地形特征值均可由此派生。
基于所述水平地表条件下的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角的角度关系表达式及所述倾斜地表的坡角数据和坡向数据,得到所述倾斜地表条件下太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角的归一化角度关系表达式。
在本实施例中,得到的倾斜地表条件下太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角的归一化角度关系表达式如下所示:
cos(it)=cos(θs)cos(θt)+sin(θs)sin(θt)cos(φst);
cos(et)=cos(θv)cos(θt)+si(θv)sin(θt)cos(φvt);
式中,θs为水平地表条件下的太阳入射天顶角;φs为水平地表条件下的太阳入射方位角;θv为水平地表条件下的卫星观测天顶角;φv为水平地表条件下的卫星观测方位角;it为倾斜地表的太阳入射天顶角;et为倾斜地表条件下的卫星观测天顶角;φt为相对方位角,相对方位角为倾斜地表条件下的太阳入射方位角φi与倾斜地表条件下的卫星观测方位角φe之间的差值。
S13,获取所述第一核驱动模型中水平地表的二向反射率;
S14,基于所述归一化的角度关系和所述水平地表的二向反射率,获取所述倾斜地表的二向反射率;
在本实施例中,第一核驱动模型的二向反射率表达式如下所示:
其中,为第一核驱动模型的体散射核;/>为几何光学散射核;Fiso为各向同性系数;Fgeo为几何散射核系数;Fvol为体散射核系数;/>为水平地表条件下太阳与卫星之间的相对方位角。
通过归一化的角度关系和所述水平地表的二向反射率得到的倾斜地表条件下的二向反射率如下所示:
ρs(it,et,δφt)=Fiso+FvolKvol(it,et,δφt)+FgeoKgeo(it,et,δφt);
式中,ρs为倾斜地表条件下的二向地表反射率,为倾斜地表条件下的半球方向反射率;/>为倾斜地表条件下的双半球反射率。
其中,
式中,gk1(et)、hk2(et)、Hk1和Hk2为相关关系式;
其中,方向半球反射率与半球方向反射率的计算相同。且半球方向反射率仅与出射方向有关,而方向半球反射率仅与入射方向有关。
S15,获取所述水平地表的大气辐射传输模型,所述大气辐射传输模型包括水平地表朗伯体假设情况的大气上界反射率、水平地表非朗伯体假设情况的大气上界反射率和水平地表非朗伯体假设情况的大气上界辐射亮度;
其中,水平地表朗伯体假设情况的大气上界反射率如下所示:
式中,ρTOA为水平地表朗伯体条件下的大气上界反射率;ρ0为大气程辐射反射率;Ts为水平地表朗伯体条件下的太阳方向的总透过率;Tv为水平地表朗伯体条件下的卫星方向的总透过率;S为大气下界反照率。
ρm为地表反射率。
水平地表非朗伯体假设情况的大气上界反射率如下所示:
式中,ts为水平地表非朗伯体条件下的太阳方向散射透过率;tv为水平地表非朗伯体条件下的卫星方向散射透过率。
S16,基于所述水平地表的大气辐射传输模型和所述倾斜地表的二向反射率,获取所述倾斜地表的地表反射率。
包括,基于所述水平地表朗伯体假设情况的大气上界反射率、水平地表非朗伯体假设情况的大气上界反射率和所述倾斜地表的二向反射率,获取所述倾斜地表的大气上界辐射表达式;
基于所述水平地表非朗伯体假设情况的大气上界辐射亮度,获取所述大气上界辐射亮度经倾斜地表重新分配后,所述倾斜地表的太阳辐射与所述水平地表的太阳辐射之间的相互关系式;所述倾斜地表的大气上界辐射亮度包括太阳直射辐射、天空散射辐射和临近地表反射辐射;
式中,是水平地表接获得的太阳直接辐射;/>是水平地表接获得的太阳散射辐射;ρadj为临近地表的平均反射率;Vt为地形可见因子系数。
基于所述相互关系式,获取所述倾斜地表的大气上界辐射亮度中太阳直射辐射、天空散射辐射和临近地表反射辐射与所述倾斜地表受到的太阳总辐射之间的比例关系式;
包括:
基于所述归一化角度关系表达式,获取所述太阳辐射的地形可见因子系数;
太阳辐射的地形可见因子系数通过下式获取:
Vt=1-Vd
其中,
获取所述倾斜地表的邻近地表反射率;倾斜地表的临近地表反射率通过多时相高分辨率山地卫星图像分析获取。
获取所述太阳辐射的大气透过率系数;通过下式获取:
式中,k为大气透过率系数;为地表法线相对太阳光线接收到的大气上界辐射。为地表法线相对太阳光线接收到的直接辐射。
基于所述地形可见因子系数、邻近地表反射率、大气透过率系数和所述相互关系式,获取所述比例关系式。获取的比例关系式如下所示:
Rdir=Edir/Eh
Rdif=Edif/Eh
R=E/Eh
式中,Rdir、Rdif和R为比例系数;Edir为倾斜地表接收到的太阳直射辐射;Edif为倾斜地表接收到的太阳散射辐射;Eh为水平地表接收到的所有太阳辐射;E为倾斜地表接收到的所有太阳辐射。
基于所述倾斜地表的大气上界辐射表达式、相互关系式和所述比例关系式,获取所述倾斜地表的地表反射率;
联立所述倾斜地表的大气上界辐射表达式、相互关系式和所述比例关系式,获取所述倾斜地表的地表反射率表达式。
地表反射率表达式如下所示:
S17,基于所述倾斜地表的地表反射率、所述第一核驱动模型和所述倾斜地表的二向反射率,得到所述第二核驱动模型中各向同性系数、体散射核系数和几何散射核系数的关系表达式;
包括:
将所述倾斜地表的二向反射率带入所述倾斜地表的地表反射率表达式中,并获取所述地表反射率的一阶近似;
将朗伯体假设下的经大气校正后的地表反射率作为真实地表反射率的一阶近似。
将所述倾斜地表的坡角数据、坡向数据和地表反射率的一阶近似带入所述地表反射率表达式中,获取所述各向同性系数、体散射核系数和几何散射核系数的关系表达式。
将倾斜地表的坡角数据、坡向数据和地表反射率的一阶近似带入地表反射率表达式中得到如下:
将该式的两端同时除以(Rdir+Rdif);得到下式:
其中,
通过化简,可将该式认作一个以Fiso为常数项,以Fvol和Fgeo为变量的二元一次方程。
S18,基于所述第二核驱动模型中各向同性系数、体散射核系数和几何散射核系数的关系表达式和所述多时相高分辨率山地卫星图像,得到所述第二核驱动模型。
获取所述多时相高分辨率山地卫星图像中的大气上界反射率、太阳入射天顶角、太阳入射方位角、卫星观测天顶角、卫星观测方位角、坡角和坡向;
将所述多时相高分辨率山地卫星图像中的大气上界反射率、太阳入射天顶角、太阳入射方位角、卫星观测天顶角、卫星观测方位角、坡角和坡向带入所述各向同性系数、体散射核系数和几何散射核系数的关系表达式,采用最小二乘法对所述各向同性系数、体散射核系数和几何散射核系数进行拟合,得到所述各向同性系数、体散射核系数和几何散射核系数;
将所述各向同性系数、体散射核系数和几何散射核系数带入所述第一核驱动模型,得到所述第二核驱动模型。
通过构建上述关于各向同性系数、体散射核系数和几何散射核系数的二元一次方程,再利用多时相高分辨率山地卫星图像中坡面多个观测角度获得的多角度地表反射率数据,再结合地形坡度和坡向数据;根据该核驱动模型可以计算获得该式左侧的式子:
和该式右侧的式子:
的解;再采用最小二乘法拟合就能够得到各向同性系数、体散射核系数和几何散射核系数
S3,采用所述第二核驱动模型对所述多时相高分辨率山地卫星图像进行地形辐射校正。
包括以下子步骤:
获取所述多时相高分辨率山地卫星图像中每个像元倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角;包括:
逐个像元获取所述多时相高分辨率山地卫星图像中水平地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角;
获取所述多时相高分辨率山地卫星图像中每个所述像元的地理纬度、地理经度、赤纬、时角和卫星经度;
基于所述多时相高分辨率山地卫星图像中每个像元水平地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角及所述多时相高分辨率山地卫星图像中每个所述像元的地理纬度、地理经度、赤纬、时角和卫星经度,获取所述多时相高分辨率山地卫星图像中每个像元倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角。
在本实施例中,计算多时相高分辨率山地卫星图像中每个像元倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角采用下述公式:
sin hs=sinψ·sinδ+cosψ·cosδ·cosΩ;
/>
式中,ψ代表像元的地理纬度,λ代表像元的地理经度,δ为赤纬,Ω为时角,λs为卫星经度,π/2-hs为太阳天顶角。
基于所述第二核驱动模型和所述每个像元中倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角,获取每个像元所对应的地表二向反射率;
基于所述每个像元所对应的地表二向反射率,对所述多时相高分辨率山地卫星图像进行同一太阳观测角度和传感器观测角度的归一化协同校正。
参照图4和图5,图4和图5示出了通过本实施例提供的方法对多时相高分辨率山地卫星图像进行校正前后的对比。
通过将密集时相卫星观测获得的大气上界反射率、每个像元中倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角,获取每个像元所对应的地表二向反射率代入,地表各向同性系数、体散射核系数以及几何散射核系数表达式,就能够利用最小二乘法拟合获得山地地表各向同性系数、体散射核系数以及几何散射核系数,得到第二核驱动模型。
再根据第二核驱动模型和山地地表各向同性系数、体散射核系数以及几何散射核系数,计算获得卫星观测、太阳高度角统一的地表二向反射率,完成大气、地形和地表方向反射的协同校正;达到消除多时相密集高分辨率山地卫星影像中的地形影响因素的效果,得到的图像更加清晰。
本实施例基于核驱动模型基本理论,根据多时相山地卫星影像形成的多角度观测数据,可反演获得适于山地地表的第二核驱动模型的系数,该系数与地表的状态,如粗糙度、导电系数、植被覆盖等信息有关,与观测角度和地形无关,可根据第二核驱动模型可将多时相山地卫星影像观测角度归一化到统一观测几何条件;从而达到消除多时相密集高分辨率山地卫星影像中的地形影响因素的效果。
在其中一个实施例中,以具有高时空分辨率的中国国产环境减灾卫星(型号:HJ-1A/B)影像为例,开展高空间分辨率的二向反射率参数提取,进而实现基于二向反射率的多时相高分辨率山地卫星图像校正。
该实施例选择的数据为位于我国横断山区的环境减灾卫星HJ-1A/B CCD数据,影像空间分辨率30米。在我国环境减灾卫星星座(HJ-1A、1B)的轨道设计中,为实现48小时内对全球的重复覆盖功能,通过多次轨控及相位的初始化,现已形成了双星相位差180度2天重访的要求,并经过地面轨迹和星座的轨道维持策略。
为了获得研究区内尽可能多的有效观测数据,本实施例采用了一种基于规则格网的HJ-1A/B影像选择方法。其中格网的大小在此处设为2500×2500个像元,选择该大小的研究区主要分别考虑了上述环境减灾卫星轨道的漂移性和区域内尽可能多的观测数据。实验区地形条件及选择的HJ-1A/B影像如图6和图7所示。可以看出,整个试验区内植被类型复杂,主要包括针叶林、阔叶林、混交林及高山草甸等,植被覆盖度较高,16天内地表反射率变化不明显。
基于选择的试验区域,经过查询统计,最终选择2013年生长季8月4日至8月19日的HJ-1A/B CCD1/CCD2影像11景。理论上,当研究区内有效观测数据(无云干扰)大于3个时,可根据有效观测数据构建数量大于3的超定方程,利用最小二乘法可以实现核驱动模型中三系数的求解。
然而,在实际操作中,一般需要有效观测远远大于未知数个数,才能获得可靠的拟合结果。在MODIS的NBAR产品算法中,设定有效观测数据大于7次时(大于有效观测数据的1倍以上),模型能够实现全反演。
基于其中一个实施例中推导得出的第二核驱动模型,利用研究区11景HJ-1A/B影像,本实施例拟合获得了研究区地表第二核驱动模型对应的各向同性散射系数、几何散射核和体散射核系数的反演结果,并给出了同期MODIS MCD BRDF产品的参数对比。
从图8反演结果对比可以看出,本研究获得的二向反射率参数空间纹理细节更加明显,而MODIS的体散射核和几何散射核系数空间分辨率较低且较为凌乱。由此可见,本发明拟合获得的二向反射率参数更适用于高空间分辨率二向反射率的校正处理。
图8给出了采用其中一个实施例提供的方法针对研究区影像地形辐射校正前后的整体及局部对比效果。从图8中可以看出,总体上,原始影像地形校正前地形效应明显,阴坡区域地表反射率受地形遮蔽影响明显低于阳坡区域,而本发明方法地形辐射校正后地形效应均得到了显著的剔除。其中,从局部放大结果也可以看出,其中一个实施例采用方法由于采用30米空间分辨率的多时相影像直接拟合获得地表的同分辨率级别二向反射率参数,因此没有尺度偏差,能够较好的校正影像中的光谱变化较大的区域。
一般而言,由于地形效应,山地区域阴坡区域反射率较低,阳坡区域反射率较高,由图9可以看出,在由坡度和坡向组成的极坐标系下,反射率强度累积表现出显著的地形特征。而地形校正后,由图10图可以看出,反射率阴坡得到了很好的增强,反射率强度累积在阴坡和阳坡更加均一。山区原始影像中地表反射率与光照系数具有很强的相关性,即相同条件下的地表,光照系数越大,直射辐射越强,经地形辐射校正后这种相关性会减小。总体而言,本发明的地形辐射校正方法有较强的地形去相关效果,相关性得到了很好的去除。
以上仅为本发明的优选实施例而已,并不用于限制本发明,对于本领域的技术人员来说,本发明可以有各种更改和变化。凡在本发明的精神和原则之内,所作的任何修改、等同替换、改进等,均应包含在本发明的保护范围之内。

Claims (9)

1.一种多时相高分辨率山地卫星影像地形辐射校正方法,其特征在于,所述方法包括,
获取第一核驱动模型和多时相高分辨率山地卫星图像,所述第一核驱动模型用于表征二向反射率与观测角度之间的关系;
基于所述多时相高分辨率山地卫星图像,对所述第一核驱动模型进行反演得到第二核驱动模型,所述第二核驱动模型用于表征倾斜地表任意观测角度下的二向反射率;具体为,获取水平地表的卫星观测角和太阳观测角,所述卫星观测角包括卫星观测天顶角和卫星观测方位角,所述太阳观测角包括太阳入射天顶角和太阳入射方位角;
基于所述第一核驱动模型,将水平地表与倾斜地表的角度关系进行归一化;
获取所述第一核驱动模型中水平地表的二向反射率;
基于所述归一化的角度关系和所述水平地表的二向反射率,获取所述倾斜地表的二向反射率;
获取所述水平地表的大气辐射传输模型,所述大气辐射传输模型包括水平地表朗伯体假设情况的大气上界反射率、水平地表非朗伯体假设情况的大气上界反射率和水平地表非朗伯体假设情况的大气上界辐射亮度;
基于所述水平地表的大气辐射传输模型和所述倾斜地表的二向反射率,获取所述倾斜地表的地表反射率;
基于所述倾斜地表的地表反射率、所述第一核驱动模型和所述倾斜地表的二向反射率,得到所述第二核驱动模型中各向同性系数、体散射核系数和几何散射核系数的关系表达式;
基于所述第二核驱动模型中各向同性系数、体散射核系数和几何散射核系数的关系表达式和所述多时相高分辨率山地卫星图像,得到所述第二核驱动模型;
采用所述第二核驱动模型对所述多时相高分辨率山地卫星图像进行地形辐射校正。
2.根据权利要求1所述的一种多时相高分辨率山地卫星影像地形辐射校正方法,其特征在于:所述基于所述第一核驱动模型,将所述水平地表与倾斜地表的角度关系进行归一化,包括:
基于所述第一核驱动模型,获取水平地表条件下的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角的角度关系表达式;
获取倾斜地表的坡角数据和坡向数据;
基于所述水平地表条件下的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角的角度关系表达式及所述倾斜地表的坡角数据和坡向数据,得到所述倾斜地表条件下太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角的归一化角度关系表达式。
3.根据权利要求2所述的一种多时相高分辨率山地卫星影像地形辐射校正方法,其特征在于:所述基于所述水平地表的大气辐射传输模型和所述倾斜地表的二向反射率,获取所述倾斜地表的地表反射率,包括:
基于所述水平地表朗伯体假设情况的大气上界反射率、水平地表非朗伯体假设情况的大气上界反射率和所述倾斜地表的二向反射率,获取所述倾斜地表的大气上界辐射表达式;
基于所述水平地表非朗伯体假设情况的大气上界辐射亮度,获取所述大气上界辐射亮度经倾斜地表重新分配后,所述倾斜地表的太阳辐射与所述水平地表的太阳辐射之间的相互关系式;所述倾斜地表的大气上界辐射亮度包括太阳直射辐射、天空散射辐射和临近地表反射辐射;基于所述相互关系式,获取所述倾斜地表的大气上界辐射亮度中太阳直射辐射、天空散射辐射和临近地表反射辐射与所述倾斜地表受到的太阳总辐射之间的比例关系式;
基于所述倾斜地表的大气上界辐射表达式、相互关系式和所述比例关系式,获取所述倾斜地表的地表反射率。
4.根据权利要求3所述的一种多时相高分辨率山地卫星影像地形辐射校正方法,其特征在于:所述基于所述相互关系式,获取所述倾斜地表的大气上界辐射亮度中太阳直射辐射、天空散射辐射和临近地表反射辐射与所述倾斜地表受到的太阳总辐射之间的比例关系式,包括:基于所述归一化角度关系表达式,获取所述太阳辐射的地形可见因子系数;
获取所述倾斜地表的邻近地表反射率;
获取所述太阳辐射的大气透过率系数;
基于所述地形可见因子系数、邻近地表反射率、大气透过率系数和所述相互关系式,获取所述比例关系式。
5.根据权利要求4所述的一种多时相高分辨率山地卫星影像地形辐射校正方法,其特征在于:所述基于所述倾斜地表的大气上界辐射表达式、相互关系式和所述比例关系式,获取所述倾斜地表的地表反射率,包括:
联立所述倾斜地表的大气上界辐射表达式、相互关系式和所述比例关系式,获取所述倾斜地表的地表反射率表达式。
6.根据权利要求5所述的一种多时相高分辨率山地卫星影像地形辐射校正方法,其特征在于:所述基于所述倾斜地表的地表反射率、所述第一核驱动模型和所述倾斜地表的二向反射率,得到所述第二核驱动模型中各向同性系数、体散射核系数和几何散射核系数的关系表达式,包括:
将所述倾斜地表的二向反射率带入所述倾斜地表的地表反射率表达式中,并获取所述地表反射率的一阶近似;
将所述倾斜地表的坡角数据、坡向数据和地表反射率的一阶近似带入所述地表反射率表达式中,获取所述各向同性系数、体散射核系数和几何散射核系数的关系表达式。
7.根据权利要求6所述的一种多时相高分辨率山地卫星影像地形辐射校正方法,其特征在于:所述基于所述第二核驱动模型中各向同性系数、体散射核系数和几何散射核系数的关系表达式和所述多时相高分辨率山地卫星图像,得到所述第二核驱动模型,包括:
获取所述多时相高分辨率山地卫星图像中的大气上界反射率、太阳入射天顶角、太阳入射方位角、卫星观测天顶角、卫星观测方位角、坡角和坡向;
将所述多时相高分辨率山地卫星图像中的大气上界反射率、太阳入射天顶角、太阳入射方位角、卫星观测天顶角、卫星观测方位角、坡角和坡向带入所述各向同性系数、体散射核系数和几何散射核系数的关系表达式,采用最小二乘法对所述各向同性系数、体散射核系数和几何散射核系数进行拟合,得到所述各向同性系数、体散射核系数和几何散射核系数;
将所述各向同性系数、体散射核系数和几何散射核系数带入所述第一核驱动模型,得到所述第二核驱动模型。
8.根据权利要求7所述的一种多时相高分辨率山地卫星影像地形辐射校正方法,其特征在于:所述采用所述第二核驱动模型对所述多时相高分辨率山地卫星图像进行地形辐射校正,包括:
获取所述多时相高分辨率山地卫星图像中每个像元倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角;基于所述第二核驱动模型和所述每个像元中倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角,获取每个像元所对应的地表二向反射率;
基于所述每个像元所对应的地表二向反射率,对所述多时相高分辨率山地卫星图像进行同一太阳观测角度和传感器观测角度的归一化协同校正。
9.根据权利要求8所述的一种多时相高分辨率山地卫星影像地形辐射校正方法,其特征在于:所述获取所述多时相高分辨率山地卫星图像中每个像元倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角,包括:
逐个像元获取所述多时相高分辨率山地卫星图像中水平地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角;获取所述多时相高分辨率山地卫星图像中每个所述像元的地理纬度、地理经度、赤纬、时角和卫星经度;
基于所述多时相高分辨率山地卫星图像中每个像元水平地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角及所述多时相高分辨率山地卫星图像中每个所述像元的地理纬度、地理经度、赤纬、时角和卫星经度,获取所述多时相高分辨率山地卫星图像中每个像元倾斜地表的太阳入射天顶角、太阳入射方位角、卫星观测天顶角和卫星观测方位角。
CN202211181020.4A 2022-09-27 2022-09-27 一种多时相高分辨率山地卫星影像地形辐射校正方法 Active CN115524763B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211181020.4A CN115524763B (zh) 2022-09-27 2022-09-27 一种多时相高分辨率山地卫星影像地形辐射校正方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211181020.4A CN115524763B (zh) 2022-09-27 2022-09-27 一种多时相高分辨率山地卫星影像地形辐射校正方法

Publications (2)

Publication Number Publication Date
CN115524763A CN115524763A (zh) 2022-12-27
CN115524763B true CN115524763B (zh) 2023-08-11

Family

ID=84700437

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211181020.4A Active CN115524763B (zh) 2022-09-27 2022-09-27 一种多时相高分辨率山地卫星影像地形辐射校正方法

Country Status (1)

Country Link
CN (1) CN115524763B (zh)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102338871A (zh) * 2010-07-22 2012-02-01 曹春香 地表反射率计算方法和装置
CN102902883A (zh) * 2012-09-24 2013-01-30 北京师范大学 一种基于多角度测量构建二向性反射分布函数(brdf)原型库的方法
CN103413014A (zh) * 2013-03-11 2013-11-27 北京师范大学 一种基于二向性反射分布函数(brdf)原型反演地表反照率的方法
CN104834814A (zh) * 2015-04-29 2015-08-12 西北师范大学 遥感影像地形标准化方法
CN105718724A (zh) * 2016-01-19 2016-06-29 北京师范大学 一种基于核驱动二向反射分布函数模型的天空散射光校正方法
CN111563228A (zh) * 2020-05-07 2020-08-21 中国科学院、水利部成都山地灾害与环境研究所 基于地表入射短波辐射的山地地表反射率地形改正方法
CN113155740A (zh) * 2020-01-07 2021-07-23 国家卫星气象中心(国家空间天气监测预警中心) 一种定标基准场brdf特性分析方法及系统
CN113672847A (zh) * 2021-08-18 2021-11-19 滁州学院 一种基于卫星遥感数据的积雪多角度二向反射率反演方法

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7593835B2 (en) * 2001-04-20 2009-09-22 Spectral Sciences, Inc. Reformulated atmospheric band model method for modeling atmospheric propagation at arbitrarily fine spectral resolution and expanded capabilities.
US10600162B2 (en) * 2016-12-29 2020-03-24 Konica Minolta Laboratory U.S.A., Inc. Method and system to compensate for bidirectional reflectance distribution function (BRDF)

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102338871A (zh) * 2010-07-22 2012-02-01 曹春香 地表反射率计算方法和装置
CN102902883A (zh) * 2012-09-24 2013-01-30 北京师范大学 一种基于多角度测量构建二向性反射分布函数(brdf)原型库的方法
CN103413014A (zh) * 2013-03-11 2013-11-27 北京师范大学 一种基于二向性反射分布函数(brdf)原型反演地表反照率的方法
CN104834814A (zh) * 2015-04-29 2015-08-12 西北师范大学 遥感影像地形标准化方法
CN105718724A (zh) * 2016-01-19 2016-06-29 北京师范大学 一种基于核驱动二向反射分布函数模型的天空散射光校正方法
CN113155740A (zh) * 2020-01-07 2021-07-23 国家卫星气象中心(国家空间天气监测预警中心) 一种定标基准场brdf特性分析方法及系统
CN111563228A (zh) * 2020-05-07 2020-08-21 中国科学院、水利部成都山地灾害与环境研究所 基于地表入射短波辐射的山地地表反射率地形改正方法
CN113672847A (zh) * 2021-08-18 2021-11-19 滁州学院 一种基于卫星遥感数据的积雪多角度二向反射率反演方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
Ainong Li等.An Improved Physics-Based Model for Topographic Correction of Landsat TM Images .Remote Sens..2015,第6296-6319页. *

Also Published As

Publication number Publication date
CN115524763A (zh) 2022-12-27

Similar Documents

Publication Publication Date Title
Delacourt et al. Remote-sensing techniques for analysing landslide kinematics: a review
Su et al. A high-precision aerosol retrieval algorithm (HiPARA) for advanced Himawari imager (AHI) data: Development and verification
Vander Jagt et al. Snow depth retrieval with UAS using photogrammetric techniques
Engdahl et al. Land-cover classification using multitemporal ERS-1/2 InSAR data
CN108132220B (zh) 林区机载推扫式高光谱影像的brdf归一化校正方法
Ni et al. Mapping three-dimensional structures of forest canopy using UAV stereo imagery: Evaluating impacts of forward overlaps and image resolutions with LiDAR data as reference
CN107917881B (zh) 基于冠层内辐射传输机理的光学遥感影像地形校正方法
CN108765488B (zh) 一种基于阴影的高分辨率遥感影像建筑物高度估测方法
CN107389617A (zh) 基于高分四号卫星的气溶胶光学厚度的反演方法及设备
Seiz et al. Cloud mapping from the ground: Use of photogrammetric methods
CN113324915B (zh) 一种支持高分辨率气溶胶光学厚度反演的城市复杂地表反射率估算方法
Rauchmiller Jr et al. Measurement of the Landsat Thematic Mapper modulation transfer function using an array of point sources
Jing et al. Sub-pixel accuracy evaluation of FY-3D MERSI-2 geolocation based on OLI reference imagery
Shen et al. Antarctic-wide annual ice flow maps from Landsat 8 imagery between 2013 and 2019
ZYLSHAL et al. CORRECTING THE TOPOGRAPHIC EFFECT ON SPOT-6/7 MULTISPECTRAL IMAGERIES: A COMPARISON OF DIFFERENT DIGITAL ELEVATION MODELS.
CN115524763B (zh) 一种多时相高分辨率山地卫星影像地形辐射校正方法
Wang et al. Removing temperature drift and temporal variation in thermal infrared images of a UAV uncooled thermal infrared imager
CA2613252A1 (en) Method and apparatus for determining a location associated with an image
Lingua et al. Iterative Refraction-Correction Method on Mvs-Sfm for Shallow Stream Bathymetry
Cariou et al. Automatic georeferencing of airborne pushbroom scanner images with missing ancillary data using mutual information
Rau et al. Landslide deformation monitoring by three-camera imaging system
Pan et al. Estimating Cloud-Free Fractional Snow Cover from Himawari-8, FY-4A and Modis Observation
Small et al. Combination of ascending/descending ERS-1 InSAR data for calibration and validation
Bostater Jr Multiresolution earth remote sensing approach
Yu et al. Atmospheric correction of HJ-1A multi-spectral and hyper-spectral images

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