CN113327197B - 遥感影像时空融合方法、智能终端及计算机可读存储介质 - Google Patents
遥感影像时空融合方法、智能终端及计算机可读存储介质 Download PDFInfo
- Publication number
- CN113327197B CN113327197B CN202110507192.5A CN202110507192A CN113327197B CN 113327197 B CN113327197 B CN 113327197B CN 202110507192 A CN202110507192 A CN 202110507192A CN 113327197 B CN113327197 B CN 113327197B
- Authority
- CN
- China
- Prior art keywords
- resolution image
- pixel
- image
- low
- fusion
- 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
Links
- 238000007500 overflow downdraw method Methods 0.000 title claims abstract description 28
- 230000004927 fusion Effects 0.000 claims abstract description 127
- 238000000034 method Methods 0.000 claims abstract description 51
- 238000001914 filtration Methods 0.000 claims abstract description 22
- 238000012216 screening Methods 0.000 claims abstract description 9
- 238000001228 spectrum Methods 0.000 claims description 34
- 230000003595 spectral effect Effects 0.000 claims description 33
- 230000006870 function Effects 0.000 claims description 12
- 238000004364 calculation method Methods 0.000 claims description 11
- 238000005070 sampling Methods 0.000 claims description 11
- 238000012937 correction Methods 0.000 claims description 10
- 230000003044 adaptive effect Effects 0.000 claims description 8
- 238000012545 processing Methods 0.000 claims description 8
- 238000004458 analytical method Methods 0.000 claims description 6
- 238000005516 engineering process Methods 0.000 claims description 6
- 238000012417 linear regression Methods 0.000 claims description 6
- 238000012935 Averaging Methods 0.000 claims description 3
- 238000001514 detection method Methods 0.000 claims description 3
- 238000007781 pre-processing Methods 0.000 claims description 3
- 230000005855 radiation Effects 0.000 claims description 3
- 238000010586 diagram Methods 0.000 description 9
- 238000004422 calculation algorithm Methods 0.000 description 4
- 238000004590 computer program Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000012821 model calculation Methods 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- OIGNJSKKLXVSLS-VWUMJDOOSA-N prednisolone Chemical compound O=C1C=C[C@]2(C)[C@H]3[C@@H](O)C[C@](C)([C@@](CC4)(O)C(=O)CO)[C@@H]4[C@@H]3CCC2=C1 OIGNJSKKLXVSLS-VWUMJDOOSA-N 0.000 description 2
- 230000001360 synchronised effect Effects 0.000 description 2
- 230000009897 systematic effect Effects 0.000 description 2
- GNFTZDOKVXKIBK-UHFFFAOYSA-N 3-(2-methoxyethoxy)benzohydrazide Chemical compound COCCOC1=CC=CC(C(=O)NN)=C1 GNFTZDOKVXKIBK-UHFFFAOYSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000000903 blocking effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 238000002474 experimental method Methods 0.000 description 1
- 239000004973 liquid crystal related substance Substances 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 230000003313 weakening effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T3/00—Geometric image transformations in the plane of the image
- G06T3/40—Scaling of whole images or parts thereof, e.g. expanding or contracting
- G06T3/4053—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution
- G06T3/4076—Scaling of whole images or parts thereof, e.g. expanding or contracting based on super-resolution, i.e. the output image resolution being higher than the sensor resolution using the original low-resolution images to iteratively correct the high-resolution images
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
本发明公开了一种遥感影像时空融合方法、智能终端及计算机可读存储存储介质,所述方法包括获取第一时相的高分辨率影像和第一低分辨率影像,以及第二时相的第二低分辨率影像;根据高分辨率影像、第一低分辨率影像和第二低分辨率影像计算像元可靠性系数;利用模糊C均值聚类法对高分辨率影像中的像元进行分类,并根据像元可靠性系数筛选出各类像元参与全局解混,以得到高分辨率影像、第一低分辨率影像和第二低分辨率影像的初始融合影像;基于自适应局部解混模型修正初始融合影像中的粗差,得到中间融合影像;对中间融合影像进行空间滤波和残差分配,得到最终融合影像。综合考虑了异源影像差异对时空融合的影响,取得更优的时空融合结果。
Description
技术领域
本发明涉及遥感领域,尤其涉及的是一种遥感影像时空融合方法、智能终端及计算机可读存储介质。
背景技术
受限于传感器技术、资金,以及云等其他因素对地表观测的影响,利用单一的遥感数据实现连续的高空间分辨率地表观测十分困难,而时空融合技术可以将拥有高空间分辨率但长回访周期的卫星影像和快速回访但低空间分辨率的卫星影像相融合,获得高时空分辨率影像,这也使得时空融合技术广泛的应用于物候分析,灾害监测,城市化分析等领域。
尽管近年来时空融合的方法发展迅速,但由于发生融合的影像之间存在影像对间的配准、异源传感器的性能、观测角度、幅宽、时间不同等异源影像的差异,时空融合的结果并不准确,为了克服这些差异带来的影像,现有技术使用线性回归的方法或通过融合两时相的差值信息来减弱这些差异,但只能一定程度上减小系统误差影响,对配准误差等无效,因此亟需一种有效新策略用于减少各种异源影像差异对时空融合的影响。
发明内容
本发明要解决的技术问题在于,针对现有技术的上述缺陷,提供一种遥感影像时空融合方法,旨在减少各种异源影像差异对时空融合的影响。
本发明解决问题所采用的技术方案如下:
第一方面,本发明实施例提供一种遥感影像时空融合方法,包括:
获取第一时相的高分辨率影像和第一低分辨率影像,以及第二时相的第二低分辨率影像;
根据所述高分辨率影像、第一低分辨率影像和第二低分辨率影像计算像元可靠性系数;
利用模糊C均值聚类法对所述高分辨率影像中的像元进行分类,并根据像元可靠性系数筛选出各类像元参与全局解混,以得到所述高分辨率影像、所述第一低分辨率影像和所述第二低分辨率影像的初始融合影像;
基于自适应局部解混模型修正所述初始融合影像中的粗差,得到中间融合影像;
对所述中间融合影像进行空间滤波和残差分配,得到目标融合影像。
在一种实施方式中,所述根据所述高分辨率影像、第一低分辨率影像和第二分辨率影像中像元光谱值计算像元可靠性系数集的步骤之前,分别对所述高分辨率影像、第一低分辨率影像和第二分辨率影像进行预处理,预处理的方法包括配准和辐射矫正。
在一种实施方式中,所述根据所述高分辨率影像、第一低分辨率影像和第二分辨率影像计算像元可靠性系数的步骤包括:
对所述高分辨率影像中各像元对应的像元光谱值进行降采样处理,得到与所述第一低分辨率影像的分辨率相同的第三低分辨率影像;
对所述第三低分辨率影像中各像元对应的像元光谱值和第一低分辨率影像中各像元对应的像元光谱值进行线性回归分析,得到回归系数;
将所述回归系数代入系统误差矫正函数中以矫正第一低分辨率影像的系统误差,得到第四低分辨率影像;
根据所述第一低分辨率影像、所述第二低分辨率影像、所述第三低分辨率影像和所述第四低分辨率影像计算像元可靠性系数集。
在一种实施方式中,所述根据所述第一低分辨率影像、所述第二低分辨率影像、所述第三低分辨率影像和所述第四低分辨率影像计算像元可靠性系数的步骤包括:
将所述第一低分辨率影像中各像元对应的像元光谱值、所述第二低分辨率影像中各像元对应的像元光谱值、所述第三低分辨率影像中各像元对应的像元光谱值和所述第四低分辨率影像中各像元对应的像元光谱值代入第一公式中,得到阈值系数,其中,所述第一公式为:所述Q(b)为所述阈值系数,所述b为波段号,所述ΔC为所述第一低分辨率影像中各像元与第二低分辨率影像中各像元的像元光谱值的差值组成的集合,所述φdown(F1)为第三低分辨率影像中各像元对应的像元光谱值组成的集合,所述φsystem(C1)为第四低分辨率影像中各像元对应的像元光谱值组成的集合,所述mean表示求平均值函数;
将所述阈值系数代入第二公式,得到像元可靠性系数,其中,所述第二公式为:其中,所述RI(xi,yi,b)表示坐标为(xi,yi)的像元对应的像元可靠性系数,所述φdown[F1(xi,yi,b)]表示所述第三低分辨率影像中坐标为(xi,yi)的像元对应的像元光谱值,所述φsystem[C1(xi,yi,b)]为第四低分辨率影像中坐标为(xi,yi)的像元对应的像元光谱值。
在一种实施方式中,所述利用模糊C均值聚类法对所述高分辨率影像中的像元进行分类,并根据像元可靠性系数筛选出各类像元参与全局解混,以得到所述高分辨率影像、所述第一低分辨率影像和所述第二低分辨率影像的初始融合影像的步骤包括:
设置最大分类值nmax和起始分类值n;
当分类数为n时,利用模糊C均值聚类法将所述高分辨率影像中的像元分成n类,根据所述像元可靠性系数从各类像元中确定用于全局解混的对象像元;
将所述对象像元的丰度输入全局解混模型,得到对象像元在所述第一低分辨率影像和所述第二低分辨率影像之间的两时相平均变化值;
将所述两时相平均变化值加入所述高分辨率影像中,得到融合影像;
将n以1为单元逐级递加,直至n=nmax,得到每一级n对应的融合影像;
计算每一级n对应的融合影像残差值,将所述残差值的最小值对应的融合影像确定为初始融合影像。
在一种实施方式中,所述基于自适应局部解混模型修正所述初始融合影像中的粗差,得到中间融合影像的步骤包括:
对所述初始融合影像进行降采样处理,得到低分辨率融合影像;
根据变化检测技术确定所述低分辨率融合影像中的带粗差像元,并根据所述带粗差像元的平均可靠性系数确定可靠带粗差像元;
以所述可靠带粗差像元为第一中心像元建立预设尺寸的局部窗口;
计算局部窗口中所述可靠带粗差像元与所述第一中心像元的时相变化值,并根据所述时相变化值从所述可靠带粗差像元中选择可靠相似变化像元;
当所述可靠相似变化像元的数目大于或等于预设数目时,对所述可靠带粗差像元依次进行分类、解混和残差计算,得到中间融合影像。
在一种实施方式中,所述对所述中间融合影像进行空间滤波和残差分配,得到目标融合影像的步骤包括:
遍历所述中间融合影像,以所述中间融合影像中的每个像元为第二中心像元建立移动窗口;
对所述移动窗口图像中的光谱相似像元进行滤波,得到滤波融合影像;
对所述滤波融合影像进行残差分配,得到目标融合影像。
在一种实施方式中,所述移动窗口的尺寸w=round[wmax-mean(RI)×(wmax-wmin)],其中,所述round为四舍五入取整函数,所述wmax为预设尺寸最大值,所述wmin为预设尺寸最小值,所述RI为所述光谱相似像元的像元可靠性系数。
第二方面,本发明实施例还提供一种智能终端,包括有存储器,以及一个或者一个以上的程序,其中一个或者一个以上程序存储于存储器中,且经配置以由一个或者一个以上处理器执行所述一个或者一个以上程序包含用于执行如上述任意一项所述的遥感影像时空融合方法。
第三方面,本发明实施例还提供一种非临时性计算机可读存储介质,当所述存储介质中的指令由电子设备的处理器执行时,使得电子设备能够执行如上述中任意一项所述的遥感影像时空融合方法。
本发明的有益效果:本发明通过获取第一时相的高分辨率影像和第一低分辨率影像,以及第二时相的第二低分辨率影像,根据高分辨率影像、第一低分辨率影像和第二低分辨率影像计算像元可靠性系数,利用模糊C均值聚类法对高分辨率影像中的像元进行分类,并根据像元可靠性系数筛选出各类像元各类像元参与全局解混,以得到高分辨率影像、第一低分辨率影像和第二低分辨率影像的初始融合影像,基于自适应局部解混模型修正初始融合影像中的粗差,得到中间融合影像,对中间融合影像进行空间滤波和残差分配,得到目标融合影像,减弱了异源影像差异对融合的影响,有效地提升了融合算法地可靠性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明中记载的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明实施例提供第一流程示意图;
图2为本发明实施例提供第二流程示意图;
图3为本发明实施例提供第三流程示意图;
图4为可靠相似变化像元求算过程示意图;
图5为局部窗口分类解混过程示意图;
图6为本发明实施例提供的智能终端的内部结构原理框图。
具体实施方式
本发明公开了一种遥感影像时空融合方法,智能终端及计算机可读存储介质,为使本发明的目的、技术方案及效果更加清楚、明确,以下参照附图并举实施例对本发明进一步详细说明。应当理解,此处所描述的具体实施例仅用以解释本发明,并不用于限定本发明。
本技术领域技术人员可以理解,除非特意声明,这里使用的单数形式“一”、“一个”、“所述”和“该”也可包括复数形式。应该进一步理解的是,本发明的说明书中使用的措辞“包括”是指存在所述特征、整数、步骤、操作、元件和/或组件,但是并不排除存在或添加一个或多个其他特征、整数、步骤、操作、元件、组件和/或它们的组。应该理解,当我们称元件被“连接”或“耦接”到另一元件时,它可以直接连接或耦接到其他元件,或者也可以存在中间元件。此外,这里使用的“连接”或“耦接”可以包括无线连接或无线耦接。这里使用的措辞“和/或”包括一个或更多个相关联的列出项的全部或任一单元和全部组合。
本技术领域技术人员可以理解,除非另外定义,这里使用的所有术语(包括技术术语和科学术语),具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是,诸如通用字典中定义的那些术语,应该被理解为具有与现有技术的上下文中的意义一致的意义,并且除非像这里一样被特定定义,否则不会用理想化或过于正式的含义来解释。
本实施例提供遥感影像时空融合方法,该方法可以应用于(数据安全保护的)智能终端。具体如图1所示,所述方法包括:
步骤S10,获取第一时相的高分辨率影像和第一低分辨率影像,以及第二时相的第二低分辨率影像;
通过遥感卫星获取不同时相的遥感影像,具体地,在时相为T1时获取目标地点的高分辨率影像和低分辨率影像,在时相为T2时获取目标地点的低分辨率影像,由于遥感卫星获取影像的周期性的原因,在短时间内无法获取在T2时相的高分辨率影像,本实施例的目的就是为了预测在T2时相的高分辨率影像。为了将两幅低分辨率影像进行区分,将在T1时相获取的低分辨率影像确定为第一低分辨率影像,将在T2时相获取的低分辨率影像确定为第二低分辨率影像。本实施例将T1时相确定为第一时相,将T2时相确定为第二时相,通过融合高分辨率影像、第一低分辨率影像以及第二低分辨率影像进行融合,预测出在第二时相的高分辨率应影像。
在一些实施例中,获取到高分辨率影像、第一低分辨率应影像和第二低分辨率影像之后,分别对高分辨率影像、第一低分辨率影像和第二分辨率影像进行预处理,预处理的方法包括配准和辐射矫正。
步骤S20,根据所述高分辨率影像、第一低分辨率影像和第二低分辨率影像计算像元可靠性系数;
本实施例用各个影像中每个像元的光谱值,即像元光谱值组成的集合表示影像,第一时相获取的高分辨率影像用F1表示,F1为高分辨率影像中各个像元对应的像元光谱值组成的集合,第一时相获取的第一低分辨率影像用C1表示,C1为第一低分辨率影像中各个像元对应的像元光谱值组成的集合,第二时相获取的第二低分辨率影像用C2表示,C2为第二低分辨率影像中各个像元对应的像元光谱值组成的集合。
关于像元可靠性系数的计算,需要获取高分辨率影像、第一低分辨率影像和第二低分辨率影像中的像元信息。
参照图2,在一些具体的实施例中,步骤S20包括:
步骤S21,对所述高分辨率影像中各像元对应的像元光谱值进行降采样处理,得到与所述第一低分辨率影像的分辨率相同的第三低分辨率影像;
步骤S22,对所述第三低分辨率影像中各像元对应的像元光谱值和第一低分辨率影像中各像元对应的像元光谱值进行线性回归分析,得到回归系数;
步骤S23,将所述回归系数代入系统误差矫正函数中以矫正第一低分辨率影像的系统误差,得到第四低分辨率影像;
步骤S24,根据所述第一低分辨率影像、所述第二低分辨率影像、所述第三低分辨率影像和所述第四低分辨率影像计算像元可靠性系数。
关于像元光谱值的计算,首先对高分辨率影像进行降采样处理,具体的,是对高分辨率影像中的各个像元对应的像元光谱值进行降采样处理,得到低分辨影像,即第三低分辨率影像,第三低分辨率影像的分辨率与第一低分辨率影像的分辨率相同,用φdown(F1)表示第三低分辨率影像,φdown为降采样函数,可以理解φdown(F1)为第三低分辨率影像中各个像元对应的像元光谱值组成的集合,需要说明的是,第二低分辨率影像的分辨率与第一低分辨率影像的也是相同的。
获得第三低分辨率影像是为了对第一低分辨率影像做系统矫正,具体的,对φdown(F1)和C1作线性回归分析,即φdown(F1)=J×C1+K+ξ,得到回归系数J和K,其中,ξ为线性回归后的残差,根据回归系数矫正第一低分辨率影像的系统误差,系统误差的矫正公式为φsystem(C1)=J×C1+K,φsystem(C1)表示系统误差矫正后的低分辨率影像,即第四低分辨率影像,φsystem(C1)为第四低分辨率影像中各个像元对应的像元光谱值组成的集合。需要说明的是,对第一低分辨率影像进行系统矫正,不改变影像的空间分辨率。
高分辨率影像、第一低分辨率影像和第二低分辨率影像中的像元信息包括每幅影像的各个像元对应的像元光谱值,根据高分辨率影像、第一低分辨率影像的像元信息获得第三低分辨率影像和第四低分辨率影像,根据第三低分辨率影像与第四低分辨率影像计算像元可靠性系数。
在一些具体的实施例中,步骤S24还包括:
步骤a,将所述第一低分辨率影像中各像元对应的像元光谱值、所述第二低分辨率影像中各像元对应的像元光谱值、所述第三低分辨率影像中各像元对应的像元光谱值和所述第四低分辨率影像中各像元对应的像元光谱值代入第一公式中,得到阈值系数,其中,所述第一公式为:所述Q(b)为所述阈值系数,所述b为波段号,所述ΔC为所述第一低分辨率影像中各像元与第二低分辨率影像中各像元的像元光谱值的差值组成的集合,所述φdown(F1)为第三低分辨率影像中各像元对应的像元光谱值组成的集合,所述φsystem(C1)为第四低分辨率影像中各像元对应的像元光谱值组成的集合,所述mean表示求平均值函数;
步骤b,将所述阈值系数代入第二公式,得到像元可靠性系数,其中,所述第二公式为:其中,所述RI(xi,yi,b)表示坐标为(xi,yi)的像元对应的像元可靠性系数,所述φdown[F1(xi,yi,b)]表示所述第三低分辨率影像中坐标为(xi,yi)的像元对应的像元光谱值,所述φsystem[C1(xi,yi,b)]为第四低分辨率影像中坐标为(xi,yi)的像元对应的像元光谱值。
φdown(F1(b))-φsystem(C1(b))表示异源影像间异源影像配准误差造成的差异,观测时间的差异等的综合作用定量表示,一般情况下φdown(F1(b))-φsystem(C1(b))符合高斯分布,ΔC为所述第一低分辨率影像的各像元与第二低分辨率影像对应的各像元的像元光谱值的差值组成的集合,用于表示低分辨率影像从T1时相到T2时相的变化值,在考虑时相变化值基础上,计算阈值系数Q(b),为保证像元可靠性系数取值范围合理,当Q(b)<2时,令Q(b)=2。进一步地计算像元可靠性系数,每一个像元对应一个像元可靠性系数,(xi,yi)为像元的坐标。为避免可靠性系数在后续计算过程中出现错误,当RI(xi,yi,b)<0.1时,令RI(xi,yi,b)=0.1。
步骤S30,利用模糊C均值聚类法对所述高分辨率影像中的像元进行分类,并根据像元可靠性系数筛选出各类像元参与全局解混,以得到所述高分辨率影像、所述第一低分辨率影像和所述第二低分辨率影像的初始融合影像;
本实施例提供一种融合方法,首先对高分辨率影像设置不同的分类数进行分类,得到不同分类数下的分类结果,进而通过全局解混模型得到每次分类后高分辨率影像、第一低分辨率影像和第二低分辨率影像的融合影像,得到的融合影像的幅数等于分类的次数,从多幅融合影像中选出最佳的一幅作为初步的融合影像,即初始融合影像。
在一些具体的实施例中,步骤S30包括:
步骤c,设置最大分类值nmax和起始分类值n;
步骤d,当分类数为n时,利用模糊C均值聚类法将所述高分辨率影像中的像元分成n类,根据所述像元可靠性系数从各类像元中确定用于全局解混的对象像元;
步骤e,将所述对象像元的丰度输入全局解混模型,得到对象像元在所述第一低分辨率影像和所述第二低分辨率影像之间的两时相平均变化值;
步骤f,将所述两时相平均变化值加入所述高分辨率影像中,得到融合影像;
步骤g,将n以1为单元逐级递加,直至n=nmax,得到每一级n对应的融合影像;
步骤h,计算每一级n对应的融合影像残差值,将所述残差值的最小值对应的融合影像确定为初步融合影像。
实现对高分辨率影像的分类,首先设置最大的分类数,即最大分类值nmax,以及最小的分类数,即起始分类值n。利用模糊C均值聚类法对高分辨率影像中的像元进行分类。首先令分类数为起始分类值n,利用模糊C均值聚类法将高分辨率影像中的像元分成n类,计算每个像元在其所在类的丰度,以及判断每类中的像元是否为非变化像元。
判断像元是否为非变化像元的方法为首先计算变化阈值的上限值和下限值 当坐标为(xi,yi)的像元在第一低分辨率影像和在第二低分辨率影像中的像元光谱值的差值ΔC(xi,yi),当且确定坐标为(xi,yi)的像元为变化像元;当确定坐标为(xi,yi)的像元为非变化像元,将每一类像元中的非变化像元的像元可靠性系数和丰度分别按照从大到小的顺序排列,得到像元可靠性系数序列和丰度序列,选取位于像元可靠性系数序列的前30%的像元可靠性系数对应的非变化像元,当选取的非变化像元数目不足预设的数目m时,继续依据丰度序列中丰度值从大到小的顺序选取非变化像元,直至选取的非变化像元的数目达到m。将m个非变化像元确定为用于全局解混的对象像元。将每类中的对象像元的丰度输入全局解混模型,得到对象像元在第一低分辨率影像和第二低分辨率影像之间各类的两时相平均变化值,全局解混模型为:
解混结果ΔF为各类的两时相平均变化值,通过约束最小二乘法可以解算得到,将两时相平均变化值加到高分辨率影像中,得到分类数为起始分类数n时高分辨率影像、第一低分辨率影像和第二低分辨率影像的融合影像。
令分类数n=n+1,使得n以1为单位逐级递加,根据上述同样的方法得到n=n+1时高分辨率影像、第一低分辨率影像和第二低分辨率影像的融合影像,直至n=nmax,得到最后一个融合影像。得到每一级n对应的融合影像后,计算每一幅影像的残差值,将融合影像输入残差模型中,式中,Rn表示残差模型计算结果,Pn表示与n对应的融合影像,I表示一幅低分辨率影像中低分辨率像元的总数,B表示总波段数,η表示融合影像中坐标为(xi,yi)的像元是否检测为变化,若是,则η=0;若不是,则η=1,比较每个融合影像的Rn,取最小的Rn对应的融合影像为初始融合影像,记为
自适应全局解混模型能够自适应选择高可靠性且未发生类别变化的像元参与全局解混计算,能够自适应选择最佳分类数下的解混结果,有效提升融合算法在复原光谱变化的能力。
步骤S40,基于自适应局部解混模型修正所述初始融合影像中的粗差,得到中间融合影像;
本实施例通过对初始融合影像进行分类和解混修正初始融合影像中的粗差,用于解混的模型为自适应局部解混模型。
参照图3,在一些具体的实施例中,步骤S40还包括:
步骤S41,对所述初始融合影像进行降采样处理,得到低分辨率融合影像;
步骤S42,根据变化检测技术确定所述低分辨率融合影像中的带粗差像元,并根据所述带粗差像元的平均可靠性系数确定可靠带粗差像元;
步骤S43,以所述可靠带粗差像元为第一中心像元建立预设尺寸的局部窗口;
步骤S44,计算局部窗口中所述可靠带粗差像元与所述第一中心像元的时相变化值,并根据所述时相变化值从所述可靠带粗差像元中选择可靠相似变化像元;
步骤S45,当所述可靠相似变化像元的数目大于或等于预设数目时,利用局部解混模型对所述可靠带粗差像元依次进行分类、解混和残差计算,得到中间融合影像。
初始融合影像是高分辨率的,对初始融合影像进行降采样处理,将初始融合影像转化成低分辨率的,得到低分辨率融合影像。利用高斯模型检测低分辨率融合影像是否存在带粗差像元,具体的,设置阈值的上限值和下限值 其中,当低分辨率融合影像中的像元的CFd大于或者小于确定该像元为带粗差像元。当带粗差像元的平均可靠性系数大于0.8时,确定该带粗差像元为可靠带粗差像元,平均可靠性系数是指带粗差像元在各波段中像元可靠性系数的平均值。利用局部窗口的思想确定可靠带粗差像元中的可靠相似变化像元,设置局部窗口尺寸的大小,即预设尺寸,一般的预设尺寸的大小为5*5低分辨率像元,以可靠带粗差像元为中心像元,即第一中心像元建立预设尺寸大小的局部窗口,在局部窗口内,可靠带粗差像元与第一中心像元有相似的时相变化值的视为可靠相似变化像元,由于第一中心像元是可靠带粗差像元,所以第一中心像元也是可靠相似变化像元。图4为可靠相似变化像元求算过程示意图。
设局部窗口内除第一中心像元以外可靠带粗差像元的坐标为(xw,yw),第一中心像元的坐标为(x0,y0),通过下述公式计算可靠带粗差像元与第一中心像元的时相变化值的相似程度,并设置阈值QDT,QDT=stddev(DT)/2,若DT(xw,yw)<QDT,确定坐标为(xw,yw)的可靠带粗差像元为可靠相似变化像元,统计每个局部窗口内可靠相似变化像元的数目,若局部窗口内可靠相似变化像元数目大于或等于9,则对局部窗口内的可靠相似变化像元进行分类和解混;若局部窗口内可靠相似变化像元数目小于9,则在局部窗内在建立3*3低分辨率像元大小的局部窗口,对3*3局部窗口内所有像元进行分类和解混。无论是对可靠相似变化像元还是对3*3局部窗口内的所有像元进行分类和解混,步骤与上述的分类解混步骤相同,同样的设置最大分类数nmax,以及起始分类数n,利用模糊C均值分类法对待分解及解混的像元进行分类,并获取待分解及解混的像元的丰度,将其代入局部解混模型中:
式中,ΔF是各类像元的两时相平均变化值,通过最小二乘法解算得到,r为用于分类和解混的像元的数目,将各类像元的两时相平均变化值加到初步融合影像中的局部窗口,每个局部窗口均具有不同级n对应的局部融合结果,可以理解,每个局部窗口均有多个局部融合结果,将每个局部窗口的每个局部融合结果代入残差模型:将具有最小残差模型计算结果的局部融合结果,即为最优的局部融合结果,作为局部窗口的影像LP,对该影像进行可靠性评估:
ra能够描述利用局部解混对初步融合结果中的粗差进行修正后的可靠性和修正率,当ra(x0,y0)<ARI(x0,y0),不对最优的局部融合结果做处理,否则利用阈值法和加权函数消除或减弱局部解混造成的错误,LP'为消除或减弱局部解混造成的错误,并影像LP'作为第一中心像元对应的局部窗口的最终的最优局部融合结果。将上述步骤遍历所有可靠带粗差像元,当每个局部窗口均获得最优的局部融合结果,得到中间融合影像
具体的,在局部窗口内可靠相似变化像元数目大于或等于9和局部窗口内可靠相似变化像元数目小于9情况下分类和解混的过程参见图5。
自适应局部解混模型能够对自适应全局解混存在粗差的像元进行再次解混,并通过加权公式修正其光谱值,有效提升融合算法在复原光谱剧烈变化并保留空间细节信息的能力。
步骤S50,对所述中间融合影像进行空间滤波和残差分配,得到目标融合影像。
在一些具体的实施例中,遍历所述中间融合影像,以所述中间融合影像中的每个像元为第二中心像元建立移动窗口;
步骤i,对所述移动窗口图像中的光谱相似像元进行滤波,得到滤波融合影像;
步骤j,对所述滤波融合影像进行残差分配,得到目标融合影像。
用移动窗口中遍历中间融合影像的每个高分辨率像元作为中心像元,,即第二中心像元,在移动窗口中根据非第二中心像元与第二中心像元的光谱距离选择至多K个光谱相似像元。移动窗口的大小w(移动窗口的长宽为2×w+1)和光谱相似像元的最大值K由所有像元所有波段的平均时空融合可靠性系数决定:w=K=round[wmax-mean(RI)×(wmax-wmin)],其中round为四舍五入取整函数,wmax和wmin为预设的最大和最小窗口大小,分别设置为60和20。
为了消除局部解混造成的块效应,根据光谱相似像元与第二中心像元的欧式距离确定权值,通过加权结合每个光谱相似像元的时相变化值作为新的第二中心像元时相变化值,对中间融合影像进行空间滤波:
本实施例通过获取第一时相的高分辨率影像和第一低分辨率影像,以及第二时相的第二低分辨率影像,根据所述高分辨率影像、第一低分辨率影像和第二低分辨率影像计算像元可靠性系数;利用模糊C均值聚类法对高分辨率影像中的像元进行分类,并根据像元可靠性系数筛选出各类像元参与全局解混,以得到高分辨率影像、第一低分辨率影像和第二低分辨率影像的初始融合影像;再基于自适应局部解混模型修正所述初始融合影像中的粗差,得到中间融合影像;对所述中间融合影像进行空间滤波和残差分配,得到目标融合影像,减弱了异源影像差异对融合的影响,有效地提升了融合算法地可靠性。
为了说明本实施例提供的遥感影像时空数据融合方法的融合效果,将本实施例提出的时空融合方法与时空自适应反射融合模型,灵活时空数据融合方法和面向对象的时空融合方法进行对比实验。用于实验对比的两期影像产生的变化的主要原因是物候生长。
表1给出了上述四种时空融合方法的定量分析指标。其中方法A为时空自适应反射融合模型,方法B为灵活时空数据融合方法,方法C为面向对象的时空融合方法,方法D为本实施例提供的遥感影像时空融合方法;选用均方根误差(RMSE)、结构相似性(SSIM)和相关系数(r)三个精度指标,分别反映出融合结果与真实影像的光谱差异、结构相似性和相关性,前者越接近0、后两者越接近1则说明融合效果越好。
从表1的对比中可得,本实施例提供的时空融合结果明显优于其他三种时空融合结果:与时空自适应反射融合模型(方法A)和面向对象的时空融合方法(方法C)相比,本实施例在所有波段拥有最小的均方根误差(RMSE),最大的结构相似性(SSIM)和相关系数(r),与灵活时空数据融合方法(方法B)相比,本实施例在除波段6外的5个波段拥有最小的均方根误差(RMSE),最大的结构相似性(SSIM)和相关系数(r)。因此,本实施例提供的遥感影像时空融合方法能够取得较优的时空融合结果。
表1
注:加粗数值表示实验最佳值。
基于上述实施例,本发明还提供了一种智能终端,其原理框图可以如图5所示。该智能终端包括通过系统总线连接的处理器、存储器、网络接口、显示屏、温度传感器。其中,该智能终端的处理器用于提供计算和控制能力。该智能终端的存储器包括非易失性存储介质、内存储器。该非易失性存储介质存储有操作系统和计算机程序。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该智能终端的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种遥感影像时空融合方法。该智能终端的显示屏可以是液晶显示屏或者电子墨水显示屏,该智能终端的温度传感器是预先在智能终端内部设置,用于检测内部设备的运行温度。
本领域技术人员可以理解,图6中的原理图,仅仅是与本发明方案相关的部分结构的框图,并不构成对本发明方案所应用于其上的智能终端的限定,具体的智能终端可以包括比图中所示更多或更少的部件,或者组合某些部件,或者具有不同的部件布置。
在一个实施例中,提供了一种智能终端,包括有存储器,以及一个或者一个以上的程序,其中一个或者一个以上程序存储于存储器中,且经配置以由一个或者一个以上处理器执行所述一个或者一个以上程序包含用于进行以下操作的指令:
获取第一时相的高分辨率影像和第一低分辨率影像,以及第二时相的第二低分辨率影像;
根据所述高分辨率影像、第一低分辨率影像和第二低分辨率影像计算像元可靠性系数;
利用模糊C均值聚类法对所述高分辨率影像中的像元进行分类,并根据像元可靠性系数筛选出各类像元参与全局解混,以得到所述高分辨率影像、所述第一低分辨率影像和所述第二低分辨率影像的初始融合影像;
基于自适应局部解混模型修正所述初始融合影像中的粗差,得到中间融合影像;
对所述中间融合影像进行空间滤波和残差分配,得到目标融合影像。
本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程,是可以通过计算机程序来指令相关的硬件来完成,所述的计算机程序可存储于一非易失性计算机可读取存储介质中,该计算机程序在执行时,可包括如上述各方法的实施例的流程。其中,本发明所提供的各实施例中所使用的对存储器、存储、数据库或其它介质的任何引用,均可包括非易失性和/或易失性存储器。非易失性存储器可包括只读存储器(ROM)、可编程ROM(PROM)、电可编程ROM(EPROM)、电可擦除可编程ROM(EEPROM)或闪存。易失性存储器可包括随机存取存储器(RAM)或者外部高速缓冲存储器。作为说明而非局限,RAM以多种形式可得,诸如静态RAM(SRAM)、动态RAM(DRAM)、同步DRAM(SDRAM)、双数据率SDRAM(DDRSDRAM)、增强型SDRAM(ESDRAM)、同步链路(Synchlink)DRAM(SLDRAM)、存储器总线(Rambus)直接RAM(RDRAM)、直接存储器总线动态RAM(DRDRAM)、以及存储器总线动态RAM(RDRAM)等。
综上所述,本发明公开了一种遥感影像时空融合方法、智能终端、存储介质,所述方法包括:
获取第一时相的高分辨率影像和第一低分辨率影像,以及第二时相的第二低分辨率影像;
根据所述高分辨率影像、第一低分辨率影像和第二低分辨率影像计算像元可靠性系数;
利用模糊C均值聚类法对所述高分辨率影像中的像元进行分类,并根据像元可靠性系数筛选出各类像元参与全局解混,以得到所述高分辨率影像、所述第一低分辨率影像和所述第二低分辨率影像的初始融合影像;
基于自适应局部解混模型修正所述初始融合影像中的粗差,得到中间融合影像;
对所述中间融合影像进行空间滤波和残差分配,得到目标融合影像。
基于上述实施例,本发明公开了一种遥感影像时空融合方法,应当理解的是,本发明的应用不限于上述的举例,对本领域普通技术人员来说,可以根据上述说明加以改进或变换,所有这些改进和变换都应属于本发明所附权利要求的保护范围。
Claims (9)
1.一种遥感影像时空融合方法,其特征在于,所述方法包括:
获取第一时相的高分辨率影像和第一低分辨率影像,以及第二时相的第二低分辨率影像;
根据所述高分辨率影像、第一低分辨率影像和第二低分辨率影像计算像元可靠性系数;
利用模糊C均值聚类法对所述高分辨率影像中的像元进行分类,并根据像元可靠性系数筛选出各类像元参与全局解混,以得到所述高分辨率影像、所述第一低分辨率影像和所述第二低分辨率影像的初始融合影像;
基于自适应局部解混模型修正所述初始融合影像中的粗差,得到中间融合影像;
对所述中间融合影像进行空间滤波和残差分配,得到目标融合影像;
所述根据所述高分辨率影像、第一低分辨率影像和第二分辨率影像计算像元可靠性系数的步骤包括:
对所述高分辨率影像中各像元对应的像元光谱值进行降采样处理,得到与所述第一低分辨率影像的分辨率相同的第三低分辨率影像;
对所述第三低分辨率影像中各像元对应的像元光谱值和第一低分辨率影像中各像元对应的像元光谱值进行线性回归分析,得到回归系数;
将所述回归系数代入系统误差矫正函数中以矫正第一低分辨率影像的系统误差,得到第四低分辨率影像;
根据所述第一低分辨率影像、所述第二低分辨率影像、所述第三低分辨率影像和所述第四低分辨率影像计算像元可靠性系数集。
2.根据权利要求1所述的遥感影像时空融合方法,其特征在于,所述根据所述高分辨率影像、第一低分辨率影像和第二分辨率影像中像元光谱值计算像元可靠性系数集的步骤之前,分别对所述高分辨率影像、第一低分辨率影像和第二分辨率影像进行预处理,预处理的方法包括配准和辐射矫正。
3.根据权利要求1所述的遥感影像时空融合方法,其特征在于,所述根据所述第一低分辨率影像、所述第二低分辨率影像、所述第三低分辨率影像和所述第四低分辨率影像计算像元可靠性系数的步骤包括:
将所述第一低分辨率影像中各像元对应的像元光谱值、所述第二低分辨率影像中各像元对应的像元光谱值、所述第三低分辨率影像中各像元对应的像元光谱值和所述第四低分辨率影像中各像元对应的像元光谱值代入第一公式中,得到阈值系数,其中,所述第一公式为:所述Q(b)为所述阈值系数,所述b为波段号,所述ΔC为所述第一低分辨率影像中各像元与第二低分辨率影像中各像元的像元光谱值的差值组成的集合,所述φdown(F1)为第三低分辨率影像中各像元对应的像元光谱值组成的集合,所述φsystem(C1)为第四低分辨率影像中各像元对应的像元光谱值组成的集合,所述mean表示求平均值函数;
将所述阈值系数代入第二公式,得到像元可靠性系数,其中,所述
4.根据权利要求3所述的遥感影像时空融合方法,其特征在于,所述利用模糊C均值聚类法对所述高分辨率影像中的像元进行分类,并根据像元可靠性系数筛选出各类像元参与全局解混,以得到所述高分辨率影像、所述第一低分辨率影像和所述第二低分辨率影像的初始融合影像的步骤包括:
设置最大分类值nmax和起始分类值n;
当分类数为n时,利用模糊C均值聚类法将所述高分辨率影像中的像元分成n类,根据所述像元可靠性系数从各类像元中确定用于全局解混的对象像元;
将所述对象像元的丰度输入全局解混模型,得到对象像元在所述第一低分辨率影像和所述第二低分辨率影像之间的两时相平均变化值;
将所述两时相平均变化值加入所述高分辨率影像中,得到融合影像;
将n以1为单元逐级递加,直至n=nmax,得到每一级n对应的融合影像;
计算每一级n对应的融合影像残差值,将所述残差值的最小值对应的融合影像确定为初始融合影像。
5.根据权利要求4所述的遥感影像时空融合方法,其特征在于,所述基于自适应局部解混模型修正所述初始融合影像中的粗差,得到中间融合影像的步骤包括:
对所述初始融合影像进行降采样处理,得到低分辨率融合影像;
根据变化检测技术确定所述低分辨率融合影像中的带粗差像元,并根据所述带粗差像元的平均可靠性系数确定可靠带粗差像元;
以所述可靠带粗差像元为第一中心像元建立预设尺寸的局部窗口;
计算局部窗口中所述可靠带粗差像元与所述第一中心像元的时相变化值,并根据所述时相变化值从所述可靠带粗差像元中选择可靠相似变化像元;
当所述可靠相似变化像元的数目大于或等于预设数目时,对所述可靠带粗差像元依次进行分类、解混和残差计算,得到中间融合影像。
6.根据权利要求1所述的遥感影像时空融合方法,其特征在于,所述对所述中间融合影像进行空间滤波和残差分配,得到目标融合影像的步骤包括:
遍历所述中间融合影像,以所述中间融合影像中的每个像元为第二中心像元建立移动窗口;
对所述移动窗口图像中的光谱相似像元进行滤波,得到滤波融合影像;
对所述滤波融合影像进行残差分配,得到目标融合影像。
7.根据权利要求6所述的遥感影像时空融合方法,其特征在于,所述移动窗口的尺寸w=round[wmax-mean(RI)×(wmax-wmin)],其中,所述round为四舍五入取整函数,所述wmax为预设尺寸最大值,所述wmin为预设尺寸最小值,所述RI为所述光谱相似像元的像元可靠性系数。
8.一种智能终端,其特征在于,包括有存储器,以及一个或者一个以上的程序,其中一个或者一个以上程序存储于存储器中,且经配置以由一个或者一个以上处理器执行所述一个或者一个以上程序包含用于执行如权利要求1-7中任意一项所述的方法。
9.一种非临时性计算机可读存储介质,其特征在于,当所述存储介质中的指令由电子设备的处理器执行时,使得电子设备能够执行如权利要求1-7中任意一项所述的方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110507192.5A CN113327197B (zh) | 2021-05-10 | 2021-05-10 | 遥感影像时空融合方法、智能终端及计算机可读存储介质 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110507192.5A CN113327197B (zh) | 2021-05-10 | 2021-05-10 | 遥感影像时空融合方法、智能终端及计算机可读存储介质 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113327197A CN113327197A (zh) | 2021-08-31 |
CN113327197B true CN113327197B (zh) | 2023-01-24 |
Family
ID=77415180
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110507192.5A Active CN113327197B (zh) | 2021-05-10 | 2021-05-10 | 遥感影像时空融合方法、智能终端及计算机可读存储介质 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113327197B (zh) |
Families Citing this family (1)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN117115679B (zh) * | 2023-10-25 | 2024-06-25 | 北京佳格天地科技有限公司 | 一种时空融合遥感影像对筛选方法 |
Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN110363246A (zh) * | 2019-07-18 | 2019-10-22 | 滨州学院 | 一种高时空分辨率植被指数ndvi的融合方法 |
CN111667437A (zh) * | 2019-03-06 | 2020-09-15 | 中国科学院微电子研究所 | 基于傅里叶域解混的遥感图像空时融合方法 |
CN111832518A (zh) * | 2020-07-22 | 2020-10-27 | 桂林电子科技大学 | 基于时空融合的tsa遥感影像土地利用方法 |
Family Cites Families (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN103177431B (zh) * | 2012-12-26 | 2015-10-14 | 中国科学院遥感与数字地球研究所 | 一种多源遥感数据时空融合方法 |
CN105046648B (zh) * | 2015-06-25 | 2019-01-22 | 北京师范大学 | 一种构建高时空遥感数据的方法 |
CN108613933A (zh) * | 2018-06-13 | 2018-10-02 | 中南林业科技大学 | 基于多源遥感数据融合的林地干旱时空动态监测方法 |
CN111666896A (zh) * | 2020-06-09 | 2020-09-15 | 中国科学院地理科学与资源研究所 | 一种基于线性融合模型的遥感影像时空融合方法 |
CN112017135B (zh) * | 2020-07-13 | 2021-09-21 | 香港理工大学深圳研究院 | 一种遥感影像数据时空融合的方法、系统及设备 |
CN112508832B (zh) * | 2020-12-03 | 2024-02-13 | 中国矿业大学 | 一种面向对象的遥感影像数据时空融合方法、系统及设备 |
-
2021
- 2021-05-10 CN CN202110507192.5A patent/CN113327197B/zh active Active
Patent Citations (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN111667437A (zh) * | 2019-03-06 | 2020-09-15 | 中国科学院微电子研究所 | 基于傅里叶域解混的遥感图像空时融合方法 |
CN110363246A (zh) * | 2019-07-18 | 2019-10-22 | 滨州学院 | 一种高时空分辨率植被指数ndvi的融合方法 |
CN111832518A (zh) * | 2020-07-22 | 2020-10-27 | 桂林电子科技大学 | 基于时空融合的tsa遥感影像土地利用方法 |
Also Published As
Publication number | Publication date |
---|---|
CN113327197A (zh) | 2021-08-31 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN111723860B (zh) | 一种目标检测方法及装置 | |
US11487966B2 (en) | Image processing method and apparatus for target recognition | |
CN109086656A (zh) | 机场异物检测方法、装置、计算机设备及存储介质 | |
CN112419155B (zh) | 一种全极化合成孔径雷达影像超分辨率重建方法 | |
CN112017135A (zh) | 一种遥感影像数据时空融合的方法、系统及设备 | |
CN111242026B (zh) | 一种基于空间层次感知模块和度量学习的遥感图像目标检测方法 | |
CN111524063A (zh) | 一种遥感图像融合方法及装置 | |
CN112990314B (zh) | 基于改进孤立森林算法的高光谱图像异常检测方法及装置 | |
CN113327197B (zh) | 遥感影像时空融合方法、智能终端及计算机可读存储介质 | |
CN111563896B (zh) | 一种用于接触网异常检测的图像处理方法 | |
CN112419202A (zh) | 基于大数据及深度学习的野生动物图像自动识别系统 | |
CN113674191A (zh) | 一种基于条件对抗网络的弱光图像增强方法和装置 | |
CN109903272A (zh) | 目标检测方法、装置、设备、计算机设备和存储介质 | |
CN116883763B (zh) | 一种基于深度学习的汽车零部件缺陷检测方法及系统 | |
CN107194351B (zh) | 基于韦伯局部对称图结构的人脸识别特征提取方法 | |
CN113537020B (zh) | 基于改进神经网络的复数sar图像目标识别方法 | |
CN116740728B (zh) | 一种用于晶圆读码器动态获取方法和系统 | |
CN113487607A (zh) | 基于多视场图像的缺陷检测方法及装置 | |
CN113066030B (zh) | 一种基于空谱融合网络的多光谱图像全色锐化方法及系统 | |
Sadeghi et al. | A new fuzzy measurement approach for automatic change detection using remotely sensed images | |
CN117710728A (zh) | Sar图像目标识别方法、装置、计算机设备和存储介质 | |
CN112419270A (zh) | 元学习下的无参考图像质量评价方法、装置及计算机设备 | |
CN109410259B (zh) | 基于置信度的结构化的双目深度图上采样方法 | |
CN110648351A (zh) | 基于稀疏表示的多外观模型融合目标跟踪方法及其装置 | |
CN114743067A (zh) | 训练数据的增强方法、装置、计算机设备及存储介质 |
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 |