CN112950492B - 基于自适应各向异性扩散的全极化sar图像去噪方法 - Google Patents
基于自适应各向异性扩散的全极化sar图像去噪方法 Download PDFInfo
- Publication number
- CN112950492B CN112950492B CN202110117113.XA CN202110117113A CN112950492B CN 112950492 B CN112950492 B CN 112950492B CN 202110117113 A CN202110117113 A CN 202110117113A CN 112950492 B CN112950492 B CN 112950492B
- Authority
- CN
- China
- Prior art keywords
- scattering
- fully
- adaptive
- gradient
- self
- 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
- 238000000034 method Methods 0.000 title claims abstract description 62
- 238000009792 diffusion process Methods 0.000 title claims abstract description 31
- 230000010287 polarization Effects 0.000 claims abstract description 91
- 238000001914 filtration Methods 0.000 claims abstract description 79
- 239000011159 matrix material Substances 0.000 claims abstract description 34
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 12
- 230000007246 mechanism Effects 0.000 claims abstract description 11
- 239000003550 marker Substances 0.000 claims abstract description 8
- 230000003044 adaptive effect Effects 0.000 claims description 63
- 238000004364 calculation method Methods 0.000 claims description 5
- 238000003064 k means clustering Methods 0.000 claims description 4
- 238000006243 chemical reaction Methods 0.000 claims description 2
- 238000010276 construction Methods 0.000 claims description 2
- 230000006978 adaptation Effects 0.000 claims 1
- 230000002401 inhibitory effect Effects 0.000 abstract description 3
- 238000010586 diagram Methods 0.000 description 21
- 230000001629 suppression Effects 0.000 description 12
- 230000001427 coherent effect Effects 0.000 description 7
- 230000005855 radiation Effects 0.000 description 7
- 239000013535 sea water Substances 0.000 description 7
- 238000011158 quantitative evaluation Methods 0.000 description 6
- 230000008859 change Effects 0.000 description 5
- 238000005388 cross polarization Methods 0.000 description 5
- 238000011156 evaluation Methods 0.000 description 5
- 230000005540 biological transmission Effects 0.000 description 4
- 230000014759 maintenance of location Effects 0.000 description 4
- 238000012545 processing Methods 0.000 description 4
- 230000000007 visual effect Effects 0.000 description 4
- 238000001514 detection method Methods 0.000 description 3
- 238000011161 development Methods 0.000 description 3
- 230000018109 developmental process Effects 0.000 description 3
- 238000000605 extraction Methods 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 3
- 230000008569 process Effects 0.000 description 3
- 238000011160 research Methods 0.000 description 3
- 239000000126 substance Substances 0.000 description 3
- 238000013459 approach Methods 0.000 description 2
- 150000001875 compounds Chemical class 0.000 description 2
- 230000021615 conjugation Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000002592 echocardiography Methods 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 230000035515 penetration Effects 0.000 description 2
- 230000000704 physical effect Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 229910052500 inorganic mineral Inorganic materials 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 239000011707 mineral Substances 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 239000003305 oil spill Substances 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 238000004451 qualitative analysis Methods 0.000 description 1
- 230000011218 segmentation Effects 0.000 description 1
- 238000010187 selection method Methods 0.000 description 1
- 239000002689 soil Substances 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/23—Clustering techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F18/00—Pattern recognition
- G06F18/20—Analysing
- G06F18/24—Classification techniques
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10032—Satellite or aerial image; Remote sensing
- G06T2207/10044—Radar image
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20004—Adaptive image processing
Landscapes
- Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Artificial Intelligence (AREA)
- Computer Vision & Pattern Recognition (AREA)
- Evolutionary Biology (AREA)
- Evolutionary Computation (AREA)
- Bioinformatics & Computational Biology (AREA)
- General Engineering & Computer Science (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
本发明公开基于自适应各向异性扩散的全极化SAR图像去噪方法，包括：S1、获取全极化SAR数据，并提取全极化SAR数据的协方差矩阵；S2、基于协方差矩阵提取极化散射熵；S3、基于极化散射熵进行梯度自适应滤波窗口的选择；S4、对协方差矩阵进行Freeman‑Durden分解，得到全极化SAR数据中各像素的布拉格散射、二次散射、体散射的功率；S5、基于各像素的布拉格散射、二次散射、体散射的功率，进行地物散射标记图的提取；S6、基于梯度自适应滤波窗口、地物散射标记图，采用散射机制的自适应各向异性扩散方法对全极化SAR图像去噪。本发明能够在抑制图像噪声的同时，充分地保留原始全极化SAR图像中的极化信息。
Description
技术领域
本发明涉及全极化SAR图像去噪技术领域，特别是涉及基于自适应各向异性扩散的全极化SAR图像去噪方法。
背景技术
极化合成孔径雷达(Polarimetric Synthetic Aperture Radar，PolSAR)利用天线的不同极化收发组合来测量地物目标，相比于合成孔径雷达可以获取更丰富的信息，具有重要研究价值和广泛的应用前景。但是由于成像机制的固有缺陷，PolSAR图像不可避免地会引入相干斑噪声，不仅导致图像视觉质量下降，而且对后续的分类、识别等处理造成了困难。因此相干斑抑制是极化SAR图像处理中必不可少的步骤。
合成孔径雷达(Synthetic Aperture Radar，SAR)诞生于20世纪60年代，是一种主动式的对地观测***，它弥补了传统雷达***中天线尺寸限制空间分辨率的缺陷。SAR是利用一个小天线沿着长线阵等速移动并辐射信号，然后把不同位置接收到的回波进行相干处理，从而获得雷达图像。目前，SAR已成为微波遥感中发展最为迅速和最有成效的传感器之一，其作为有源***，不仅具有全天候、全天时工作能力，而且可以在不同的微波频段、不同极化状态下对地观测，甚至可以透过地表和植被获取地表下信息。
极化合成孔径雷达(Polarimetric Synthetic Aperture Radar，PolSAR)是可以测量辐射信号极化特征的新型合成孔径雷达。PolSAR通过交替发射并同时接收极化状态正交的电磁波脉冲，不仅能测量地物回波的振幅，还能记录不同极化状态组合下回波的相位差，可用于分析任意极化状态的地物回波信息。PolSAR以其能够获得目标完全极化信息的独特优越性，极大的拓展了SAR的应用领域和价值。因此，PolSAR在农业(农作物分辨、生长状态估计等)、森林(生物估计、火灾救灾等)、地质(结构观测、矿物资源检测等)、水文(土壤温湿估计、反演等)、海洋学(波特性估计、溢油检测等)、搜救(失事飞机、轮船定位等)以及军事(军事目标识别等)等很多方面都得到了普遍的应用。针对PolSAR的研究也因此受到了普遍的关注和重视，已经成为国际雷达***与技术发展的热点之一。
由于***成像机理，相干斑(Speckle)是SAR(包括PolSAR)图像数据固有的一种确定性的干涉现象。相干斑在SAR图像上表现为像素值的剧烈随机起伏，每个像素点的强度值不再直接反映地面目标的雷达反射系数，导致图像的视觉质量变差，使利用SAR图像进行目标识别变得十分复杂，增加了目标识别的虚警概率和漏警概率。因此，极化SAR图像去噪成为了极化SAR图像解译中必不可少的重要步骤。在极化SAR图像去噪过程中，除了抑制相干斑噪声，图像极化信息的保持也具有十分重要的意义。
极化SAR相干斑抑制算法的发展经历了单极化SAR图像去噪、多极化SAR图像去噪与全极化SAR图像去噪三个阶段。目前国内外针对全极化SAR相干斑抑制的研究，最初主要基于空间域、频率域和极化域三种方向进行相干斑抑制，各方法的优缺点如下：
基于空间域的方法。基于空间域的方法主要根据数据的统计特征，计算邻域像元和目标像元之间的相似性，以此获得去噪后的估计影像。这类方法的优点是可减少由多视情况造成的分辨率损失，缺点是滤波后的影像较为模糊，丢失了许多细节信息，且滤波效果的好坏取决于滤波窗口尺寸的选择。
基于频率域的方法。基于频率域的方法主要利用小波变换理论与贝叶斯方法来对相干斑进行抑制。该类方法的优点是可以根据数据的特点来选择合适的小波基，缺点是难以区分高频信息，且滤波后的影像丢失极化信息，对后续的影像分类分割任务造成一定的影响。
基于极化域的方法。基于极化域的方法主要利用散射矩阵的最优组合，计算各个极化通道之间的线性相关性，最终得到去噪后的影像。该类方法的优点是能够处理多极化及多时相的数据，缺点是引入极化通道之间的相互串扰，且难以保持数据的极化信息。
随着技术的不断发展，并且极化方式从单极化向多极化的转变，更多新的相干斑抑制方法被提出来，如：基于散射模型的相干斑抑制方法、基于偏微分方程的相干斑抑制方法等。这些方法在相干斑抑制方面取得了不错的效果，但也存在着一些明显的缺点：在抑制相干斑噪声的同时，不能够很好的保留全极化SAR数据特有的极化信息，而SAR极化信息是地物目标识别的重要信息，丢失极化信息将造成数据信息的损失，给后续目标检测、分类等任务带来不便。
因此，提供一种基于自适应各向异性扩散的全极化SAR图像去噪方法，以实现在抑制噪声的同时，保持SAR图像的极化信息显得尤为重要。
发明内容
本发明的目的是提供基于自适应各向异性扩散的全极化SAR图像去噪方法，以解决现有技术中存在的技术问题，能够在抑制图像噪声的同时，充分地保留原始全极化SAR图像中的极化信息。
为实现上述目的，本发明提供了如下方案：本发明提供基于自适应各向异性扩散的全极化SAR图像去噪方法，包括：
S1、获取全极化SAR数据，并提取所述全极化SAR数据的协方差矩阵；
S2、基于所述全极化SAR数据的协方差矩阵，提取极化散射熵；
S3、基于所述极化散射熵，进行梯度自适应滤波窗口的选择；
S4、对所述全极化SAR数据的协方差矩阵进行Freeman-Durden分解，得到全极化SAR数据中各像素的布拉格散射、二次散射、体散射的功率；
S5、基于全极化SAR数据中各像素的布拉格散射、二次散射、体散射的功率，进行地物散射标记图的提取；
S6、基于梯度自适应滤波窗口、地物散射标记图，采用散射机制的自适应各向异性扩散方法对全极化SAR图像去噪。
优选地，所述步骤S3具体包括：
计算极化散射熵的垂直梯度、水平梯度，基于极化散射熵的垂直梯度、水平梯度，计算自适应滤波窗口的初始方向和尺寸；
基于自适应滤波窗口的初始尺寸构建水平滤波窗口，基于自适应滤波窗口的初始方向对所述水平滤波窗口进行旋转，并计算自适应滤波窗口各顶点的坐标，基于自适应滤波窗口的各顶点坐标，计算搜索窗口的顶点坐标；
遍历搜索窗口，获取搜索窗口、自适应滤波窗口中重叠的像素，计算所有重叠的像素的垂直梯度代数和、水平梯度代数和，利用垂直梯度代数和、水平梯度代数和更新自适应滤波窗口的方向；
计算搜索窗口、自适应滤波窗口中所有重叠的像素的垂直梯度绝对值之和、水平梯度绝对值之和，利用垂直梯度绝对值之和、水平梯度绝对值之和更新自适应滤波窗口的大小，直至自适应滤波窗口的尺寸满足预设阈值。
优选地，自适应滤波窗口的初始方向α0(i,j)如式7所示：
α0(i,j)＝arctan(Iy(i,j)/Ix(i,j))……………………7
其中，Ix(i,j)、Iy(i,j)分别为像素(i，j)处的极化散射熵的垂直梯度、水平梯度，i,j分别为像素的x轴、y轴坐标。
优选地，搜索窗口的顶点坐标如下式所示：
式中，Xmin、Xmax分别为搜索窗口顶点的横坐标的最小值、最大值；Ymin、Ymax分别为搜索窗口顶点的纵坐标的最小值、最大值；x1、x2、x3、x4分别为自适应滤波窗口四个顶点的横坐标；y1、y2、y3、y4分别为自适应滤波窗口四个顶点的纵坐标；ceil()代表向上取整，floor()代表向下取整。
优选地，所述步骤S5具体包括：
S5.1、基于全极化SAR数据中各像素的布拉格散射、二次散射、体散射的功率，对全极化SAR数据中的各像素进行初始聚类；
S5.2、基于初始聚类结果，采用Wishart距离对全极化SAR数据中的所有像素进行合并聚类；
S5.3、基于合并聚类结果，采用Wishart距离对全极化SAR数据中的各像素进行迭代分类；
S5.4、基于迭代分类结果，采用K-means聚类算法对全极化SAR数据中的各像素进行重新分配，生成地物散射标记图。
优选地，所述步骤S5.1中，将全极化SAR数据中所有像素依据布拉格散射、二次散射、体散射的功率的大小进行初始聚类。
优选地，所述步骤S5.3中，分别计算全极化SAR数据中每个像素与合并聚类结果中若干个聚类中心之间的Wishart距离，并根据距离计算结果对全极化SAR数据中的各像素进行重新分类。
优选地，所述步骤S6具体包括：
S6.1、基于梯度自适应滤波窗口、地物散射标记图，剔除梯度自适应滤波窗口中，与中心像素散射类型不同的邻域像素；
S6.2、基于梯度自适应滤波窗口中剩余的像素，计算梯度算子；
S6.3、基于梯度算子，计算自适应梯度阈值；
S6.4、基于自适应梯度阈值，构建自适应扩散模型；
S6.5、基于自适应扩散模型进行迭代计算，完成全极化SAR图像的去噪。
本发明公开了以下技术效果：
本发明提供一种基于自适应各向异性扩散的全极化SAR图像去噪方法，利用梯度计算滤波自适应滤波窗口，同时利用目标分解理论得到不同散射机制类型，从而得到地物散射标记图，基于自适应滤波窗口和地物散射标记图，通过改进的自适应各向异性扩散方程对全极化SAR图像进行相干斑抑制，从而能够在抑制图像噪声的同时，充分地保留原始全极化SAR图像中的极化信息。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。
图1为本发明基于自适应各向异性扩散的全极化SAR图像去噪方法流程图；
图2为本发明实施例中基于自适应各向异性扩散的全极化SAR图像去噪方法示意图；
图3为本发明实施例中极化散射熵提取结果示意图；
图4为本发明实施例中全极化SAR数据的协方差矩阵的Freeman-Durden分解结果示意图；其中，图4(a)、图4(b)、图4(c)分别为Ps、Pd、Pv分量的分解结果示意图；
图5为本发明实施例中地物散射标记图示意图；
图6为本发明实施例中全极化SAR图像的去噪结果示意图；其中，图6(a)、图6(b)分别为HH通道、VV通道的去噪结果示意图；
图7为本发明实施例中不同去噪方法的去噪结果对比图；其中，图7(a)为原始图像，图7(b)、图7(c)、图7(d)、图7(e)、图7(f)、图7(g)、图7(h)分别为本发明方法、SARD、Re-Lee、Boxcar、Gauss、IDAN、Mean-shift去噪结果示意图；
图8为本发明实施例中原始图像与通过本发明方法去噪后数据在不同计划状态下的极化特征图；其中，图8(a)、图8(b)分别为原始图像共极化通道(油膜区)、去噪后图像共极化通道(油膜区)极化特征图；图8(c)、图8(d)分别为原始图像交叉极化通道(油膜区)、去噪后图像交叉极化通道(油膜区)极化特征图；图8(e)、图8(f)分别为原始图像共极化通道(海水区)、去噪后图像共极化通道(海水区)极化特征图；图8(e)、图8(f)分别为原始图像共极化通道(海水区)、去噪后图像共极化通道(海水区)极化特征图；图8(g)、图8(h)分别为原始图像交叉极化通道(海水区)、去噪后图像交叉极化通道(海水区)极化特征图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
为使本发明的上述目的、特征和优点能够更加明显易懂，下面结合附图和具体实施方式对本发明作进一步详细的说明。
本发明中涉及到的专业术语解释如下：
全极化SAR：包含HH、HV、VH、VV四种极化方式的合成孔径雷达传感器，具有全天时、全天候、可穿云透雾的特点。
相干斑噪声：由合成孔径雷达(SAR)包括极化SAR成像机制产生的，由大量散射单元反射波的相干叠加引起的。相干斑使相邻像素间的信号强度发生变化，视觉上表现为颗粒状的噪声。
各向异性：是指物质的全部或部分化学、物理等性质随着方向的改变而有所变化，在不同的方向上呈现出差异的性质。表现在图像上就是以某点为中心，该点各个方向上图像的特性均有所不同。
参照图1-2所示，本实施例提供基于自适应各向异性扩散的全极化SAR图像去噪方法，包括如下步骤：
S1、获取全极化SAR数据，并提取所述全极化SAR数据的协方差矩阵；
所述全极化SAR数据的协方差矩阵的提取包括极化散射矩阵S与极化协方差矩阵C3的提取。
极化散射矩阵S能完整描述SAR目标的电磁散射特性，其定义如式(1)所示：
其中，SHV代表水平极化发射，垂直极化接收；SVH代表垂直极化发射，水平极化接收；SHH代表水平极化发射与接收；SVV代表垂直极化发射与接收。极化散射矩阵的对角元素代表相同极化方式的入射场和散射场关系，称为“同极化”项。非对角元素代表正交极化方式的入射角和散射场关系，称为“交叉极化”项。
当满足互易性定理时，有SHV＝SVH，将散射矩阵S利用Lexicographer基分解得到极化散射矢量
其定义如式(2)所示：
式中，因子
可以保证目标矢量的范数与散射总功率(Span)相等，即满足“总功率守恒”。
当SHV＝SVH时，目标满足互易性定理，此时极化散射矢量
可简化为三维极化散射矢量，其定义如式(3)所示：
极化协方差矩阵C3由极化散射矢量
与极化散射矢量
的共轭转置矢量
的外积生成，C3定义如式(4)所示：
其中，*T，*分别代表共轭转置和共轭，<·>表示整体平均值。
S2、基于所述全极化SAR数据的协方差矩阵，提取极化散射熵；
通过对3×3埃尔米特极化协方差矩阵C3进行对角化，计算得到极化协方差矩阵C3的特征值和特征矢量，如式(5)所示：
式中，∑为C3的特征值矩阵，∑是3×3的非负实对角阵，∑的主对角线元素为C3的特征值λ1,λ2,λ3，有特征值λ1≥λ2≥λ3≥0，UC为C3的特征向量矩阵，UC＝[e1e2e3]是3×3特殊酉矩阵SU(3)，其中e1,e2,e3分别是C3的特征向量。
极化散射熵H定义如式(6)所示：
式中，
极化散射熵描述分布式散射体的随机散射程度，极化散射熵的值较低时表现为单次散射，极化散射熵的值较高时表现为随机散射。例如，当雷达波束照射到清洁海面时，主要发生的是表面布拉格散射，散射机制相对单一，故H的值较小；当雷达波束照射到油膜覆盖的海水时，产生表面布拉格散射的同时，还存在镜面散射，散射过程相对复杂，散射随机性较强，故H的值较大。提取的极化散射熵如图3所示。
S3、基于所述极化散射熵，进行梯度自适应滤波窗口的选择；
本发明提出一种梯度自适应窗口选择方法，该方法利用极化散射熵的梯度变化来反映图像细节变化，依据梯度自适应的选择滤波窗口的方向和大小，使梯度较小的均质区域窗口较大，梯度较大的非均质区域窗口较小，边缘区域窗口较狭窄。
梯度自适应滤波窗口选择的具体方法包括：
S3.1、计算极化散射熵的垂直梯度Ix(i,j)、水平梯度Iy(i,j)，其中，i,j分别为像素的x轴、y轴坐标。
S3.2、基于极化散射熵的垂直梯度Ix(i,j)、水平梯度Iy(i,j)，计算自适应滤波窗口Qchoose的初始方向和尺寸。
设a为最大的方形窗口尺寸，α0(i,j)为自适应滤波窗口的初始方向，如式(7)所示，w0(i,j)为自适应滤波窗口的初始宽度，如式(8)所示，l0(i,j)为自适应滤波窗口的初始长度，如公式(9)所示。
α0(i,j)＝arctan(Iy(i,j)/Ix(i,j))……………………(7)
w0(i,j)＝a/(|Ix(i,j)|+1)………………………(8)
l0(i,j)＝a/(|Iy(i,j)|+1)………………………(9)
其中，Ix(i,j)、Iy(i,j)分别为像素(i，j)处的极化散射熵的垂直梯度、水平梯度。
S3.3、构建长lo、宽wo的水平滤波窗口，基于自适应滤波窗口的初始方向对水平滤波窗口进行旋转，并计算自适应滤波窗口Qchoose各顶点的坐标；给定以单个像素为中心Qchoose的初始宽wo、长lo与方向αo，则转换公式如公式(10)所示：
式中，αo为自适应滤波窗口的初始方向；x1、x2、x3、x4分别为自适应滤波窗口四个顶点的横坐标，y1、y2、y3、y4分别为自适应滤波窗口四个顶点的纵坐标。
S3.4、基于自适应滤波窗口Qchoose的各顶点坐标，计算搜索窗口Qsearch的顶点坐标，其定义如式(11)所示：
式中，ceil()代表向上取整，floor()代表向下取整。
S3.5、遍历搜索窗口Qsearch，获取搜索窗口、自适应滤波窗口中重叠的像素，计算所有重叠的像素的垂直梯度代数和Iys(i,j)与水平梯度代数和Ixs(i,j)，其定义如式(12)所示。
式中，N表示搜索窗口、自适应滤波窗口中重叠的像素的数量。
S3.6、利用垂直梯度代数和Iys(i,j)、水平梯度代数和Ixs(i,j)更新自适应滤波窗口的方向θm，定义如式(13)所示：
θm(i,j)＝arctan(Iys(i,j)/Ixs(i,j))…………………(13)。
S3.7、计算搜索窗口、自适应滤波窗口中所有重叠的像素的垂直梯度绝对值之和Iyc(i,j)、水平梯度绝对值之和Ixc(i,j)，其定义如式(14)所示：
S3.8、利用垂直梯度绝对值之和Iyc(i,j)、水平梯度绝对值之和Ixc(i,j)更新自适应滤波窗口的大小，设更新后的选择窗口宽为wm(i,j)，长为lm(i,j)，则其定义如式(15)所示：
式中，wm(i,j)、lm(i,j)分别为更新后的自适应滤波窗口的宽、长。
S3.9、当自适应滤波窗口的尺寸满足公式(16)时，停止更新，自适应滤波窗口宽为wm+1(i,j)，长为lm+1(i,j)，方向为θm+1：
式中，T1和T2为预设阈值。
S4、对所述全极化SAR数据的协方差矩阵进行Freeman-Durden分解，得到全极化SAR数据中各像素的布拉格散射、二次散射、体散射的功率；
Freeman-Durden分解是将三分量散射机制模型用于极化SAR观测的技术，它不需要使用任何地面测量数据。在满足互易性定理的情况下，使用Freeman-Durden分解将目标全极化SAR数据的协方差矩阵进行重新表示，其定义如式(17)所示：
式中，fv、fd、fs均为系数，与fv、fd、fs相乘的三个矩阵分别代表体散射、二次散射与布拉格散射的协方差矩阵，α、β分别为二次散射参数、布拉格散射参数，α、β的值由
的正负确定。
在Freeman-Durden分解中，fv、fd、fs系数满足如式(18)所示的等式关系。其中，当
时，α＝-1；当
时，β＝1。
根据极化总功率守恒定理，求得Freeman-Durden分解的最终表达式，其定义如式(19)所示，Freeman-Durden分解的结果如图4所示。
式中，Ps、Pd、Pv分别表示布拉格散射、二次散射、体散射的功率。
S5、基于全极化SAR数据中各像素的布拉格散射、二次散射、体散射的功率，进行地物散射标记图的提取；
在极化SAR图像中，不同的地物具有不同的散射特征，而不同的地物散射特性可以用来描述不同类别地物或相同类别地物在极化散射方面的差异性，因此对极化SAR图像的去噪应该基于相同散射机制进行。
提取地物散射标记图具体包括如下步骤：
S5.1、基于全极化SAR数据中各像素的布拉格散射、二次散射、体散射的功率，对全极化SAR数据中的各像素进行初始聚类；
在不同散射类型中，将全极化SAR数据中所有像素依据布拉格散射、二次散射、体散射的功率的大小进行初始聚类，得到30个初始类别。
S5.2、基于初始聚类结果，采用Wishart距离对全极化SAR数据中的所有像素进行合并聚类；
距离度量公式如式(20)所示：
式中，Vp与Vq分别代表初始聚类结果中第p类与第q类的协方差矩阵的类均值，Dpq代表初始聚类结果中第p类与第q类之间的距离，若初始聚类结果中某两类之间的Dpq最短，则将这两类进行合并；本实施例中，最终不同散射类型的所有初始聚类将被合并为9类。
S5.3、基于合并聚类结果，采用Wishart距离对全极化SAR数据中的各像素进行迭代分类；
分别计算全极化SAR数据中每个像素与合并聚类结果中的9个聚类中心之间的Wishart距离，并根据距离计算结果对全极化SAR数据中的各像素进行重新分类，为了使分类结果更加收敛，本实施例利用迭代Wishart分类器迭代计算5-10次。Wishart距离度量公式如式(21)所示：
式中，C代表像素的极化协方差矩阵，Vf代表合并聚类结果中的第f类聚类中心。
S5.4、基于迭代分类结果，采用K-means聚类算法对全极化SAR数据中的各像素进行重新分配，生成地物散射标记图。
K-means聚类算法是一种迭代求解的聚类分析算法，具体包括：
首先，基于迭代分类结果，随机选取若干个像素作为初始的K-means聚类中心；本实施例中，随机选取3个像素。
然后，计算全极化SAR数据中的各像素与各个K-means聚类中心之间的绝对差之和，基于绝对差之和，将全极化SAR数据中的每个像素分配给距离最近的K-means聚类中心。
K-means聚类中心以及分配给该K-means聚类中心的像素代表一个聚类。每分配一个像素，K-means聚类中心根据聚类中现有的像素重新计算，这个过程不断重复直到没有对象被重新分配给不同的聚类为止，最终生成的标记后的地物散射特征图，即地物散射标记图，如图5所示。K-means距离度量公式如式(22)所示：
式中，P为数据的维度，此处P＝1；xj为像素；cj为聚类中心。
S6、基于梯度自适应滤波窗口、地物散射标记图，采用散射机制的自适应各向异性扩散方法对全极化SAR图像去噪。具体包括：
S6.1、基于梯度自适应滤波窗口、地物散射标记图，剔除梯度自适应滤波窗口中，与中心像素散射类型不同的邻域像素；
S6.2、基于梯度自适应滤波窗口中剩余的像素，计算梯度算子
如式(23)所示：
式中，b∈[1,B]且b为整数；I(xb,yb)代表邻域像素，I(x,y)代表中心像素；B为梯度自适应滤波窗口中与中心像素散射类型相同的邻域像素的个数；η为邻域像素与中心像素的距离。
S6.3、基于梯度算子，计算自适应梯度阈值k，如式(24)所示：
式中，
为梯度自适应滤波窗口内梯度的平均值。
S6.4、基于自适应梯度阈值，构建自适应扩散模型，如式(25)所示：
式中，g()为自适应扩散函数，s表示自适应扩散函数变量。
S6.5、基于自适应扩散模型进行迭代计算，完成全极化SAR图像的去噪，如式(26)所示：
式中，m为迭代次数，Δt为时间间隔，λ(I0-Im(x,y))为保真项，λ为保真系数；通过保真项，能够使去噪后的全极化SAR图像尽可能接近原始影像，有效保持了图像中的极化信息。全极化SAR图像的去噪结果如图6所示。
为进一步验证本发明基于自适应各向异性扩散的全极化SAR图像去噪方法的有效性和可靠性，本实施例分别利用本发明方法和SARD、Re-Lee、Boxcar、Gauss、IDAN、Mean-shift 6种经典的去噪方法获得的极化SAR图像去噪结果进行对比，去噪结果对比如图7所示。其中，分别使用极化特征图、等效视数ENL、辐射分辨率、均值、方差和协方差矩阵差异度6个去噪评价指标对相干斑抑制结果进行定量与定性地评估，图8是极化特征图定性分析结果，表1是等效视数定量评定结果，表2是辐射分辨率定量评定结果，表3是图像均值定量评定结果，表4是图像方差定量评价结果，表5是协方差矩阵差异度定量评价结果。其中，极化特征图是收发天线在不同极化状态下的回波功率图，是表征目标极化特性的有利工具，通常采用对比极化滤波前后极化特征图一致性程度的方法，直观地评估极化滤波算法的极化特征保持性，极化特征图越相近，表明算法的极化特征保持性越强，算法越可靠；等效视数ENL能够衡量滤波后图像的平滑程度，如果某均质区域的ENL越大，则该区域越平滑，滤波性能也就越好；辐射分辨率表示区分极化SAR目标后向散射系数的能力，是衡量极化SAR***区分相邻分布目标的能力的一种量度。去噪的好坏直接影响极化SAR图像解译和定量化应用，辐射分辨率的大小由消除相干斑噪声的多少决定，因此好的极化SAR去噪方法能提高辐射分辨率；图像均值是整个图像的平均强度，反映了图像的平均灰度，即图像所包含目标的平均后向散射系数；图像方差代表了图像区域中所有点偏离均值的程度，反映了图像的不均匀性；图像的均值和方差是反映图像整体特征的指标，所以应当尽量保持图像的均值，同时减少图像的方差。协方差矩阵差异度Dod能够衡量滤波前后影像的协方差矩阵的保真程度，通过比较滤波前后影像的协方差矩阵差异度可以对极化SAR去噪方法的极化信息保持性能进行定量评价，如果Dod越小，则协方差矩阵的变化就越小，极化信息保持程度也就越高。
表1
表2
表3
表4
表5
首先从定性的角度评价，从图7可以直观的看到本发明方法的滤波结果提取的不同极化状态的极化特征图与原始图像提取的极化特征图几乎没有区别，说明本发明方法的极化信息保持能力较强。另外从定量的角度评价，从表1与表2看出本发明方法的等效视数和辐射分辨率均明显优于其他方法，从表3中可以看出本发明方法滤波图像的均值与原始图像均值保持一致，从表4看出Re-Lee和IDAN方法滤波结果的方差更小一些，但是从表3可知，这两种方法的滤波影像均值与原始图像相差较大，故综合表3和表4可以看出，本发明方法的滤波影像要优于其他方法。从表5可以看出，本发明方法的协方差矩阵差异度最小，故本发明方法的极化信息保持性优于其它方法。这表明基于自适应各向异性扩散方法相对于其他滤波方法能更有效的保持极化SAR数据的极化信息。因此，本发明在相干斑抑制和极化信息保持方面均优于其它滤波方法。
以上所述的实施例仅是对本发明的优选方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形和改进，均应落入本发明权利要求书确定的保护范围内。
Claims (5)
1.基于自适应各向异性扩散的全极化SAR图像去噪方法，其特征在于，包括：
S1、获取全极化SAR数据，并提取所述全极化SAR数据的协方差矩阵；
S2、基于所述全极化SAR数据的协方差矩阵，提取极化散射熵；
S3、基于所述极化散射熵，进行梯度自适应滤波窗口的选择；
S4、对所述全极化SAR数据的协方差矩阵进行Freeman-Durden分解，得到全极化SAR数据中各像素的布拉格散射、二次散射、体散射的功率；
S5、基于全极化SAR数据中各像素的布拉格散射、二次散射、体散射的功率，进行地物散射标记图的提取；
S6、基于梯度自适应滤波窗口、地物散射标记图，采用散射机制的自适应各向异性扩散方法对全极化SAR图像去噪；
所述步骤S3具体包括：
计算极化散射熵的垂直梯度、水平梯度，基于极化散射熵的垂直梯度、水平梯度，计算自适应滤波窗口的初始方向和尺寸；
构建长lo、宽wo的水平滤波窗口，转换公式如下式所示：
式中，αo为自适应滤波窗口的初始方向；x1、x2、x3、x4分别为自适应滤波窗口四个顶点的横坐标，y1、y2、y3、y4分别为自适应滤波窗口四个顶点的纵坐标；
基于自适应滤波窗口的初始尺寸构建水平滤波窗口，基于自适应滤波窗口的初始方向对所述水平滤波窗口进行旋转，并计算自适应滤波窗口各顶点的坐标，基于自适应滤波窗口的各顶点坐标，计算搜索窗口的顶点坐标；
遍历搜索窗口，获取搜索窗口、自适应滤波窗口中重叠的像素，计算所有重叠的像素的垂直梯度代数和、水平梯度代数和，利用垂直梯度代数和、水平梯度代数和更新自适应滤波窗口的方向；
计算搜索窗口、自适应滤波窗口中所有重叠的像素的垂直梯度绝对值之和、水平梯度绝对值之和，利用垂直梯度绝对值之和、水平梯度绝对值之和更新自适应滤波窗口的大小，直至自适应滤波窗口的尺寸满足预设阈值；
更新自适应滤波窗口的大小，设更新后的选择窗口宽为wm(i,j)，长为lm(i,j)，则其定义如下式所示：
式中，wm(i,j)、lm(i,j)分别为更新后的自适应滤波窗口的宽、长；
自适应滤波窗口的初始方向α0(i,j)如式7所示：
α0(i,j)＝arctan(Iy(i,j)/Ix(i,j))……………………7
其中，Ix(i,j)、Iy(i,j)分别为像素(i，j)处的极化散射熵的垂直梯度、水平梯度，i,j分别为像素的x轴、y轴坐标；
搜索窗口的顶点坐标如下式所示：
式中，Xmin、Xmax分别为搜索窗口顶点的横坐标的最小值、最大值；Ymin、Ymax分别为搜索窗口顶点的纵坐标的最小值、最大值；x1、x2、x3、x4分别为自适应滤波窗口四个顶点的横坐标；y1、y2、y3、y4分别为自适应滤波窗口四个顶点的纵坐标；ceil()代表向上取整，floor()代表向下取整。
2.根据权利要求1所述的基于自适应各向异性扩散的全极化SAR图像去噪方法，其特征在于，所述步骤S5具体包括：
S5.1、基于全极化SAR数据中各像素的布拉格散射、二次散射、体散射的功率，对全极化SAR数据中的各像素进行初始聚类；
S5.2、基于初始聚类结果，采用Wishart距离对全极化SAR数据中的所有像素进行合并聚类；
S5.3、基于合并聚类结果，采用Wishart距离对全极化SAR数据中的各像素进行迭代分类；
S5.4、基于迭代分类结果，采用K-means聚类算法对全极化SAR数据中的各像素进行重新分配，生成地物散射标记图。
3.根据权利要求2所述的基于自适应各向异性扩散的全极化SAR图像去噪方法，其特征在于，所述步骤S5.1中，将全极化SAR数据中所有像素依据布拉格散射、二次散射、体散射的功率的大小进行初始聚类。
4.根据权利要求2所述的基于自适应各向异性扩散的全极化SAR图像去噪方法，其特征在于，所述步骤S5.3中，分别计算全极化SAR数据中每个像素与合并聚类结果中若干个聚类中心之间的Wishart距离，并根据距离计算结果对全极化SAR数据中的各像素进行重新分类。
5.根据权利要求1所述的基于自适应各向异性扩散的全极化SAR图像去噪方法，其特征在于，所述步骤S6具体包括：
S6.1、基于梯度自适应滤波窗口、地物散射标记图，剔除梯度自适应滤波窗口中，与中心像素散射类型不同的邻域像素；
S6.2、基于梯度自适应滤波窗口中剩余的像素，计算梯度算子；
S6.3、基于梯度算子，计算自适应梯度阈值；
S6.4、基于自适应梯度阈值，构建自适应扩散模型；
S6.5、基于自适应扩散模型进行迭代计算，完成全极化SAR图像的去噪。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110117113.XA CN112950492B (zh) | 2021-01-28 | 2021-01-28 | 基于自适应各向异性扩散的全极化sar图像去噪方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110117113.XA CN112950492B (zh) | 2021-01-28 | 2021-01-28 | 基于自适应各向异性扩散的全极化sar图像去噪方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN112950492A CN112950492A (zh) | 2021-06-11 |
| CN112950492B true CN112950492B (zh) | 2022-04-29 |
Family
ID=76238497
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110117113.XA Expired - Fee Related CN112950492B (zh) | 2021-01-28 | 2021-01-28 | 基于自适应各向异性扩散的全极化sar图像去噪方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN112950492B (zh) |
Families Citing this family (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN116091345B (zh) * | 2022-12-23 | 2023-10-03 | 重庆大学 | 一种基于局部熵和保真项的各向异性扩散医学图像去噪方法、***及存储介质 |
Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102096921A (zh) * | 2011-01-10 | 2011-06-15 | 西安电子科技大学 | 基于邻域对数比值及各向异性扩散的sar图像变化检测方法 |
| CN103778613A (zh) * | 2014-02-21 | 2014-05-07 | 武汉大学 | 一种窗口自适应的极化sar影像滤波方法 |
| CN104318246A (zh) * | 2014-10-20 | 2015-01-28 | 西安电子科技大学 | 基于深度自适应脊波网络的极化sar图像分类 |
| CN104951789A (zh) * | 2015-07-15 | 2015-09-30 | 电子科技大学 | 一种基于全极化sar图像的快速滑坡提取方法 |
| CN105160353A (zh) * | 2015-08-18 | 2015-12-16 | 西安电子科技大学 | 基于多特征集的极化sar数据地物分类方法 |
| CN106778884A (zh) * | 2016-12-23 | 2017-05-31 | 辽宁工程技术大学 | 基于散射熵和三分量分解平面的极化sar图像分类方法 |
| CN107229933A (zh) * | 2017-05-11 | 2017-10-03 | 西安电子科技大学 | 自适应体散射模型的freeman/特征值分解方法 |
| CN107644408A (zh) * | 2017-09-27 | 2018-01-30 | 安徽大学 | 一种基于各项异性扩散的全极化雷达影像去噪方法 |
| CN108387898A (zh) * | 2018-02-27 | 2018-08-10 | 河海大学 | 基于同质点选取的极化特征参数提取方法 |
| CN108876746A (zh) * | 2018-06-28 | 2018-11-23 | 广东工业大学 | 一种基于各向异性扩散方程的sar图像降斑方法及*** |
| CN110363105A (zh) * | 2019-06-25 | 2019-10-22 | 电子科技大学 | 一种全极化sar图像相干斑抑制方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103824084B (zh) * | 2014-03-12 | 2016-11-02 | 西安电子科技大学 | 基于sdit和svm的极化sar图像分类方法 |
| CN109740475B (zh) * | 2018-12-25 | 2020-04-21 | 杭州世平信息科技有限公司 | 一种遥感图像地面场景分类方法 |
-
2021
- 2021-01-28 CN CN202110117113.XA patent/CN112950492B/zh not_active Expired - Fee Related
Patent Citations (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN102096921A (zh) * | 2011-01-10 | 2011-06-15 | 西安电子科技大学 | 基于邻域对数比值及各向异性扩散的sar图像变化检测方法 |
| CN103778613A (zh) * | 2014-02-21 | 2014-05-07 | 武汉大学 | 一种窗口自适应的极化sar影像滤波方法 |
| CN104318246A (zh) * | 2014-10-20 | 2015-01-28 | 西安电子科技大学 | 基于深度自适应脊波网络的极化sar图像分类 |
| CN104951789A (zh) * | 2015-07-15 | 2015-09-30 | 电子科技大学 | 一种基于全极化sar图像的快速滑坡提取方法 |
| CN105160353A (zh) * | 2015-08-18 | 2015-12-16 | 西安电子科技大学 | 基于多特征集的极化sar数据地物分类方法 |
| CN106778884A (zh) * | 2016-12-23 | 2017-05-31 | 辽宁工程技术大学 | 基于散射熵和三分量分解平面的极化sar图像分类方法 |
| CN107229933A (zh) * | 2017-05-11 | 2017-10-03 | 西安电子科技大学 | 自适应体散射模型的freeman/特征值分解方法 |
| CN107644408A (zh) * | 2017-09-27 | 2018-01-30 | 安徽大学 | 一种基于各项异性扩散的全极化雷达影像去噪方法 |
| CN108387898A (zh) * | 2018-02-27 | 2018-08-10 | 河海大学 | 基于同质点选取的极化特征参数提取方法 |
| CN108876746A (zh) * | 2018-06-28 | 2018-11-23 | 广东工业大学 | 一种基于各向异性扩散方程的sar图像降斑方法及*** |
| CN110363105A (zh) * | 2019-06-25 | 2019-10-22 | 电子科技大学 | 一种全极化sar图像相干斑抑制方法 |
Non-Patent Citations (5)
| Title |
|---|
| 一种自适应的混合Freeman /Eigenvalue 极化分解模型;何连，秦其明，任华忠;《国土资源遥感》;20170630;第29卷(第2期);全文 * |
| 基于GA-WNN 的极化SAR 海洋溢油检测方法研究;陈伟民, 丁亚雄, 宋冬梅, 王 斌等;《海洋科学》;20181231;第42卷(第1期);全文 * |
| 基于场景异质性的极化SAR相干斑滤波方法;行晓黎;《中国博士学位论文全文数据库信息科技辑》;20180715;第2018年卷(第7期);全文 * |
| 基于改进扩散系数的各向异性扩散相干斑抑制算法;李艳华;《中国优秀硕士学位论文数据库信息科技辑》;20160315;第2016年卷(第3期);全文 * |
| 基于新极化特征参数的SAR 海洋溢油检测;任慧敏，宋冬梅等;《遥感技术与应用》;20200831;第35卷(第4期);全文 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN112950492A (zh) | 2021-06-11 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| Huang et al. | Under-foliage object imaging using SAR tomography and polarimetric spectral estimators | |
| Frey et al. | Analyzing tomographic SAR data of a forest with respect to frequency, polarization, and focusing technique | |
| Ertin et al. | GOTCHA experience report: Three-dimensional SAR imaging with complete circular apertures | |
| CN112558066B (zh) | 双极化sar图像*** | |
| Cloude et al. | Target detection beneath foliage using polarimetric synthetic aperture radar interferometry | |
| CN101369019A (zh) | 基于极化数据融合的极化干涉合成孔径雷达三维成像方法 | |
| Chen et al. | Mitigation of azimuth ambiguities in spaceborne stripmap SAR images using selective restoration | |
| Shahrezaei et al. | Fractal analysis and texture classification of high-frequency multiplicative noise in SAR sea-ice images based on a transform-domain image decomposition method | |
| Rambour et al. | Introducing spatial regularization in SAR tomography reconstruction | |
| Stojanovic et al. | Joint space aspect reconstruction of wide-angle SAR exploiting sparsity | |
| Iqbal et al. | SAR image despeckling by selective 3D filtering of multiple compressive reconstructed images | |
| Nielsen et al. | Direction-of-arrival estimation for radar ice sounding surface clutter suppression | |
| CN112950492B (zh) | 基于自适应各向异性扩散的全极化sar图像去噪方法 | |
| Volosyuk et al. | Phenomenological description of coherent radar images based on the concepts of the measure of set and stochastic integral | |
| Omati et al. | The comparative study of three nonparametric methods of SAR tomography for building reconstruction | |
| Stagliano et al. | Ship detection from SAR images based on CFAR and wavelet transform | |
| Demirci et al. | An investigation of the performances of polarimetric target decompositions using GB-SAR imaging | |
| Brigui et al. | New SAR target imaging algorithm based on oblique projection for clutter reduction | |
| Dos Santos et al. | Polarimetric responses and scattering mechanisms of tropical forest in the Brazilian Amazon | |
| Karachristos et al. | A Review on PolSAR Decompositions for Feature Extraction | |
| Azzaz et al. | Classification of radar echoes using fractal geometry | |
| Aghababaee et al. | Full 3D DEM generation in urban area by improving estimation from SAR tomography | |
| Wecklich et al. | Production of a global forest/non-forest map utilizing TanDEM-X interferometric SAR data | |
| Akbari | Multitemporal analysis of multipolarization synthetic aperture radar images for robust surface change detection | |
| Prakash et al. | Emerging techniques of polarimetric interferometric synthetic aperture radar for scattering-based characterization |
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: 20220429 |