CN114355348B - Sar干涉图小波降噪处理方法及其处理装置 - Google Patents
Sar干涉图小波降噪处理方法及其处理装置 Download PDFInfo
- Publication number
- CN114355348B CN114355348B CN202210020853.6A CN202210020853A CN114355348B CN 114355348 B CN114355348 B CN 114355348B CN 202210020853 A CN202210020853 A CN 202210020853A CN 114355348 B CN114355348 B CN 114355348B
- Authority
- CN
- China
- Prior art keywords
- wavelet
- noise
- filtering
- denoising
- interferogram
- 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
- 238000012545 processing Methods 0.000 title claims abstract description 42
- 238000003672 processing method Methods 0.000 title claims abstract description 16
- 238000001914 filtration Methods 0.000 claims abstract description 134
- 238000000354 decomposition reaction Methods 0.000 claims abstract description 55
- 239000011159 matrix material Substances 0.000 claims abstract description 17
- 230000002087 whitening effect Effects 0.000 claims abstract description 17
- 238000012360 testing method Methods 0.000 claims abstract description 6
- 238000000034 method Methods 0.000 claims description 144
- 230000006870 function Effects 0.000 claims description 49
- 230000009467 reduction Effects 0.000 claims description 34
- 230000008569 process Effects 0.000 claims description 15
- 238000004364 calculation method Methods 0.000 claims description 12
- 238000007689 inspection Methods 0.000 claims description 5
- 230000014509 gene expression Effects 0.000 claims description 4
- 238000004590 computer program Methods 0.000 claims description 3
- 238000003860 storage Methods 0.000 claims description 3
- 238000010276 construction Methods 0.000 claims 1
- 238000002474 experimental method Methods 0.000 description 55
- 238000004422 calculation algorithm Methods 0.000 description 20
- 238000004088 simulation Methods 0.000 description 20
- 230000003044 adaptive effect Effects 0.000 description 17
- 230000000694 effects Effects 0.000 description 13
- 238000011160 research Methods 0.000 description 10
- 230000009466 transformation Effects 0.000 description 10
- 238000010586 diagram Methods 0.000 description 9
- 230000008901 benefit Effects 0.000 description 8
- 238000003384 imaging method Methods 0.000 description 8
- 230000002860 competitive effect Effects 0.000 description 7
- 238000005065 mining Methods 0.000 description 7
- 230000000007 visual effect Effects 0.000 description 7
- 238000001514 detection method Methods 0.000 description 6
- 239000000654 additive Substances 0.000 description 5
- 230000000996 additive effect Effects 0.000 description 5
- 238000004458 analytical method Methods 0.000 description 4
- 238000009499 grossing Methods 0.000 description 4
- 230000014759 maintenance of location Effects 0.000 description 4
- 238000010587 phase diagram Methods 0.000 description 4
- XLYOFNOQVPJJNP-UHFFFAOYSA-N water Substances O XLYOFNOQVPJJNP-UHFFFAOYSA-N 0.000 description 4
- 238000010835 comparative analysis Methods 0.000 description 3
- 238000009826 distribution Methods 0.000 description 3
- 238000011156 evaluation Methods 0.000 description 3
- 238000005457 optimization Methods 0.000 description 3
- 238000012795 verification Methods 0.000 description 3
- 230000006399 behavior Effects 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 238000009792 diffusion process Methods 0.000 description 2
- 238000005516 engineering process Methods 0.000 description 2
- 230000000717 retained effect Effects 0.000 description 2
- 238000013519 translation Methods 0.000 description 2
- 238000012935 Averaging Methods 0.000 description 1
- 241000223025 Caragana microphylla Species 0.000 description 1
- 230000002146 bilateral effect Effects 0.000 description 1
- 230000008602 contraction Effects 0.000 description 1
- 230000001186 cumulative effect Effects 0.000 description 1
- 238000013016 damping Methods 0.000 description 1
- 230000007423 decrease Effects 0.000 description 1
- 230000003247 decreasing effect Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 230000009977 dual effect Effects 0.000 description 1
- 238000007429 general method Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 230000007774 longterm Effects 0.000 description 1
- 238000012423 maintenance Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000012544 monitoring process Methods 0.000 description 1
- 239000002245 particle Substances 0.000 description 1
- 238000004321 preservation Methods 0.000 description 1
- 238000011158 quantitative evaluation Methods 0.000 description 1
- 238000011946 reduction process Methods 0.000 description 1
- 239000007787 solid Substances 0.000 description 1
- 238000001228 spectrum Methods 0.000 description 1
- 230000001629 suppression Effects 0.000 description 1
- 238000010998 test method Methods 0.000 description 1
Images
Landscapes
- Image Processing (AREA)
Abstract
本发明提供了一种SAR干涉图小波降噪处理方法,包括:确定SAR干涉图中的信号点和噪声点,构建噪点识别矩阵,仅对噪声点进行滤波降噪处理;提取出含有噪声的干涉图数据,分解为实部和虚部,利用Biorthogonal小波基对实部和虚部分别进行小波分解,并采用白化检验的方式使小波分解的层数达到最优;利用改进的小波阈值函数和Birge‑Massart惩罚函数自适应获取的分层阈值对小波系数进行降噪处理;对滤波降噪后的实部和虚部进行小波重构,获取滤波降噪后的干涉图,并结合SAR干涉图中的信号点,获取最终的降噪干涉图。本发明还提供了一种SAR干涉图小波降噪处理装置。该发明既能够有效降噪,又能够保持边缘细节信息。
Description
技术领域
本发明涉及图像处理领域。更具体地说,本发明涉及一种SAR干涉图小波降噪处理方法及其处理装置。
背景技术
由于散射体与传感器之间距离不同,影响回波信号的相干性,使SAR(SyntheticAperture Radar,合成孔径雷达)成像产生相干斑噪声,这在一定程度上降低了SAR影像的质量,可导致SAR影像获取的地物特征存在失真现象。因此,对SAR影像进行滤波降噪,可提升SAR影像的判读能力,获取更加真实可靠的相位信息。
目前,SAR降噪技术主要有两类:一类是成像前的多视处理技术,一类是成像后的滤波处理技术。成像前的多视处理技术一般使用多视平均的方法,其原理是对返回信号的方位向或距离向降低处理器的带宽,再将信号频谱分割成若干部分,分别成像后进行强度简单叠加得到最终影像。该类方法达到了抑制噪声的目的,但在一定程度上影响了SAR影像的空间分辨率,且对相干斑噪声的滤除不是很完美。成像后的滤波处理技术又可分为两种:(1)空间域滤波方法,主要包括均值、中值和统计类滤波方法。该类方法基于影像像素间空间上的相关性对SAR影像的噪声进行去除,以达到降噪目的,应用较为广泛,但常会导致影像异质区域边缘模糊、纹理损失;(2)频率域滤波方法,主要包含离散小波变换和二维离散傅立叶变换下的低通、高通和高频滤波等方法。该类方法是将影像从空间域通过离散变换转换到频率域,在频率域进行降噪处理,再通过逆变换转回到空间域,从而实现降噪的目的。
小波变换属于典型的成像后滤波处理技术,其可以将时域和频域相结合来处理信号,克服了以往降噪方法的局限性,是目前成像后滤波处理研究的热点。小波变换最初的思想形成于20世纪初,法国的Morlet在1974年首次提出了小波变换的方法。1989年,Mallat将多尺度分析的思想引入到小波变换中,给出了构造正交小波基的一般方法和与FFT相对应的快速小波算法(Mallat算法),并将分散在各应用领域中的小波成果统一在同一个理论框架下。1995年,斯坦福大学的Donoho团队提出了小波域阈值滤波算法,取得了大量理论及应用成果,为降噪研究奠定了坚实基础。在小波域阈值滤波中,阈值的选取尤为重要,会直接影响到最终的滤波结果。目前,确定阈值的方法有模极大值法、平移不变量法、硬阈值法、软阈值法和半软阈值法等。模极大值法降噪后信噪比较高,但算法复杂,难于在实际处理中广泛应用;平移不变量法可以处理不连续点,但计算速度较慢;硬阈值法效果较好,但存在伪Gibbs现象,导致视觉失真;软阈值法滤波结果相对平滑,但会造成边缘模糊等失真现象;半软阈值法提高了信噪比,降噪效果得到了明显的提升,但需要确定两个阈值,增加了不可靠性。相关研究表明,阈值过大会造成有用信息的丢失,阈值过小会造成噪声去除不彻底。
近20年来,小波滤波方法不断得到改进与提升。Gao等人提出了一种基于平稳小波变换结合方向滤波器组的SAR图像去噪方法。该方法使用贝叶斯最大后验估计来替代传统阈值方法,降噪效果明显,但计算量太大。Zhang等人在马尔可夫随机场模型的基础上,结合改进后的粒子滤波器,提出了一种新的基于小波的SAR影像去斑方法,实验结果表明该方法有效地保留了SAR图像的边缘、纹理信息和结构特征,但算法过程较为繁琐。李金伦等人针对混合脉冲和高斯噪声,改进中值和提升小波变换去噪方法,实现了对混合噪声的去除。查显杰等人基于Daubechies小波对SAR干涉图进行小波包去噪,采用多级分解降噪,结果显示比直接使用软阈值法去噪效果更有效,同时证实了分层阈值的优势。Bhateja等人结合小波和各向异性扩散滤波器,考虑了不同滤波器的组合情况,在环境检测中有效地对SAR图像的斑点噪声进行了抑制。江虹等人通过一个变量μ结合软、硬阈值重构了小波阈值函数,降噪效果明显。Saravani等人提出了一种新的小波变换迭代算法,该算法通过聚合平稳小波变换、双边滤波、贝叶斯估计和各向异性扩散滤波对SAR影像降噪,有效降低了斑点噪声,保留了图像的边缘细节信息;但迭代算法计算量过大,可能造成低效率的问题。张玉叶等人利用平稳小波对信号进行分解,采用增强Lee滤波和自适应阈值滤波方法对不同信号分别处理,最终实现了降噪的目的;但结合两种滤波处理过程相对复杂,并没有对小波滤波方面进行细化和改进。Yan等人应用快速非局部均值滤波获取了小波包系数,对SAR影像实部和虚部分别进行小波包变换,避免了相位跳跃对后续滤波的影响,有效提高了降噪的质量;但增加了该算法的复杂性。
综上,SAR影像降噪作为雷达遥感领域长期的研究热点,以往传统滤波方法很难兼顾有效降噪的同时保持边缘细节信息和特征。近年来,小波降噪作为一种应用广泛的新型滤波方法,许多学者虽然使用不同研究思路提升了小波降噪效果,但大部分算法往往存在以下局限性:(1)算法较为复杂,导致计算量过大;(2)算法引入其它变量,增加了不可靠性;(3)在小波降噪方法中对分解、降噪和重构等细节归纳与研究不够。
发明内容
本发明的一个目的是解决至少上述问题和/或缺陷,并提供至少后面将说明的优点。
本发明还有一个目的是提供一种SAR干涉图小波降噪处理方法及其处理装置,其既能够有效降噪,又能够保持边缘细节信息。
为了实现本发明的这些目的和其它优点,第一方面,本发明提供了一种SAR干涉图小波降噪处理方法,包括:
步骤一、确定SAR干涉图中的信号点和噪声点,构建噪点识别矩阵,以便于仅对噪声点进行滤波降噪处理;
步骤二、提取出含有噪声的干涉图数据,分解为实部和虚部,利用Biorthogonal小波基对实部和虚部分别进行小波分解,并采用白化检验的方式使小波分解的层数达到最优;
步骤三、利用改进的小波阈值函数和Birge-Massart惩罚函数自适应获取的分层阈值对小波分解后得到的小波系数进行处理,即对实部和虚部分别进行滤波降噪;
步骤四、对滤波降噪后的实部和虚部进行小波重构,获取滤波降噪后的干涉图,并结合步骤一中SAR干涉图中的信号点,获取最终的降噪干涉图。
优选的是,所述的SAR干涉图小波降噪处理方法,对实部和虚部分别进行小波分解时采用的二维Mallat快速分解公式为:
其中,d为小波系数,c为尺度系数,V、H和D分别为垂直方向、水平方向和对角线方向,j为尺度,h(N)为低通滤波器系数,g(N)为高通滤波器系数,m、n为二维信号的序列,k=2m+N,l=2n+N,m,n,k,l∈Z;
对滤波降噪后的实部和虚部进行小波重构,采用的二维Mallat重构公式为:
优选的是,所述的SAR干涉图小波降噪处理方法,所述步骤三中,改进的小波阈值函数为:
其满足:
优选的是,所述的SAR干涉图小波降噪处理方法,所述步骤一具体包括:
计算SAR干涉图中的每一个点的相邻像素点的相位梯度;
判断相位梯度累计值是否为零;
若不为零,则将该点设为噪声点,反之,将该点设为信号点;
设置噪声点的值为1,信号点的值为0,构造元素为0/1的噪点识别矩阵。
优选的是,所述的SAR干涉图小波降噪处理方法,SAR干涉图中每一个点(x,y)的相邻像素点的相位梯度计算方式为:
Δ1=px,y-px+1,y
Δ2=px+1,y-px+1,y+1
Δ3=px+1,y+1-px,y+1
Δ4=px,y+1-px,y
优选的是,所述的SAR干涉图小波降噪处理方法,所述采用白化检验的方式使小波分解的层数达到最优的具体过程为:
S1、预设分解层数f为1;
S2、获取小波分解后的小波系数df,j;
S3、计算自相关性系数ρf,计算方式为:
S5、若小于等于,则f=f+1,继续执行S2、S3和S4;
S6、若不小于等于,则确定最优分解层数fbest=f-1。
优选的是,所述的SAR干涉图小波降噪处理方法,所述Birge-Massart惩罚函数的表达式为:
crit(r)=-sum(w(v)2,v≤r)+2sigma2r(ALPHA+log(M/r));
thr=|w(r*)|,
其中,w(v)为按照绝对值递减的顺序存储的小波系数,ALPHA为稀疏参数,M是小波系数的个数,r*为Birge-Massart惩罚函数的极小值,thr为分层阈值。
第二方面,本发明提供了一种SAR干涉图小波降噪处理装置,包括:
噪点识别矩阵构建模块,其用于确定SAR干涉图中的信号点和噪声点,构建噪点识别矩阵;
实部和虚部分解模块,其用于对提取出的含有噪声的干涉图数据进行实部和虚部分解;
小波分解模块,其利用Biorthogonal小波基对实部和虚部分别进行小波分解,并采用白化检验的方式使小波分解的层数达到最优;
小波系数处理模块,其利用改进的小波阈值函数和Birge-Massart惩罚函数自适应获取的分层阈值对小波分解后得到的小波系数进行处理;
小波重构模块,其用于对滤波降噪后的实部和虚部进行小波重构,获取滤波降噪后的干涉图。
第三方面,本发明提供了一种电子设备,包括:至少一个处理器,以及与所述至少一个处理器通信连接的存储器,其中,所述存储器存储有可被所述至少一个处理器执行的指令,所述指令被所述至少一个处理器执行,以使所述至少一个处理器执行所述的方法。
第四方面,本发明提供了一种存储介质,其上存储有计算机程序,该程序被处理器执行时,实现所述的方法。
本发明至少包括以下有益效果:一、本发明先对SAR干涉图中的每一个点进行分析判断,确定SAR干涉图中有哪些点是信号点,哪些点是噪声点,从而构建噪点识别矩阵,以便于对信号点不做滤波处理,仅对噪声点进行滤波降噪处理,避免了干涉图出现过平滑现象而导致边缘细节信息丢失。二、本发明对分解的实部和虚部进行小波分解时,采用的Biorthogonal小波基可以有效地消除斑点噪声,同时计算效率也高。三、本发明使用的改进后的阈值函数,在阈值t处连续,克服了硬阈值法的缺点,同时也解决了软阈值法存在的较大偏差的问题。还由于考虑到软硬阈值函数同时只有一个阈值需要确定,该方法比半软阈值法更加简单可靠。四、本发明采用白化检验方法自适应地获取最优分解层数,既能够达到降噪的目的,同时又能够极大程度上保留细节信息。五、在最优分层数的基础上,本发明使用Birge-Massart惩罚函数方法可实现自适应确定不同尺度下不同层的阈值。因此,本发明提供的方法,能够实现对噪声的有效去除,且经过模拟数据和真实数据的实验验证,也说明本发明提供的方法能够获取质量较好的降噪干涉图,并且优于实验中其它竞争滤波方法。
本发明的其它优点、目标和特征将部分通过下面的说明体现,部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。
附图说明
图1为本发明其中一个实施例中SAR干涉图小波降噪处理方法流程示意图。
图2为本发明其中一个实施例中构建噪点识别矩阵流程示意图。
图3为本发明其中一个实施例中采用白化检验的方式使小波分解的层数达到最优的流程示意图。
图4为本发明其中一个实施例中SAR干涉图小波降噪处理装置的关系示意图。
图5为模拟数据实验中的模拟无噪声干涉相位图。
图6为模拟实验s2=0.5时的加噪图;图7为模拟实验s2=0.5时的自适应均值滤波图;图8为模拟实验s2=0.5时的维纳滤波图;图9为模拟实验s2=0.5时的lee滤波图;图10为模拟实验s2=0.5时的frost滤波图;图11为模拟实验s2=0.5时的本发明实施例滤波方法图。
图12为模拟实验s2=0.8时的加噪图;图13为模拟实验s2=0.8时的自适应均值滤波图;图14为模拟实验s2=0.8时的维纳滤波图;图15为模拟实验s2=0.8时的lee滤波图;图16为模拟实验s2=0.8时的frost滤波图;图17为模拟实验s2=0.8时的本发明实施例滤波方法图。
图18为本发明其中一个实施例真实数据实验中蒙古地震的原图;图19为本发明其中一个实施例真实数据实验中蒙古地震的自适应均值滤波图;图20为本发明其中一个实施例真实数据实验中蒙古地震的维纳滤波图;图21为本发明其中一个实施例真实数据实验中蒙古地震的lee滤波图;图22为本发明其中一个实施例真实数据实验中蒙古地震的frost滤波图;图23为本发明其中一个实施例真实数据实验中蒙古地震的本发明实施例滤波方法图。
图24为本发明其中一个实施例真实数据实验中伽师地震的原图;图25为本发明其中一个实施例真实数据实验中伽师地震的自适应均值滤波图;图26为本发明其中一个实施例真实数据实验中伽师地震的维纳滤波图;图27为本发明其中一个实施例真实数据实验中伽师地震的lee滤波图;图28为本发明其中一个实施例真实数据实验中伽师地震的frost滤波图;图29为本发明其中一个实施例真实数据实验中伽师地震的本发明实施例滤波方法图。
图30为本发明其中一个实施例真实数据实验中陕西矿区的原图;图31为本发明其中一个实施例真实数据实验中陕西矿区的自适应均值滤波图;图32为本发明其中一个实施例真实数据实验中陕西矿区的维纳滤波图;图33为本发明其中一个实施例真实数据实验中陕西矿区的lee滤波图;图34为本发明其中一个实施例真实数据实验中陕西矿区的frost滤波图;图35为本发明其中一个实施例真实数据实验中陕西矿区的本发明实施例滤波方法图。
图36为本发明其中一个实施例真实数据实验中唐山矿区的原图;图37为本发明其中一个实施例真实数据实验中唐山矿区的自适应均值滤波图;图38为本发明其中一个实施例真实数据实验中唐山矿区的维纳滤波图;图39为本发明其中一个实施例真实数据实验中唐山矿区的lee滤波图;图40为本发明其中一个实施例真实数据实验中唐山矿区的frost滤波图;图41为本发明其中一个实施例真实数据实验中唐山矿区的本发明实施例滤波方法图。
具体实施方式
下面结合附图对本发明做进一步的详细说明,以令本领域技术人员参照说明书文字能够据以实施。
如图1所示,本发明实施例提供了一种SAR干涉图小波降噪处理方法,包括下列步骤:
S10、确定SAR干涉图中的信号点和噪声点,构建噪点识别矩阵,以便于仅对噪声点进行滤波降噪处理。
其中,如图2所示,构建噪点识别矩阵,进行噪点先验的具体过程为:
S11、计算SAR干涉图中的每一个点(x,y)的相邻像素点的相位梯度,
S13、若不为零,则将该点设为噪声点,反之,将该点设为信号点;
S14、设置噪声点的值为1,信号点的值为0,构造元素为0/1的噪点识别矩阵,构造的噪点识别矩阵与研究对象大小一样;
S15、保留信号点不做滤波处理,仅对噪点进行滤波处理。
在上述步骤中,由于斑点噪声在空间域中表现为较强的不连续性,为了避免干涉图出现过平滑现象而导致边缘细节信息丢失,因此采用了对噪点进行先验处理的方法。
S20、提取出含有噪声的干涉图数据,分解为实部和虚部,利用Biorthogonal小波基对实部和虚部分别进行小波分解,并采用白化检验的方式使小波分解的层数达到最优。
其中,SAR干涉图是复数图像,包含幅度和相位双重信息,降噪是基于复数的降噪。Lee等人在1998年证明了干涉图中噪声满足加性噪声模型,模型为:ψz=ψx+υ,式中,Ψz为真实测量的干涉图,Ψx为无噪声的干涉图,υ为零均值的噪声。Ψx与υ被假定为独立的随机变量。加性噪声的特点是噪声方差与均值无关。在加性噪声模型的基础上,2002年López-Martínez和Fàbregas进一步证明了干涉图的实部和虚部中噪声满足该模型,并定义了干涉相位噪声的模型。实部和虚部的表达式分别为:
cos(ψz)=Nc cos(ψx)+υc
υc=υ1′cos(ψx)-υ2′sin(ψx)
sin(ψz)=Nc sin(ψx)+υs
υs=υ1′sin(ψx)+υ2′cos(ψx)
式中,Nc为干涉相位质量参数,υc为复干涉图实部零均值的噪声,υs为干涉图虚部零均值的噪声,υ1’和υ2’为零均值的随机变量。因此,噪声具有由Nc给出的乘法行为以及由υc和υs给出的加法行为,但Nc与υc、υs不是独立的,取决于其相关性ρ。随着ρ的增加,Nc增加,υc和υs的方差减少。因此,在对SAR干涉图进行降噪时,需要将其分解为实部和虚部,分别进行小波变换进行小波降噪处理。
其中,在进行小波变换时,小波基的选取是小波变换过程中的非常关键的步骤。因此,需要对小波基的数学特性进行对比分析。主要评价标准包含:(1)是否具有紧支性和适宜的支撑长度,保证小波函数具有良好时频局部特征;(2)是否具有正则性,保证小波函数可微性和重构的稳定性;(3)是否具有对称性,控制重构时发生的畸变;(4)是否具有适当的消失矩,利于奇异监测;(5)是否具有线性相位,防止小波分解和重构中发生失真情况。七种常用的小波基性质对比如下表一所示:
表一 小波基的性质
由表一可知,以上七种常用小波基在正交性、紧支性、对称性、正则性和消失矩方面的性质比较相似,但在线性相位特性上,Biorhogonal小波基具有相对优势。本发明实施例采用Biorhogonal小波基可以有效地消除斑点噪声,同时计算效率较高。
本发明实施例对实部和虚部分别进行小波变换分解时采用的二维Mallat快速分解公式为:
其中,d为小波系数,c为尺度系数,V、H和D分别为垂直方向、水平方向和对角线方向,j为尺度,h(N)为低通滤波器系数,g(N)为高通滤波器系数,m、n为二维信号的序列,k=2m+N,l=2n+N,m,n,k,l∈Z。
通过小波变换,含有噪声的影像可得到尺度信息,因为信号和噪声的小波系数在不同尺度上具有不同性质,所以尺度信息可作为相关部分,降噪是将相关部分进行阈值处理,即序列大于设定阈值时置零,小于设定阈值时不变,这一过程称为小波去相关。若小波分解层数过少,尺度信息占主导地位,噪声去除不明显;而分解层数过多会造成细节信息丢失。因此,本发明实施例采用白化验证的方式使小波分解的层数达到最优。
其中,如图3所示,所述采用白化检验的方式使小波分解的层数达到最优的具体过程为:
S21、预设分解层数f为1;
S22、获取小波分解后的小波系数df,j;
S23、计算自相关性系数ρf,计算方式为:
S25、若小于等于,则f=f+1,继续执行S22、S23和S24;
S26、若不小于等于,则确定最优分解层数fbest=f-1。
对小波分解后的小波系数进行白化检验,若满足白噪声检验条件,说明该分量表现为白噪声特性,否则表现为非白噪声特性,小波分解层数达到最佳。因此,采用白化检验方法自适应地获取最优分解层数,可达到降噪目的,同时极大程度上保留细节信息。
S30、利用改进的小波阈值函数和Birge-Massart惩罚函数自适应获取的分层阈值对小波分解后得到的小波系数进行处理,即对实部和虚部分别进行滤波降噪。
其中,在对小波系数进行阈值降噪处理过程中,小波阈值函数体现了对不同小波系数的不同处理策略,软阈值函数和硬阈值函数是Donoho提出的两种阈值处理函数。采用硬阈值函数降噪较好地保留影像的边缘特征,但降噪结果方差较大,会出现由于信号不连续点位置导致的伪Gibbs现象,造成视觉失真;采用软阈值函数降噪的效果比硬阈值法会更加平滑,但可能会由于过度平滑造成边缘信息丢失,导致失真现象;硬阈值法主要因为不连续性造成结果有较大方差,软阈值法主要是因为对所有大于阈值的系数共同做了收缩造成结果有较大偏差,基于此,半软阈值法集合上述两种方法的优势,克服了其局限性,但该方法需要确定两个阈值,增加了算法的复杂性,使不可靠因素增多。本发明实施例采用改进的小波阈值函数对小波系数进行处理。改进的小波阈值函数为:
其满足:
其中,为施加阈值改进后的小波系数,t为阈值,sgn(·)为符号函数。改进的小波阈值函数在阈值t处连续,克服了硬阈值法的缺点;同时,以θj(m,n)=dj(m,n)为渐近线,解决了软阈值法存在较大偏差的问题;考虑到软硬阈值函数同时只有一个阈值需要确定,因此,该方法比半软阈值法更加简单可靠。
为了获得不同尺度下不同层的分层阈值,在最优分层的基础上,本发明实施例使用Birge-Massart惩罚函数方法可实现自适应确定不同尺度下不同层阈值的目的。所述Birge-Massart惩罚函数的表达式为:
crit(r)=-sum(w(v)2,v≤r)+2sigma2r(ALPHA+log(M/r));
thr=|w(r*)|,
其中,w(v)为按照绝对值递减的顺序存储的小波系数,ALPHA为稀疏参数,M是小波系数的个数,r*为Birge-Massart惩罚函数的极小值,thr为分层阈值。因此,获得分层阈值后,利用改进的小波阈值函数可以对不同尺度下不同层的小波系数进行滤波降噪处理,此时,改进的小波阈值函数中的阈值t即为分层阈值thr。当该层的小波系数的绝对值小于等于该层的分层阈值时,令其为零;小波系数的绝对值大于该层的分层阈值时,根据改进的小波阈值函数的规则做相应变换。
S40、对滤波降噪后的实部和虚部进行小波重构,获取滤波降噪后的干涉图,并结合步骤一中SAR干涉图中的信号点,获取最终的降噪干涉图。
其中,对滤波降噪后的实部和虚部进行小波重构,即对滤波降噪后的实部和虚部进行小波逆变换,将实部和虚部结合后,构成滤波后的干涉图。然后再结合步骤S10中SAR干涉图的信号点构成最终的降噪干涉图。进行小波重构时采用二维Mallat重构算法公式:
本发明又一实施例提供了一种SAR干涉图小波降噪处理装置,如图4所示,包括:
噪点识别矩阵构建模块,其用于确定SAR干涉图中的信号点和噪声点,构建噪点识别矩阵;
实部和虚部分解模块,其用于对提取出的含有噪声的干涉图数据进行实部和虚部分解;
小波分解模块,其利用Biorthogonal小波基对实部和虚部分别进行小波分解,并采用白化检验的方式使小波分解的层数达到最优;
小波系数处理模块,其利用改进的小波阈值函数和Birge-Massart惩罚函数自适应获取的分层阈值对小波分解后得到的小波系数进行处理;
小波重构模块,其用于对滤波降噪后的实部和虚部进行小波重构,获取滤波降噪后的干涉图。
本发明其中一实施例还提供了一种电子设备,包括:至少一个处理器,以及与所述至少一个处理器通信连接的存储器,其中,所述存储器存储有可被所述至少一个处理器执行的指令,所述指令被所述至少一个处理器执行,以使所述至少一个处理器执行所述的方法。
本发明其中一实施例还提供了一种存储介质,其上存储有计算机程序,该程序被处理器执行时,实现所述的方法。
如上所述,本发明实施例提供的SAR干涉图小波降噪处理方法及其处理装置,既能够有效降噪,又能够保持边缘细节信息。
下面将结合具体的实施例来验证本发明实施例所述的SAR干涉图小波降噪处理方法的有效性和优化性,分别进行模拟数据实验和真实数据实验。为了进行算法比较,使用竞争滤波方法,包括自适应均值滤波、维纳滤波、lee滤波和frost滤波。为了评价降噪效果,分别从定性和定量两方面展开评价,定性评价采取目视解译方法,定量评价采取残差点数量、PSNR(峰值信噪比)、RMSE(均方根误差)和EPI(边缘保持指数)指标。
一、模拟数据实验
实验采用模拟大小为500*500的无噪声干涉相位,如图5所示。
为了创建噪声图像,在原图的基础上,分别添加方差为0.5和0.8的加性斑点噪声。通过竞争滤波方法和本发明实施例滤波方法分别进行对含有噪声的模拟干涉相位图进行降噪处理。结果如图6加噪图、图7自适应均值滤波、图8维纳滤波、图9lee滤波、图10frost滤波、图11本发明实施例方法,图12加噪图、图13自适应均值滤波、图14维纳滤波、图15lee滤波、图16frost滤波、图17本发明实施例方法和表二、表三所示,其中,图6、图7、图8、图9、图10、图11和表2是方差s2=0.5噪声水平下的结果,图12、图13、图14、图15、图16、图17和表3是方差s2=0.8噪声水平下的结果。
表二 模拟实验(s2=0.5)的定量指标结果表
由s2=0.5噪声水平下的实验结果可知,各滤波方法都在不同程度上滤除了噪声。但是经过目视解译,可直观地观察到,各竞争滤波方法的降噪结果仍存在较多颗粒感的噪点,尤其是自适应均值滤波结果失真严重。本发明实施例采用的方法的结果更具平滑性和连续性。通过表2的各定量指标对比分析可知,本发明实施例滤波方法与无噪声原图更具相似性,降噪后干涉图质量较好,保真性最佳,同时在边缘信息的保持能力也最优。
表三 模拟实验(s2=0.8)的定量指标结果表
s2=0.8噪声水平下的实验结果与s2=0.5噪声水平下的实验结果近似一致。在s2=0.5噪声水平下本发明实施例方法比竞争方法中最优方法的指标分别在残差点数量、PSNR、RMSE、EPI指标上优化了79.95%、12.51%、51.50%和9.99%,而在s2=0.8噪声水平下本发明实施例方法比竞争方法中最优方法的指标分别在残差点数量、PSNR、RMSE和EPI指标上优化了84.23%、13.49%、52.54%和17.72%,实验表明在噪声水平越大的情况下,本发明实施例提出的滤波方法的优势越明显。
总体而言,本发明实施例提出的滤波算法较上述传统竞争方法在一定程度上有所提升,且在噪声水平越高的情况下,本发明实施例滤波算法的优化性越强。
二、真实数据实验
为了进一步验证该算法的有效性,本发明实施例选取两次地震数据和两次矿区数据,共进行四次真实数据的降噪实验,从多角度证明本发明实施例算法的有效性。
实验一:蒙古地震
实验一采用2021年1月12日蒙古地震的干涉图(2000*2000)进行滤波处理,结果如图18原图、图19自适应均值滤波、图20维纳滤波、图21lee滤波、图22frost滤波、图23本发明实施例滤波方法和表四所示。
表四 真实数据实验一的定量指标结果表
本次实验原始干涉图噪声分布广而多,左下部分属于水域,整体干涉条纹明显。由表四可知,在4种竞争方法中,lee滤波结果残差点数量最少,这是由于lee滤波针对乘性噪声基于标准差进行滤波,滤波质量较高;维纳滤波PSNR和RMSE指标最优,是由于维纳滤波的本质遵循最小均方差的原则,但在保持边缘信息同时可能会产生伪边缘的错误信息;frost滤波边缘保持能力最强,这是由于frost滤波基于阻尼指数,一定程度能够较好地保持影像的细节特征和信息。而本发明实施例提出的滤波方法无论是在滤波后干涉图质量方面,还是降噪效果和保持边缘信息方面都表现出更优水平,且在残差点数量、PSNR、RMSE和EPI指标上都比竞争性方法分别优化了20.98%、5.96%、15.72%和17.58%。
实验二、伽师地震
实验二采用2020年1月19日伽师地震的干涉图(1000*1500)进行滤波处理,结果如图24原图、图25自适应均值滤波、图26维纳滤波、图27lee滤波、图28frost滤波、图29本发明实施例滤波方法和表五所示。
表五 真实数据实验二的定量指标结果表
本次实验原始干涉图噪声分布较少,无水域部分,整体干涉条纹相对更加明显和连续。在残差点数量中,本发明实施例滤波方法并没有达到最优,仅次于lee滤波,结合其他三个指标进一步分析,得出lee滤波可能在此次真实数据噪声水平下造成过度降噪,导致残差点数量最少和边缘信息保持能力过低。最终,由滤波结果和多层次分析可知,本次实验与实验一的结论相似,且本发明实施例滤波方法比竞争方法中的最优方法在PSNR、RMSE和EPI指标上分别优化了4.76%、23.21%和14.09%。结合实验一结果,在噪声水平越大的情况下,本发明实施例滤波方法的优势越明显。这个结论与模拟实验的结论一致。
实验三、陕西矿区
实验三采用陕西柠条塔南翼矿区2019年12月15日至2019年12月27日的干涉图(300*300)进行滤波处理。矿区范围相比地震范围偏小,所以选取以该实验数据的矿区为中心,扩大研究范围取300*300大小的干涉相位图作为研究对象。结果如图30原图、图31自适应均值滤波、图32维纳滤波、图33lee滤波、图34frost滤波、图35本发明实施例滤波方法和表六所示。
表六 真实数据实验三的定量指标结果表
本次实验数据噪声分布广而多,无水域部分,该干涉图范围内分布着不同大小的矿区,矿区周边基岩裸露,危崖陡峭,整体干涉条纹不明显。由如图30、图31、图32、图33、图34、图35目视解译可知,维纳滤波存在过度滤波现象,由表六指标可知,自适应均值滤波和lee滤波也存在过平滑情况。单一追求平滑导致干涉图模糊,边缘特征和信息丢失。在残差点数量指标下,各类竞争方法与本文滤波方法分别优化了75.93%、63.18%、74.50%、55.25%和75.15%。虽然本发明实施例滤波方法在残差点数量指标中并非最佳并略次于自适应均值滤波,但在降噪过程中还需要较高程度保持边缘细节信息的综合考虑下,本发明实施例滤波方法优于其余竞争方法,且本发明实施例滤波方法比竞争方法中的最优方法在PSNR、RMSE和EPI指标上分别优化了2.87%、16.08%和60.74%。本次实验与实验一对比分析,在不同干涉条纹表现的相位图中,各指标优化程度不尽相同。目视观察结果表明,干涉条纹明显的区域降噪效果更优;定量指标结果表明,干涉条纹不明显区域边缘保持能力优化更佳。
实验四、唐山矿区
实验四采用唐山矿区2016年01月02日至2016年01月14日的干涉图(300*300)进行滤波处理。结果如图36原图、图37自适应均值滤波、图38维纳滤波、图39lee滤波、图40frost滤波、图41本发明实施例滤波方法和表七所示。
表七 真实数据实验四的定量指标结果表
本次实验数据噪声分布广而多,该干涉图范围内分布着不同大小的矿区,矿区周边有水域部分,整体干涉条纹不明显,单个矿区范围内干涉条纹相对清晰。由图36、图37、图38、图39、图40、图41目视解译可知,维纳滤波、lee滤波和frost滤波存在明显相位值不连续的马赛克效应[36],自适应均值滤波比本文滤波方法的结果更平滑一些,但由表七指标可知,自适应均值滤波边缘保持指数较低,说明存在过平滑现象。综合而言,本发明实施例滤波方法优于其余竞争方法,且本发明实施例滤波方法比竞争方法中的最优方法在残差点数量、PSNR、RMSE和EPI指标上分别优化了33.06%、1.54%、8.50%和17.29%。
通过上述实验结果的对比分析,证明了本发明实施例提出的滤波算法的有效性和优越性,该算法在有效去除噪声的同时保持了干涉图的边缘细节信息,并且在噪声水平越高的情况下,滤波效果越理想。
综上,在本发明实施例中,提出了基于噪点先验的小波降噪改进处理方法,应用于SAR影像干涉图去噪研究。由于单一小波降噪算法可能造成影像边缘信息丢失,为了高效去除噪声,先详细地对SAR干涉图中的噪声进行了分析,结合噪点先验方法对影像信号进行判断,然后再进行对应的滤波处理,以保证影像边缘信息的完整。同时,将传统小波域阈值降噪分为三步进行优化,分别为小波基选取、传统阈值函数改进和基于白化检验的自适应分层阈值确定。在小波基选取时,通过对正交性、紧支性、对称性、正则性、消失矩和线性相位特性上进行比较分析,从而选取最合适的小波基;在改进阈值函数时,结合硬、软和半软阈值法优缺点进行选取;在分层阈值确定时,引入白化检验获取最优分层,借助惩罚函数自适应获取阈值。在此基础上,构建小波函数中的规则模型进行降噪处理,实现对噪声的有效去除。经过模拟数据和真实数据的实验验证,表明本发明实施例提出的方法能够获取质量较好的降噪干涉图,并且优于实验中其他竞争滤波方法,为后期SAR影像反演形变信息等研究提供了更为可靠的干涉图。同时,该方法在越高的噪声水平中,滤波效果越理想。
尽管本发明的实施方案已公开如上,但其并不仅仅限于说明书和实施方式中所列运用。它完全可以被适用于各种适合本发明的领域。对于熟悉本领域的人员而言,可容易地实现另外的修改。因此在不背离权利要求及等同范围所限定的一般概念下,本发明并不限于特定的细节和这里示出与描述的图例。
Claims (9)
1.SAR干涉图小波降噪处理方法,其特征在于,包括:
步骤一、确定SAR干涉图中的信号点和噪声点,构建噪点识别矩阵,以便于仅对噪声点进行滤波降噪处理;
步骤二、提取出含有噪声的干涉图数据,分解为实部和虚部,利用Biorthogonal小波基对实部和虚部分别进行小波分解,并采用白化检验的方式使小波分解的层数达到最优;
步骤三、利用改进的小波阈值函数和Birge-Massart惩罚函数自适应获取的分层阈值对小波分解后得到的小波系数进行处理,即对实部和虚部分别进行滤波降噪;
步骤四、对滤波降噪后的实部和虚部进行小波重构,获取滤波降噪后的干涉图,并结合步骤一中SAR干涉图中的信号点,获取最终的降噪干涉图;
其中,所述步骤三中,改进的小波阈值函数为:
其满足:
3.如权利要求1所述的SAR干涉图小波降噪处理方法,其特征在于,所述步骤一具体包括:
计算SAR干涉图中的每一个点的相邻像素点的相位梯度;
判断相位梯度累计值是否为零;
若不为零,则将该点设为噪声点,反之,将该点设为信号点;
设置噪声点的值为1,信号点的值为0,构造元素为0/1的噪点识别矩阵。
6.如权利要求1所述的SAR干涉图小波降噪处理方法,其特征在于,所述Birge-Massart惩罚函数的表达式为:
crit(r)=-sum(w(v)2,v≤r)+2sigma2r(ALPHA+log(M/r));
thr=|w(r*)|,
其中,w(v)为按照绝对值递减的顺序存储的小波系数,ALPHA为稀疏参数,M是小波系数的个数,r*为Birge-Massart惩罚函数的极小值,thr为分层阈值。
7.用于处理权利要求1~6任一所述的SAR干涉图小波降噪处理方法的SAR干涉图小波降噪处理装置,其特征在于,包括:
噪点识别矩阵构建模块,其用于确定SAR干涉图中的信号点和噪声点,构建噪点识别矩阵;
实部和虚部分解模块,其用于对提取出的含有噪声的干涉图数据进行实部和虚部分解;
小波分解模块,其利用Biorthogonal小波基对实部和虚部分别进行小波分解,并采用白化检验的方式使小波分解的层数达到最优;
小波系数处理模块,其利用改进的小波阈值函数和Birge-Massart惩罚函数自适应获取的分层阈值对小波分解后得到的小波系数进行处理;
小波重构模块,其用于对滤波降噪后的实部和虚部进行小波重构,获取滤波降噪后的干涉图。
8.一种电子设备,其特征在于,包括:至少一个处理器,以及与所述至少一个处理器通信连接的存储器,其中,所述存储器存储有可被所述至少一个处理器执行的指令,所述指令被所述至少一个处理器执行,以使所述至少一个处理器执行权利要求1-6中任一项所述的方法。
9.一种存储介质,其上存储有计算机程序,其特征在于,该程序被处理器执行时,实现权利要求1-6中任一项所述的方法。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210020853.6A CN114355348B (zh) | 2022-01-10 | 2022-01-10 | Sar干涉图小波降噪处理方法及其处理装置 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202210020853.6A CN114355348B (zh) | 2022-01-10 | 2022-01-10 | Sar干涉图小波降噪处理方法及其处理装置 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN114355348A CN114355348A (zh) | 2022-04-15 |
CN114355348B true CN114355348B (zh) | 2022-11-08 |
Family
ID=81107856
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202210020853.6A Active CN114355348B (zh) | 2022-01-10 | 2022-01-10 | Sar干涉图小波降噪处理方法及其处理装置 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN114355348B (zh) |
Families Citing this family (2)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN116068521B (zh) * | 2023-03-15 | 2023-06-23 | 长沙东玛克信息科技有限公司 | 一种雷达探测信号主动降噪方法 |
CN117473233B (zh) * | 2023-12-27 | 2024-03-01 | 华东交通大学 | 样品成分分析方法、***、存储介质及计算机 |
Family Cites Families (10)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US7260272B2 (en) * | 2003-07-10 | 2007-08-21 | Samsung Electronics Co.. Ltd. | Method and apparatus for noise reduction using discrete wavelet transform |
CN102928517A (zh) * | 2012-11-15 | 2013-02-13 | 河北省电力公司电力科学研究院 | 一种基于小波分解阈值去噪的瓷绝缘子振动声学检测数据降噪的方法 |
CN103761719B (zh) * | 2014-01-06 | 2017-01-04 | 暨南大学 | 一种基于邻域相关性的自适应小波阈值去噪方法 |
CN103854264B (zh) * | 2014-03-28 | 2016-05-11 | 中国石油大学(华东) | 一种基于改进型阈值函数的小波变换图像去噪方法 |
CN104156918B (zh) * | 2014-08-01 | 2017-02-15 | 西安电子科技大学 | 基于联合稀疏表示与残差融合的sar图像噪声抑制方法 |
CN104200440A (zh) * | 2014-09-16 | 2014-12-10 | 哈尔滨恒誉名翔科技有限公司 | 一种基于多尺度小波变换的斑点图像处理算法 |
CN104459633B (zh) * | 2014-12-01 | 2016-08-17 | 中国科学院电子学研究所 | 结合局部频率估计的小波域InSAR干涉相位滤波方法 |
CN105913393B (zh) * | 2016-04-08 | 2018-08-17 | 暨南大学 | 一种自适应小波阈值图像去噪方法及装置 |
CN107895354A (zh) * | 2017-11-29 | 2018-04-10 | 陕西师范大学 | 一种基于非下采样Shearlet变换的高分三SAR图像降斑方法 |
CN113378661B (zh) * | 2021-05-25 | 2024-03-22 | 浙江工业大学 | 一种基于改进小波阈值和相关检测的直流电能信号去噪方法 |
-
2022
- 2022-01-10 CN CN202210020853.6A patent/CN114355348B/zh active Active
Also Published As
Publication number | Publication date |
---|---|
CN114355348A (zh) | 2022-04-15 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN114355348B (zh) | Sar干涉图小波降噪处理方法及其处理装置 | |
Yu et al. | Image denoising using trivariate shrinkage filter in the wavelet domain and joint bilateral filter in the spatial domain | |
Corner et al. | Noise estimation in remote sensing imagery using data masking | |
CN101847257B (zh) | 基于非局部均值与多级定向图像的图像降噪方法 | |
CN105205788B (zh) | 一种针对高通量基因测序图像的去噪方法 | |
Ranjbarzadeh et al. | LNPSS: SAR image despeckling based on local and non-local features using patch shape selection and edges linking | |
Bhateja et al. | An improved local statistics filter for denoising of SAR images | |
Liu et al. | Synthetic aperture radar image de-noising based on Shearlet transform using the context-based model | |
Li et al. | Research on wavelet-based contourlet transform algorithm for adaptive optics image denoising | |
Li et al. | A method to improve the accuracy of SAR image change detection by using an image enhancement method | |
CN107085826B (zh) | 基于加权重叠非局部回归先验的图像超分辨率重建方法 | |
CN115082336A (zh) | 一种基于机器学习的sar图像相干斑抑制方法 | |
CN105976340A (zh) | 改进的基于小波分解的旋滤波算法 | |
Joseph et al. | Hybrid spatio-frequency domain global thresholding filter (HSFGTF) model for SAR image enhancement | |
CN111461999B (zh) | 一种基于超像素相似性测量的sar图像相干斑抑制方法 | |
Makandar et al. | Computation pre-processing techniques for image restoration | |
Gir et al. | Speckle reduction of synthetic aperture radar images using median filter and savitzky-golay filter | |
Abu-Ein | A novel methodology for digital removal of periodic noise using 2D fast Fourier transforms | |
Rahim et al. | Efficient Contourlet Transformation Technique for Despeckling of Polarimetric Synthetic Aperture Radar Image | |
Kaur et al. | Image de-noising techniques: a review paper | |
CN113506212B (zh) | 一种改进的基于pocs的高光谱图像超分辨率重建方法 | |
Fang et al. | De-noising of SAR images based on Wavelet-Contourlet domain and PCA | |
Bhargava et al. | An Effective Method for Image Denoising Using Non-local Means and Statistics based Guided Filter in Nonsubsampled Contourlet Domain. | |
Nisha et al. | Wavelet coefficients thresholding techniques for denoising MRI images | |
Kaur et al. | Study of Image Denoising and Its Techniques |
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 |