CN115861629A - 一种高分耕地影像提取方法 - Google Patents

一种高分耕地影像提取方法 Download PDF

Info

Publication number
CN115861629A
CN115861629A CN202211563055.4A CN202211563055A CN115861629A CN 115861629 A CN115861629 A CN 115861629A CN 202211563055 A CN202211563055 A CN 202211563055A CN 115861629 A CN115861629 A CN 115861629A
Authority
CN
China
Prior art keywords
image
farmland
ndvi
resolution
time sequence
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
Application number
CN202211563055.4A
Other languages
English (en)
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.)
Nanjing Surveying And Mapping Research Institute Co ltd
Original Assignee
Nanjing Surveying And Mapping Research Institute Co ltd
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 Nanjing Surveying And Mapping Research Institute Co ltd filed Critical Nanjing Surveying And Mapping Research Institute Co ltd
Priority to CN202211563055.4A priority Critical patent/CN115861629A/zh
Publication of CN115861629A publication Critical patent/CN115861629A/zh
Pending legal-status Critical Current

Links

Images

Abstract

一种高分耕地影像提取方法,包括(1)合并现有中等空间分辨率耕地遥感数据产品作为初始耕地空间范围,综合植被指数和谐波时序分析因子对产品优化,去除显著非耕地区域;(2)采用前一步成果为样本库,基于时序中分影像构建空间‑光谱‑时序特征集,并通过随机森林算法重新提取耕地;(3)对高空间分辨率影像进行多尺度分割,获取分割矢量斑块,将上一步成果嵌套至矢量斑块中,提取斑块均值和周长面积比,对其自然断点法处理获取高置信度的高空间分辨率样本;(4)将迁移获取的高分样本输入至语义分割模型以提取耕地;本方法利用现有遥感产品为先验信息,有效去除产品中错误信息,降低人工标注样本成本,自动化程度高。

Description

一种高分耕地影像提取方法
技术领域
本发明涉及遥感信息提取技术领域,特别涉及一种高分耕地影像提取方法。
背景技术
耕地是指种植农作物的土地,包括熟地、新开发、复垦、整理地,休闲地(含轮歇地、休耕地);以种植农作物(含蔬菜)为主,间有零星果树、桑树或其他树木的土地;平均每年能保证收获一季的已垦滩地和海涂。耕地与经济发展和人类生活密切相关,及时可靠的耕地信息在作物产量估算、生态环境等研究中尤为重要,已引起世界各国广泛关注,而我国对耕地的管理也极为重视,因此,获取具有时效性强、精度高的耕地分布对于耕地保护和可持续发展意义重大。近年来,由于经济的快速发展和工业化、城市化进程的快速推进,使得耕地在部分区域呈现非农化、非粮化现象。随着遥感技术的不断进步,遥感影像的空间、光谱、时间分辨率均在不断提升,特别是Landsat、Sentinel等数据的免费共享,为耕地信息的准确提取提供了数据基础。
当前,多种包含耕地的遥感土地覆盖产品主要集中于中等或粗空间分辨率,对于高空分辨率的耕地产品较少。监督分类是常用的一种可靠地物信息提取方案,充足和可信的训练样本是其提高识别准确性的关键。在大多数研究中,常常通过目视解译和实地调查来收集地物样本。然而,大范围获取足够样本通常耗时费力。因此,迫切需要开发用于自动获取样本的方法,以支持大规模耕地绘图,尤其是用于支撑高空间分辨率图像耕地识别的研究。
发明内容
本发明针对现有技术中的不足,提供一种高分耕地影像提取方法;以解决目前遥感测绘领域中对高空分辨率的耕地产品识别较少的问题,以及减少大范围获取样本耗时费力的问题,提供自动获取样本的方法以支持大规模耕地绘图,尤其是用于支撑高空间分辨率图像耕地识别的研究。
为实现上述目的,本发明采用以下技术方案:
一种高分耕地影像提取方法,包括以下步骤:
S1、获取现有的10米空间分辨率遥感数据产品,提取耕地附近区域的影像后取并集,得到初始的先验耕地范围;
S2、获取10米空间分辨率时序遥感影像,计算时序归一化植被指数NDVI;
S3、基于步骤S2中获取的时序归一化植被指数NDVI,分别计算其标准差NDVIstd和谐波时序分量系数;
S4、对步骤S3中的NDVIstd标准差和谐波时序分量系数分别进行高斯分布聚类,得到聚类后的耕地区域和非耕地区域,以此校正步骤S1中的先验耕地范围,得到10米优化耕地产品,该产品中包含有校正后的耕地区域影像和校正后的非耕地区域影像;
S5、基于步骤S2中的10米空间分辨率时序遥感影像,分别提取该影像中的纹理特征、光谱特征、时序特征,组成集合,进而构建多类型影像特征集;
S6、对步骤S4中10米优化耕地产品进行分层随机抽样获取耕地区域影像的样本和非耕地区域影像的样本,整合步骤S5中多类型影像特征集,将两者输入随机森林分类器,进而分类提取10米重分类耕地产品,该产品中包含有具备多类型影像特征的耕地区域影像和非耕地区域影像;
S7、获取高空间分辨率遥感影像,并以步骤S6中所获的10米重分类耕地产品为源样本,对高空间分辨率遥感影像进行多尺度分割获取矢量斑块,将两者进行叠置分析分别求取矢量斑块的均值和周长面积比这两个指标,利用自然断点法对指标进行阈值分割,得到高空间分辨率的耕地影像样本和非耕地影像样本;
S8、将步骤S7中高空间分辨率的耕地样本和非耕地样本输入至深度学习语义分割模型中进行训练,训练完成后用于对目标区范围进行识别,以提取高空间分辨率耕地影像。
为优化上述技术方案,采取的具体措施还包括:
进一步地,在步骤S1中所用10米空间分辨率遥感数据产品是ESAWorldCover和Esri Land Cover。
进一步地,
在步骤S2中,所用10米空间分辨率时序遥感影像是基于Sentinel-2卫星获取的,且除Coastal aerosol、water vapor和SWIR-Cirrus外的影像波段,其余波段均使用;
在步骤S2中,时序归一化植被指数NDVI计算公式如下:
Figure BDA0003985330910000021
式中,ρNIR表示近红外波段反射率;ρR表示红光波段反射率。
进一步地,步骤S3的具体内容为:
标准差NDVIstd的计算公式如下:
Figure BDA0003985330910000022
式中,NDVImean表示多次NDVI观测值的平均值;N表示多次NDVI观测值的有效数;NDVIn表示第n次的NDVI观测值;
谐波时序分量系数ai和bi的计算公式如下:
Figure BDA0003985330910000031
式中,YNDVI(tj)表示NDVI的时间序列,ε(tj)表示NDVI的误差序列,其中tj表示获取NDVI的时间节点;a0表示主分量,为常量;i表示NDVI时间序列中的周期项个数;f表示频率;ai和bi表示为谐波时序分量系数。
进一步地,在步骤S4中,所述校正步骤S1中的先验耕地范围的具体内容为:以聚类后的耕地区域范围为基准,对比先验耕地范围中与之不匹配的区域,以此确定校正后的耕地区域影像和校正后的非耕地区域影像。
进一步地,在步骤S5中,所述分别提取影像中的纹理特征、光谱特征、时序特征,组成集合,构建多类型影像特征集的具体内容为:
构建的影像纹理特征包括:方差、对比度、相关性,采用灰度共生矩阵提取方法;
构建的影像光谱特征包括:归一化植被指数、归一化水体指数、影像中红、绿、蓝、近红外、红边、短波红外的波段;
构建的时序特征包括:NDVI标准差、NDVI时序谐波分量系数、NDVI极差、NDVI最小值;
以此组成集合形成多类型影像特征集。
进一步地,在步骤S7中,
所述高空间分辨率遥感影像的空间分辨率为1米;
所述多尺度分割算法参数设置包括:尺度100、形状0.1、紧致度0.5;
所述执行阈值分割时阈值由自然断点法确定,当斑块均值大于阈值时为耕地,周长面积比大于阈值时为田埂或道路。
进一步地,在步骤S8中,
所述深度学习语义分割模型采用的是DeepLab V3+;
所述深度学习语义分割模型在训练之前,将高空间分辨率的样本划分至训练集、验证集和测试集中,且划分比例为5:1:4,以用于深度学习语义分割模型的训练、验证和测试;
所述深度学习语义分割模型的优化器设置为SGD,批量大小BatchSize设置为50,训练迭代次数设置为50,学习率lr设置为10-3,损失函数为交叉熵损失;
所述深度学习语义分割模型训练结束后对目标区范围进行识别时,其批量大小BatchSize设置为50。
一种计算机可读存储介质,存储有计算机程序,所述计算机程序使计算机执行如上述任一项所述的高分耕地影像提取方法。
一种电子设备,包括:存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行计算机程序时,实现如上述任一项所述的高分耕地影像提取方法。
本发明的有益效果是:
(1)提出了借用并优化当前耕地产品的解决方法。充分利用现有耕地产品的先验信息,通过作物时序信息和耕地特征,排除其中的错误信息,提升耕地产品的质量,大大降低了人工收集地物样本的成本。
(2)设计了样本跨尺度迁移的模块。以中等空间分辨率耕地提取结果为源样本,利用图像分割思想获取矢量斑块来对其进行规范,设计指标提高其在高空间分辨率影像上的置信度,进而完成样本的跨尺度迁移。
(3)本申请在高空间分辨率耕地的自动化提取方面取得了良好的效果,相比于传统的基于目视解译获取样本进行耕地识别的深度学习方案而言,运行效率高,为基于高空间分辨率遥感影像的耕地识别提供了新的解决思路。
附图说明
图1是本申请整体技术方案流程示意图。
图2是本申请实施例中的区域影像以及人工标注后影像的示意图。
图3是本申请实施中训练后的深度学习语义分割模型所实际预测识别的影像示意图。
具体实施方式
现在结合附图对本发明作进一步详细的说明。
本发明实例围绕获取自山东淄博的高空间分辨率遥感影像,其空间分辨率为1米,包含RGB三个波段,影像大小为12951×9462,获取日期为2020年5月12日。以此为基准,获取该区域相应包含耕地的10米空间分辨率遥感公共产品,分别为ESAWorldCover 2020和EsriLand Cover 2020。同时,下载获取覆盖该区域的2020年全年的Sentinel-2时序影像,空间分辨率也为10米。
本实施案例总体流程图如图1所示,包括以下步骤:
第一步:数据准备及预处理。获取山东淄博研究区的数据集,具体包括:遥感公共产品、Sentinel-2时序影像和高空间分辨率遥感影像。使用地理信息分析软件ArcGIS分别对影像进行投影、配准和裁剪,以使不同数据在地理空间进行对齐。同时人工标注该区域耕地和非耕地样本,用于耕地识别结果的精度验证。该区域RGB影像及其对应的地表真实标注如图2所示(其中(a)是试验区RGB影像,(b)是真值标签)。此外,提取10米空间分辨率遥感公共产品耕地类取并集,获得初始的先验耕地范围。
第二步:计算时序归一化植被指数NDVI。基于10米分辨率的Sentinel-2时序遥感影像,首先进行形态特征的提取,计算时序归一化植被指数NDVI,计算公式如下:
Figure BDA0003985330910000051
式中,ρR为红光波段反射率,ρNIR为近红外波段反射率。
第三步:计算NDVI标准差和NDVI谐波时序分量系数。基于第二步获取的时序NDVI,分别计算获取其标准差NDVIstd和谐波时序分量系数ai和bi
第四步:高斯二项分布聚类获取优化耕地产品。将第三步获取的结果分别进行高斯分布聚类,获取聚类后的耕地和非耕地区,以聚类后的耕地范围为基准,掩膜去除第一步所获先验耕地范围中与之不匹配的区域,进而获取10米优化耕地产品。
第五步:构建多类型影像特征集。分别计算提取影像多类型特征,采用灰度共生矩阵提取纹理特征:方差、对比度、相关性;构建影像的光谱特征:归一化植被指数、改进的归一化水体指数、中值和四分之三分位数合成的影像红、绿、蓝、近红外、红边、短波红外1和短波红外2波段;构建的时序特征包括:NDVI标准差、NDVI时序谐波分量系数、NDVI极值差、NDVI最小值。在此基础上,组建多类型影像特征集。
第六步:重分类获取10米耕地。对第四步中10米优化耕地产品进行分层随机抽样获取耕地和非耕地样本,整合第五步中多类型影像特征集,将两者输入随机森林分类器,随机森林算法参数设置包括:分类树的数量为500,子集中的特征数量设置为特征数量值的平方根。最终分类提取10米重分类耕地。
第七步:样本的跨尺度迁移。以第六步中所获10米重分类耕地为源样本,对1米空间分辨率遥感影像进行多尺度分割获取矢量斑块,其中,多尺度分割算法参数设置包括:尺度100、形状0.1、紧致度0.5。将两者进行叠置分析分别求取斑块均值和周长面积比两个指标,利用自然断点法对指标进行阈值分割,当斑块均值大于自然断点法识别的阈值时为耕地,周长面积比大于阈值时为田埂或道路,将其用于对耕地掩膜,去除其中的非耕地部分,进而以此迁移获取高空间分辨率耕地样本。
第八步:高空间分辨率耕地提取。将第七步中高空间分辨率样本输入至DeepLabV3+中,在模型训练阶段,输入数据为带有标签的1米空间分辨率遥感影像,输入影像大小为128×128,训练集、验证集和测试集的划分比例为5:1:4,优化器设置为SGD,批量大小BatchSize设置为50,训练迭代次数设置为50,学习率lr设置为10-3,损失函数为交叉熵损失(CrossEntropy Loss)。预测阶段:输入影像大小128×128,批量大小BatchSize设置为50。在完成预测后,测试集的总体精度为93.80%,F1-score 94.16%,识别结果局部细节如图3所示,其中图3中的(a)和(b)是试验区的部分实际影像,(a1)和(b1)是分别对应(a)和(b)的真值标签,即人工标注后的示意图,(a2)和(b2)是训练后深度学习语义分割模型对应(a)和(b)所识别后的耕地影像。
需要注意的是,发明中所引用的如“上”、“下”、“左”、“右”、“前”、“后”等的用语,亦仅为便于叙述的明了,而非用以限定本发明可实施的范围,其相对关系的改变或调整,在无实质变更技术内容下,当亦视为本发明可实施的范畴。
以上仅是本发明的优选实施方式,本发明的保护范围并不仅局限于上述实施例,凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出,对于本技术领域的普通技术人员来说,在不脱离本发明原理前提下的若干改进和润饰,应视为本发明的保护范围。

Claims (10)

1.一种高分耕地影像提取方法,其特征在于,包括以下步骤:
S1、获取现有的10米空间分辨率遥感数据产品,提取耕地附近区域的影像后取并集,得到初始的先验耕地范围;
S2、获取10米空间分辨率时序遥感影像,计算时序归一化植被指数NDVI;
S3、基于步骤S2中获取的时序归一化植被指数NDVI,分别计算其标准差NDVIstd和谐波时序分量系数;
S4、对步骤S3中的NDVIstd标准差和谐波时序分量系数分别进行高斯分布聚类,得到聚类后的耕地区域和非耕地区域,以此校正步骤S1中的先验耕地范围,得到10米优化耕地产品,该产品中包含有校正后的耕地区域影像和校正后的非耕地区域影像;
S5、基于步骤S2中的10米空间分辨率时序遥感影像,分别提取该影像中的纹理特征、光谱特征、时序特征,组成集合,进而构建多类型影像特征集;
S6、对步骤S4中10米优化耕地产品进行分层随机抽样获取耕地区域影像的样本和非耕地区域影像的样本,整合步骤S5中多类型影像特征集,将两者输入随机森林分类器,进而分类提取10米重分类耕地产品,该产品中包含有具备多类型影像特征的耕地区域影像和非耕地区域影像;
S7、获取高空间分辨率遥感影像,并以步骤S6中所获的10米重分类耕地产品为源样本,对高空间分辨率遥感影像进行多尺度分割获取矢量斑块,将两者进行叠置分析分别求取矢量斑块的均值和周长面积比这两个指标,利用自然断点法对指标进行阈值分割,得到高空间分辨率的耕地影像样本和非耕地影像样本;
S8、将步骤S7中高空间分辨率的耕地样本和非耕地样本输入至深度学习语义分割模型中进行训练,训练完成后用于对目标区范围进行识别,以提取高空间分辨率耕地影像。
2.根据权利要求1所述的一种高分耕地影像提取方法,其特征在于,在步骤S1中所用10米空间分辨率遥感数据产品是ESAWorldCover和Esri Land Cover。
3.根据权利要求1所述的一种高分耕地影像提取方法,其特征在于,
在步骤S2中,所用10米空间分辨率时序遥感影像是基于Sentinel-2卫星获取的,且除Coastal aerosol、water vapor和SWIR-Cirrus外的影像波段,其余波段均使用;
在步骤S2中,时序归一化植被指数NDVI计算公式如下:
Figure FDA0003985330900000011
式中,ρNIR表示近红外波段反射率;ρR表示红光波段反射率。
4.根据权利要求1所述的一种高分耕地影像提取方法,其特征在于,步骤S3的具体内容为:
标准差NDVIstd的计算公式如下:
Figure FDA0003985330900000021
式中,NDVImean表示多次NDVI观测值的平均值;N表示多次NDVI观测值的有效数;NDVIn表示第n次的NDVI观测值;
谐波时序分量系数ai和bi的计算公式如下:
Figure FDA0003985330900000022
式中,YNDVI(tj)表示NDVI的时间序列,ε(tj)表示NDVI的误差序列,其中tj表示获取NDVI的时间节点;a0表示主分量,为常量;i表示NDVI时间序列中的周期项个数;f表示频率;ai和bi表示为谐波时序分量系数。
5.根据权利要求1所述的一种高分耕地影像提取方法,其特征在于,在步骤S4中,所述校正步骤S1中的先验耕地范围的具体内容为:以聚类后的耕地区域范围为基准,对比先验耕地范围中与之不匹配的区域,以此确定校正后的耕地区域影像和校正后的非耕地区域影像。
6.根据权利要求1所述的一种高分耕地影像提取方法,其特征在于,在步骤S5中,所述分别提取影像中的纹理特征、光谱特征、时序特征,组成集合,构建多类型影像特征集的具体内容为:
构建的影像纹理特征包括:方差、对比度、相关性,采用灰度共生矩阵提取方法;
构建的影像光谱特征包括:归一化植被指数、归一化水体指数、影像中红、绿、蓝、近红外、红边、短波红外的波段;
构建的时序特征包括:NDVI标准差、NDVI时序谐波分量系数、NDVI极差、NDVI最小值;
以此组成集合形成多类型影像特征集。
7.根据权利要求1所述的一种高分耕地影像提取方法,其特征在于,在步骤S7中,
所述高空间分辨率遥感影像的空间分辨率为1米;
所述多尺度分割算法参数设置包括:尺度100、形状0.1、紧致度0.5;
所述执行阈值分割时阈值由自然断点法确定,当斑块均值大于阈值时为耕地,周长面积比大于阈值时为田埂或道路。
8.根据权利要求1所述的一种高分耕地影像提取方法,其特征在于,在步骤S8中,
所述深度学习语义分割模型采用的是DeepLab V3+;
所述深度学习语义分割模型在训练之前,将高空间分辨率的样本划分至训练集、验证集和测试集中,且划分比例为5:1:4,以用于深度学习语义分割模型的训练、验证和测试;
所述深度学习语义分割模型的优化器设置为SGD,批量大小BatchSize设置为50,训练迭代次数设置为50,学习率lr设置为10-3,损失函数为交叉熵损失;
所述深度学习语义分割模型训练结束后对目标区范围进行识别时,其批量大小BatchSize设置为50。
9.一种计算机可读存储介质,存储有计算机程序,其特征在于,所述计算机程序使计算机执行如权利要求1-8任一项所述的高分耕地影像提取方法。
10.一种电子设备,其特征在于,包括:存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序,所述处理器执行计算机程序时,实现如权利要求1-8任一项所述的高分耕地影像提取方法。
CN202211563055.4A 2022-12-07 2022-12-07 一种高分耕地影像提取方法 Pending CN115861629A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202211563055.4A CN115861629A (zh) 2022-12-07 2022-12-07 一种高分耕地影像提取方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202211563055.4A CN115861629A (zh) 2022-12-07 2022-12-07 一种高分耕地影像提取方法

Publications (1)

Publication Number Publication Date
CN115861629A true CN115861629A (zh) 2023-03-28

Family

ID=85670659

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202211563055.4A Pending CN115861629A (zh) 2022-12-07 2022-12-07 一种高分耕地影像提取方法

Country Status (1)

Country Link
CN (1) CN115861629A (zh)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117132897A (zh) * 2023-07-14 2023-11-28 西北大学 一种基于Landsat和MODIS的自动选取棉花样本的提取方法
CN117690017A (zh) * 2023-11-16 2024-03-12 宁波大学 一种顾及物候时序特征的单、双季水稻提取方法

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN117132897A (zh) * 2023-07-14 2023-11-28 西北大学 一种基于Landsat和MODIS的自动选取棉花样本的提取方法
CN117690017A (zh) * 2023-11-16 2024-03-12 宁波大学 一种顾及物候时序特征的单、双季水稻提取方法

Similar Documents

Publication Publication Date Title
Chen et al. Mapping croplands, cropping patterns, and crop types using MODIS time-series data
CN111598019B (zh) 基于多源遥感数据的作物类型与种植模式识别方法
CN111368736B (zh) 基于sar和光学遥感数据的水稻精细化估产方法
CN115861629A (zh) 一种高分耕地影像提取方法
CN112164062A (zh) 一种基于遥感时序分析的抛荒地信息提取方法及装置
CN111462223B (zh) 基于Sentinel-2影像的江淮地区大豆和玉米种植面积识别方法
Chen et al. Predicting individual apple tree yield using UAV multi-source remote sensing data and ensemble learning
CN114821362A (zh) 一种基于多源数据的水稻种植面积提取方法
CN116543316B (zh) 一种利用多时相高分辨率卫星影像识别稻田内草皮的方法
CN114926748A (zh) Sentinel-1/2微波与光学多光谱影像结合的大豆遥感识别方法
WO2023109652A1 (zh) 水稻种植提取及复种指数监测方法、系统、终端及存储介质
CN107680098A (zh) 一种甘蔗蔗节特征的识别方法
Zeng et al. A phenology-based vegetation index classification (PVC) algorithm for coastal salt marshes using Landsat 8 images
CN107121681B (zh) 基于高分卫星遥感数据的居民地提取系统
Chen et al. Effect of crop spectra purification on plant nitrogen concentration estimations performed using high-spatial-resolution images obtained with unmanned aerial vehicles
Qiu et al. Automated cropping intensity extraction from isolines of wavelet spectra
Agustin et al. Deep Learning-based Method for Multi-Class Classification of Oil Palm Planted Area on Plant Ages Using Ikonos Panchromatic Imagery
CN116310864A (zh) 一种作物倒伏自动识别方法、系统、电子设备及介质
CN115049900A (zh) 一种基于光学与雷达数据的水稻种植区识别方法及系统
Zhang et al. A Mapping Approach for Eucalyptus Plantations Canopy and Single-Tree Using High-Resolution Satellite Images in Liuzhou, China
CN113096114B (zh) 结合形态与指数的高分辨率城市水体图斑遥感提取方法
Song et al. Hierarchical extraction of cropland boundaries using Sentinel-2 time-series data in fragmented agricultural landscapes
Zulkafli et al. Tree crown density analysis from hyperspectral image
Song et al. Crop Classification Using Multitemporal Landsat 8 Images
CN115063610B (zh) 基于Sentinel-1、2影像的大豆种植区识别方法

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