CN103645476A - 一种合成孔径雷达差分干涉图序列的时空同质滤波方法 - Google Patents
一种合成孔径雷达差分干涉图序列的时空同质滤波方法 Download PDFInfo
- Publication number
- CN103645476A CN103645476A CN201310699869.5A CN201310699869A CN103645476A CN 103645476 A CN103645476 A CN 103645476A CN 201310699869 A CN201310699869 A CN 201310699869A CN 103645476 A CN103645476 A CN 103645476A
- Authority
- CN
- China
- Prior art keywords
- pixel
- intensity
- radar
- image
- differential
- 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.)
- Granted
Links
- 238000000034 method Methods 0.000 title claims abstract description 19
- 238000001914 filtration Methods 0.000 claims abstract description 50
- 238000009826 distribution Methods 0.000 claims abstract description 38
- 230000001186 cumulative effect Effects 0.000 claims abstract description 35
- 238000005315 distribution function Methods 0.000 claims abstract description 32
- 238000004364 calculation method Methods 0.000 claims abstract description 5
- 238000005305 interferometry Methods 0.000 claims description 74
- 241001347978 Major minor Species 0.000 claims description 19
- 238000006073 displacement reaction Methods 0.000 claims description 7
- 239000004576 sand Substances 0.000 claims description 6
- 238000012952 Resampling Methods 0.000 claims description 4
- 230000015572 biosynthetic process Effects 0.000 claims description 3
- 238000003384 imaging method Methods 0.000 claims description 3
- 230000001788 irregular Effects 0.000 claims description 3
- 238000004088 simulation Methods 0.000 claims description 3
- 238000000528 statistical test Methods 0.000 claims description 3
- 238000006467 substitution reaction Methods 0.000 claims description 2
- 239000000654 additive Substances 0.000 abstract description 5
- 230000000996 additive effect Effects 0.000 abstract description 5
- 238000012300 Sequence Analysis Methods 0.000 abstract 1
- 230000003044 adaptive effect Effects 0.000 abstract 1
- 238000005259 measurement Methods 0.000 abstract 1
- 238000005070 sampling Methods 0.000 abstract 1
- 238000010586 diagram Methods 0.000 description 6
- 238000012731 temporal analysis Methods 0.000 description 4
- 238000000700 time series analysis Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 3
- 238000005516 engineering process Methods 0.000 description 2
- 230000002452 interceptive effect Effects 0.000 description 2
- 238000012544 monitoring process Methods 0.000 description 2
- 238000004064 recycling Methods 0.000 description 2
- 230000000717 retained effect Effects 0.000 description 2
- 240000007594 Oryza sativa Species 0.000 description 1
- 235000007164 Oryza sativa Nutrition 0.000 description 1
- 238000009825 accumulation Methods 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 230000007935 neutral effect Effects 0.000 description 1
- 238000012545 processing Methods 0.000 description 1
- 235000009566 rice Nutrition 0.000 description 1
- 230000009897 systematic effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
- G01S13/88—Radar or analogous systems specially adapted for specific applications
- G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
- G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
- G01S13/9021—SAR image post-processing techniques
- G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
- G01S7/00—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
- G01S7/02—Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
- G01S7/023—Interference mitigation, e.g. reducing or avoiding non-intentional interference with other HF-transmitters, base station transmitters for mobile communication or other radar systems, e.g. using electro-magnetic interference [EMI] reduction techniques
-
- 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
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- General Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Theoretical Computer Science (AREA)
- Image Processing (AREA)
Abstract
一种合成孔径雷达差分干涉图序列的时空同质滤波方法,它有八大步骤:步骤一:主副图像的自适应精密配准及副图像的重采样;步骤二:生成每张雷达图像的强度图;步骤三:生成序列差分干涉图;步骤四:确定滤波搜索窗口和统计邻域窗口的大小;步骤五:计算差分干涉图中每个像元的强度累积分布函数;步骤六:确定差分干涉图中每个像元在其滤波搜索窗口内的强度分布意义下的同质像元;步骤七:相位同质的度量;步骤八:差分干涉图的空间域滤波。本发明能有效地同时滤除差分干涉图中的乘性斑点噪声和加性白噪声,并能很好地保留雷达图像的原始分辨率,从而可以大大地提高差分干涉图时序分析的有效性和精度。
Description
技术领域
本发明涉及一种合成孔径雷达差分干涉图序列的时空同质滤波方法,属于合成孔径雷达干涉技术领域。它可以有效地同时滤除差分干涉图中的乘性斑点噪声和加性白噪声,并能很好地保留雷达图像的原始分辨率,从而可以大大地提高差分干涉图时序分析的有效性和精度。
背景技术
同所有的其它图像数据一样,合成孔径雷达差分干涉图也会受到各种噪声的污染。为了消除噪声,从被噪声污染了的图像中识别出真实图像,需要用滤波器对图像进行滤波。人们常用低通滤波器对图像进行滤波,在时域里,其实就是事先定义一个适当大小的小窗口(boxcar),让小窗口在图像上滑动,对这个小窗口中所有像元的数值进行平均运算,将平均值赋予该窗口所在位置的中心点像元。用数学表达式来表示的话,可以写成:
上式中,表示在窗口内所有(N个)像元(值)组成的N维向量(u1,u2,……,uN)T的均值。显然,这种滤波器只适用于常数图像的加性噪声的滤波,其期望值就等于其均值,也即图像常数。但在实际情况下,图像都是非常数的。这时,方差可以被减小了N倍:
上式中,用随机变量Uk表示窗口内序号为k的像元,Var表示方差。
方差虽然被减小了N倍,但均值却是有偏的,且可以是非常大:
这种有偏估计显然会极大地破坏图像的原始分辨率。
雷达差分干涉图可以表示为:
上式中,Im和Is分别表示主副图像的雷达反射强度(Intensity),φ表示主副图像去除参考面及高程相位后的相位差,即差分干涉相位。如果把雷达反射强度和差分干涉相位都看成是随机变量V和φ的话,那么雷达图像的强度将服从伽玛分布(GammaDistribution):
方差为:
显然,雷达强度信号的噪声是一种与信号相关的乘性噪声,要去除这类噪声,线性滤波器是无能为力的。
雷达差分干涉图的相位服从高斯分布:
为了将雷达差分干涉图中由强度乘性噪声,雷达专业中称之为斑点噪声(speckle),和相位高斯噪声消除,本发明设计出了一种基于噪声特性的滤波方法。在雷达差分干涉技术的实际运用中,常常是需要对一序列差分干涉图进行时序分析,因此,本发明利用分布在时间域上的序列雷达图像的强度信息去除斑点噪声,再利用空间域上的相位信息去除高斯噪声,故称之为雷达差分干涉图的时空同质滤波。所谓图像滤波,就是对图像中每一个像元要用一个新的值去替代它原来的值,这个新的值来源于对以该像元为中心的一个适度大小的窗口内所有其它像元的值的某种运算的结果。这个窗口叫滤波搜索窗口,如图(1)中粗黑线框所示。对搜索窗口中心像元s滤波,即是用搜索窗口内所有像元的加权平均来替代其原来的值:
发明内容
1.目的:本发明的目的是提供一种合成孔径雷达差分干涉图序列的时空同质滤波方法。它克服了现有滤波技术的不足,能在保持住原始雷达图像分辨率的情况下,滤除雷达斑点噪声和差分干涉相位中的高斯噪声,从而可获取高分辨率的雷达差分干涉相位,可以大大提高后续差分相位时序分析的有效性和精度
2.技术方案:本发明是一种合成孔径雷达差分干涉图序列的时空同质滤波方法,该方法具体步骤如下:
步骤一:主副图像的自适应精密配准及副图像的重采样
在主副图像进行配准时,不仅要考虑雷达干涉系统(成像几何,干涉几何,雷达运行轨道等等)引起的地面同一目标在主副图像上的位置的偏移,还要考虑地面目标因较大的高程或较大的形变位移引起的主副图像上的位置偏移。在本发明中,副图像利用Sinc函数进行重采样,Sinc函数的位移量由下式决定:
Offset(i,j)=(△i,△j)=((△i,△j)i,j)System+((△i,△j)i,j)Object (10)上式中,offset(i,j)表示了副图像中位于第i行第j列的像元相对于主图像上同一像元位置的偏移量,它是系统偏移和目标本身因高程或形变引起的偏移量的和。offset(i,j)由副图像s中位于(i,j)的像元的邻域与相应的主图像p中像元的邻域的协方差函数取最大值时的位移量决定,见下面的(11)和(12)式
上两式中,p和s分别表示主副图像,c表示在主副图像中以第i行第j列为中心的N行M列的两个窗口的协方差函数,n和m分别表示行和列偏移量,k和l表示窗口中参与相关计算像元的行列位置,(△i,△j)表示当c取得最大值时的行列号,即(nmax,mmax)。
步骤二:生成每张雷达图像的强度图
单视复数雷达图像的每一个像元都是复数,它包含了两个信息,即雷达目标的后向散射强度,简称强度信息,和表征雷达与目标距离的相位信息。为了计算像元强度的统计特征,要先计算出每一个像元的强度:
上式中,Ii,j表示雷达图像中位于第i行第j列像元的强度值,re表示像元的实部值,im表示像元的虚部值。
步骤三:生成序列差分干涉图
上式中,fm,fs和fms分别表示主图像,副图像和它们进行差分干涉后生成的差分干涉图,φsim是雷达干涉系统选定的参考面及参考面上地形地物在雷达覆盖区域上的模拟干涉相位,Im和Is分别表示主副图像的强度,φ表示差分干涉相位。当雷达图像较多,多于30张时,采用一张主图像的模式,否则采用多主图像模式,也就是所谓的短基线模式。
步骤四:确定滤波搜索窗口和统计邻域窗口的大小
所谓图像滤波,就是对图像中每一个像元要用一个新的值去替代它原来的值,这个新的值来源于对以该像元为中心的一个适度大小的窗口内所有其它像元的值的某种运算的结果。这个窗口叫滤波搜索窗口,通常取21行21列,或31行31列,如图1.中黑线框所示。为了计算差分相位的统计特性,需要定义一个像元的统计邻域窗口,如图1.中阴影所示,通常取5x5,或7x7个像元。
步骤五:计算差分干涉图中每个像元的强度累积分布函数
首先将在步骤二中生成的所有N张雷达强度图做成一个堆栈,如图2.所示。对每一个像元构造一个强度向量:
Us=(us1,us2,……usN)T (15)上式中,Us表示像元s的强度向量,us1,us2,……usN表示像元s在N张雷达强度图中的强度值的一个从小到大的排序。然后计算每一个像元的强度累积分布函数:
上式中,Fs(x)表示像元s在N张雷达强度图中的累积分布函数。理论上,雷达强度服从Gamma分布,但在大数据应用中,这种参数分布,即条件似然函数,计算复杂,耗时长,故本发明中,用无参数累积分布函数来近似替代其理论分布函数。
步骤六:确定差分干涉图中每个像元在其滤波搜索窗口内的强度分布意义下的同质像元这一步,是要确定,在一个搜索窗口内,哪些像元与该搜索窗口中心像元属于在强度概率分布意义下的同质像元,换句话说,是要确定,在一个搜索窗口内,哪些像元与该搜索窗口中心像元具有相同,至少是相似的强度累积分布函数。如果Fs(x)和Ft(x)分别表示搜索窗口中心像元s和搜索窗口内某一像元t的累积分布函数,它们之间的最大距离定义为这两个分布的相似性的量度,见下式,并如图3.所示:
由统计检验理论得知,两个累积分布的距离D服从如下式所示的KS分布:
KS分布Qks(λ)具有如下的性质:
Qks(0)=1及Qks(∞)=0 (19)可见,Qks(λ)是一个单调函数,如果两个累积分布函数的距离D对应的Qks(λ)大于某一个阈值,即D是足够的小,则认为被比较的累积分布是相似的,用到两个像元的强度的比较上,则认为这两个像元在强度上是同质的。图4.是搜索强度同质像元后其结果的一个示意图,在以像元i为中心的搜索窗口内有j,k,l,m,n,o,p等6个像元与它是强度意义下的同质像元,这些像元构成搜索窗口内以像元i相关联的强度同质区域,它是搜索窗口Ω的一个子集Ω’,这个子集可以是不规则排列的。每一个像元都具有它自己的强度同质子集,并且对堆栈中所有差分干涉图都有效。
步骤七:相位同质的度量
通过前面的步骤,差分干涉图中每一个像元都确定了一个在预先规定的搜索窗口内与该像元强度意义下的同质的子集,斑点噪声在这个子集内被去除了。步骤七,是要确定,在这个子集内,在差分相位概率分布的意义下,哪些像元与搜索窗口中心像元是同质的呢?在强度同质子集内,差分干涉相位服从高斯分布,它的噪声是一个均值为零的高斯噪声。居于此,本发明选用高斯函数来衡量两个像元差分相位的差的分布相似性的度量,用W来表示:
上式中,Ws,t可以看作是像元s和t之间差分相位分布的相似性,ω是每一个像元的统计邻域,如图4.中同质像元周边的阴影所示,通常是一个5x5或7x7的小窗口。φs,k和φt,k分别是像元s和t在其各自邻域ω内对应像元k的差分相位,σ是邻域内噪声的方差。可见:
0<Ws,t≤1 (21)
步骤八:差分干涉图的空间域滤波
上一步骤中,像元间差分相位分布相似的程度可以用(20)式中的Ws,t来度量,亦可以理解为它们在差分相位意义上同质性的度量,Ws,t还具有(21)式的性质,因此,本发明利用这个差分相位同质性的度量作为差分干涉图的空间域上的滤波参数—权(weight):
上式中,us是差分干涉图中像元s的复数值,Ω's是属于像元s的强度同质子集,ut是像元s的强度同质子集Ω's内的像元t的复数值,Weight(s,t)是像元t对于像元s所具有的权,它等于(20)式中的Ws,t。
作为总结,可以说,本发明提出的一个序列的差分干涉图的滤波同其它图像滤波方法的基础是一样的,都是利用了被滤波的像元的周边像元的加权平均。但本发明的特点在于,它没有将这个周边的形状固定,而是自适应地由周边那些在强度和差分相位概率分布意义下的同质像元所确定。因此,本发明提出的滤波方法,不仅可以同时消除差分干涉图中的乘性和加性噪声,同时还在很大程度上保留了雷达图像的原始分辨率。
附图说明
图1.搜索窗口和统计邻域窗口示意图,每一个小方格表示一个像元。搜索窗口是定义在以某一个像元,比如s,为中心的一个规则区域(M行,N列),统计邻域窗口是指在搜索窗口内某一像元,比如k,其周边一个较小的规则区域(m行,n列),如阴影所示。
图2.差分干涉图的堆栈及像元强度向量的结构示意图。
图3.累积分布函数间的距离定义示意图。
图4.在搜索窗口内的强度同质子集分布示意图。
图5.范家坪滑坡TerraSAR-X序列雷达图像强度均值图
图6.a范家坪滑坡20120107-20120129原始差分干涉图
图6.b范家坪滑坡20120107-20120129时空滤波后的差分干涉图
图7.本发明数据处理流程框图。
图中符号说明如下:
图1.中:s表示搜索窗口中心像元,k表示位于搜索窗口中的像元,它将被检验是否与像元s强度同质。
图2.中:N表示构成堆栈的差分干涉图数量,time表示按时间顺序堆栈。
图3.中:D表示两个累积分布函数之间的最大距离。
图4.中:i表示搜索窗口中的中心像元,j,k,l,m,n,o,及p表示搜索窗口内与i强度同质的像元。阴影表示强度同质像元的为下一步估计相位同质统计的邻域窗口。
具体实施方式
以长江三峡大坝地区范家坪滑坡地区INSAR形变监测中的差分干涉图的滤波为例,说明本发明在实际工程应用中的具体操作步骤。见图7,本发明是一种合成孔径雷达差分干涉图序列的时空同质滤波方法,该方法具体步骤如下:
步骤一:主副图像的自适应精密配准及副图像的重采样
为利用INSAR技术监测长江三峡大坝地区范家坪滑坡的非稳定特性,搜集了德国TerraSAR-X卫星X波段3米分辨率的21景单视复数(Single Look Complex)数据,时间跨距是从2012年1月至2013年2月。首先要选取一幅图像为主图像,作为参考图像,原则上它应该位于时域(时间基线集)和空间域(空间基线集)的中心。为了后续时序分析方便,这里选取第一张图像,即2012年1月7号的数据为主图像。以(1)式为模型基础,利用(2)式和(3)式表示的协方差方法计算其余各副图像中所有像元相对于主图像的相对偏移向量,用移位的Sinc函数对各副图像进行自适应二维重采样,以达到精密配准的目的,精度为十分之一个像元。
Offset(i,j)=(△i,△j)=((△i,△j)i,j)System+((△i,△j)i,j)Object (1)
式(1)中的第一项表示主图像中位于第i行第j列的像元与副图像中的同名点的位置的偏差量,该偏差量是一个二维向量,即式(1)中第二项。偏差向量是由系统产生的分量和主图像中该像元对应的地面目标本身形变产生的分量构成,分别为式(1)中第三和第四项所示。
在卫星轨道数和雷达成像时间等参数已知的情况下,可以先计算出系统产生的偏移分量。以该系统分量为基础,可以大致确定副图像中对应的像元位置。这样,对于主图像的中每一个像元(i,j),皆可以找到一个以该像元为中心的一个适当大小的窗口,比如512行和512列,同时亦可以确定副图像中对应的该窗口的位置。计算这两个窗口的协方差函数,如式(2)所示。M,N为窗口的行列数,p和s分别表示主副图像。使得协方差函数取得最大值时的m和n即为主图像中像元(i,j)在副图像中的偏移量,如式(3)所示。求出了所有副图像所有像元相对与主图像的偏移向量后再利用移位Sinc函数对所有幅图像进行重采样,使得所有幅图像皆与主图像精密配准。在该例中,感兴趣的地区仅是范家坪滑坡,故将所有21张雷达图像剪裁出以范家坪滑坡为中心的小区域图像,其包含2040(列)x1100(行)个像元。
步骤二:生成每张雷达图像的强度图
单视复数雷达图像的每一个像元都是复数,它包含了两个信息,即雷达目标的后向散射强度,简称强度信息,和表征雷达与目标距离的相位信息。为了计算像元强度的统计特征,要计算出每张图像中每一个像元的强度:
上式中,Ii,j表示雷达图像中,位于第i行第j列像元的强度值,re表示像元的实部值,im表示像元的虚部值。为了检查21张雷达图像配准后的效果,可以计算21张强度图像的均值图像,如图5.,图像清晰程度可以大致显示配准效果。
步骤三:以2012年1月7日图像为主图像,其余为副图像,生成序列差分干涉图
上式中,fm,fs和fms分别表示主图像,副图像和它们进行差分干涉后生成的差分干涉图,φsim是雷达干涉系统选定的参考面及参考面上地形地物在雷达覆盖区域上的模拟干涉相位。利用TerraSAR-X卫星的WGS84精密轨道数据和地球WGS84椭球体(其表面即为INSAR系统参考面)参数,根据雷达成像时间及成像参数,可以计算出参考面及参考面上地形的(DEM)引起的干涉相位φsim。根据(5)式即可计算出差分干涉图。Im和Is分别表示主副图像的强度,φ表示差分干涉相位。
步骤四:确定滤波搜索窗口和统计邻域窗口的大小
所谓图像滤波,就是对图像中每一个像元要用一个新的值去替代它原来的值,这个新的值来源于对以该像元为中心的一个适度大小的窗口内所有其它像元的值的某种运算的结果。这个窗口叫滤波搜索窗口,通常取21行21列,或31行31列。现实采用的TerraSAR-X数据像元大小为3米x3米,搜索窗口取21x21已够大,搜索范围已达3600m2.为了计算差分相位的统计特性,需要定义一个像元的统计邻域窗口,如图1.中阴影所示,这里取7x7个像元,以保证统计特性有较多样本单元。
步骤五:计算差分干涉图中每个像元的强度累积分布函数
像元位置在强度图和差分干涉图中是相同的,故差分干涉图中像元的累积分别函数是在序列强度图中进行的。首先将在步骤二中生成的所有N张雷达强度图做成一个堆栈,如图2.所示。对每一个像元构造一个强度向量:
Us=(us1,us2,……usN)T (6)上式中,Us表示像元s的强度向量,us1,us2,……usN表示像元s在N张雷达强度图中的强度值的一个从小到大的排序。然后计算每一个像元的强度累积分布函数:
上式中,Fs(x)表示像元s在N张雷达强度图中的累积分布函数,x是累积分布函数的自变量,如图3.所示。理论上,雷达强度服从Gamma分布,但在大数据应用中,这种参数分布,即条件似然函数,计算复杂,耗时长,故本发明中,用无参数累积分布函数以近似的方式替代其理论分布函数。这一步的输出,是差分干涉图中的每一个像元的累积分布函数Fs(x)。
步骤六:确定差分干涉图中每个像元在其滤波搜索窗口内的强度分布意义下的同质像元
这一步,是要确定,在一个搜索窗口内,哪些像元与该搜索窗口中心像元属于在强度概率分布意义下的同质像元,换句话说,是要确定,在一个搜索窗口内,哪些像元与该搜索窗口中心像元具有相同,至少是相似的强度累积分布函数。如果Fs(x)和Ft(x)分别表示搜索窗口中心像元s和搜索窗口内某一像元t的累积分布函数,它们之间的最大距离定义为这两个分布的相似性的量度,见下式,并如图3.所示:
在差分干涉图中,对每一个像元定义了一个以该像元为中心的大小相同的搜索窗口,其大小在步骤四中选定为21x21,共有441个像元,那么,除了该像元本身外,在其余440个像元中,哪些像元与该像元强度同质呢?这就要比较这两个像元的累积分布函数的距离。
由统计检验理论得知,两个累积分布的距离D服从如下式所示的KS分布:
KS分布Qks(λ)具有如下的性质:
Qks(0)=1及Qks(∞)=0 (10)可见,Qks(λ)是一个单调函数,如果两个累积分布函数的距离D对应的Qks(λ)大于某一个阈值,即D是足够的小,则认为被比较的累积分布是相似的,用到两个像元的强度的比较上,则认为这两个像元在强度上是同质的。图4.是搜索强度同质像元后其结果的一个示意图,在以像元i为中心的搜索窗口内有j,k,l,m,n,o,p等6个像元与它是强度意义下的同质像元,这些像元构成搜索窗口内以像元i相关联的强度同质区域,它是搜索窗口Ω的一个子集Ω’,这个子集可以是不规则排列的。每一个像元都具有它自己的强度同质子集,并且对堆栈中所有差分干涉图都有效。
根据以上描述,对于搜索窗口内某一像元t,搜索出它的累积分布函数Ft(x)与中心像元s的累积分布函数Fs(x)的距离D,如图3.所示,将其带入(9)式,计算出对应的Qks(λ),如果其值大于或等于阈值,就确定像元t是中性像元s在强度概率分布意义下的同质像元。阈值通常取0.95或更高。
这一步的输出,是确定了差分干涉图中每一个像元具有的在搜索窗口Ω内的强度同质像元的集合Ω'。
步骤七:相位同质的度量
通过前面的步骤,差分干涉图中每一个像元都确定了一个在预先规定的搜索窗口内与该像元强度意义下的同质的子集,斑点噪声在这个子集内被去除了。接下来,是要确定,在这个子集内,在差分相位概率分布的意义下,哪些像元与搜索窗口中心像元是同质的呢?在强度同质子集内,差分干涉相位服从高斯分布,它的噪声是一个均值为零的高斯噪声。居于此,选用高斯函数来衡量两个像元差分相位的差的分布相似性的度量,用W来表示:
在式(11)中,Ws,t可以看作是像元s和t之间差分相位分布的相似性,ω是每一个像元的统计邻域,如图4.中同质像元周边的阴影所示,通常是一个5x5或7x7的小窗口。φs,k和φt,k分别是像元s和t在其各自邻域ω内对应像元k的差分相位。需要强调的是,式(11)中求和运算是在邻域窗口内进行的。另外,与式(20)不同,因邻域内噪声的方差σ未知,在实际应用时,用滤波参数h来替代,它起到一个控制滤波程度的作用。由(11)式的指数函数性质可知:
0<Ws,t≤1 (12)
根据Ws,t在式(11)和(12)的性质,选用它作为对像元s进行相位滤波时像元t的权。这一步的输出,即是利用(11)式,计算在强度同质子集内,所有强度同质像元与搜索窗口中心像元间的相位分布的相似性。
步骤八:差分干涉图的空间域滤波
上一步骤中,像元间差分相位分布相似的程度可以用(11)式中的Ws,t来度量,亦可以理解为它们在差分相位意义上同质性的度量,Ws,t还具有(12)式的性质,因此,利用这个差分相位同质性的度量作为差分干涉图的空间域上的滤波参数—权(weight):
上式中,是差分干涉图中像元s的滤波后的复数值,Ω's是属于像元s的强度同质子集,ut是像元s的强度同质子集Ω's内的像元t的复数值,Weight(s,t)是像元t对于像元s所具有的权,它等于(11)式中的Ws,t。
图6.a是以范家坪滑坡TerraSAR-X的2012年1月7日图像和2012年1月29日图像干涉后的差分干涉图,差分干涉相位几乎被各种噪声淹没。图6.b是图6.a经过本发明的时空同质滤波以后的结果,可以见到十分清晰的高分辨率的雷达差分干涉相位。从黑到白的灰度变化显示了一个以2π为周期的差分相位变化。
可以看出,本发明提出的一个序列的差分干涉图的滤波同其它图像滤波方法的基础是一样的,都是利用了被滤波的像元的周边像元的加权平均。但本发明的特点在于,它没有将这个周边的形状固定,而是自适应地由周边那些在强度和差分相位概率分布意义下的同质像元所确定。因此,本发明提出的滤波方法,不仅可以同时消除差分干涉图中的乘性和加性噪声,同时还在很大程度上保留了雷达图像的原始分辨率。
Claims (1)
1.一种合成孔径雷达差分干涉图序列的时空同质滤波方法,其特征在于:该方法具体步骤如下:
步骤一:主副图像的自适应精密配准及副图像的重采样
在主副图像进行配准时,不仅要考虑雷达干涉系统的成像几何、干涉几何、雷达运行轨道引起的地面同一目标在主副图像上的位置的偏移,还要考虑地面目标因较大的高程或较大的形变位移引起的主副图像上的位置偏移;在此,副图像利用Sinc函数进行重采样,Sinc函数的位移量由下式决定:
Offset(i,j)=(△i,△j)=((△i,△j)i,j)System+((△i,△j)i,j)Object (10)上式中,offset(i,j)表示副图像中位于第i行第j列的像元相对于主图像上同一像元位置的偏移量,它是系统偏移和目标本身因高程或形变引起的偏移量的和;offset(i,j)由副图像s中位于(i,j)的像元的邻域与相应的主图像p中像元的邻域的协方差函数取最大值时的位移量决定,见下面的(11)和(12)式
上两式中,p和s分别表示主副图像,c表示在主副图像中以第i行第j列为中心的N行M列的两个窗口的协方差函数,n和m分别表示行和列偏移量,k和l表示窗口中参与相关计算像元的行列位置,(△i,△j)表示当c取得最大值时的行列号,即(nmax,mmax);
步骤二:生成每张雷达图像的强度图
单视复数雷达图像的每一个像元都是复数,它包含了两个信息,即雷达目标的后向散射强度,简称强度信息,和表征雷达与目标距离的相位信息;为了计算像元强度的统计特征,要先计算出每一个像元的强度:
步骤三:生成序列差分干涉图
上式中,fm,fs和fms分别表示主图像,副图像和它们进行差分干涉后生成的差分干涉图,φsim是雷达干涉系统选定的参考面及参考面上地形地物在雷达覆盖区域上的模拟干涉相位,Im和Is分别表示主副图像的强度,φ表示差分干涉相位;当雷达图像较多,多于30张时,采用一张主图像的模式,否则采用多主图像模式,也就是所谓的短基线模式;
步骤四:确定滤波搜索窗口和统计邻域窗口的大小
所谓图像滤波,就是对图像中每一个像元要用一个新的值去替代它原来的值,这个新的值来源于对以该像元为中心的一个适度大小的窗口内所有其它像元的值的某种运算的结果,这个窗口叫滤波搜索窗口,通常取21行21列,或31行31列;为了计算差分相位的统计特性,需要定义一个像元的统计邻域窗口,通常取5x5,或7x7个像元;
步骤五:计算差分干涉图中每个像元的强度累积分布函数
首先将在步骤二中生成的所有N张雷达强度图做成一个堆栈,对每一个像元构造一个强度向量:
Us=(us1,us2,……usN)T (15)上式中,Us表示像元s的强度向量,us1,us2,……usN表示像元s在N张雷达强度图中的强度值的一个从小到大的排序,然后计算每一个像元的强度累积分布函数:
上式中,Fs(x)表示像元s在N张雷达强度图中的累积分布函数,理论上,雷达强度服从Gamma分布,但在大数据应用中,这种参数分布,即条件似然函数,计算复杂,耗时长,故用无参数累积分布函数来近似替代其理论分布函数;
步骤六:确定差分干涉图中每个像元在其滤波搜索窗口内的强度分布意义下的同质像元
这一步,是要确定,在一个搜索窗口内,哪些像元与该搜索窗口中心像元属于在强度概率分布意义下的同质像元,换句话说,是要确定,在一个搜索窗口内,哪些像元与该搜索窗口中心像元具有相同,至少是相似的强度累积分布函数;如果Fs(x)和Ft(x)分别表示搜索窗口中心像元s和搜索窗口内某一像元t的累积分布函数,它们之间的最大距离定义为这两个分布的相似性的量度,见下式,
Qks(0)=1及Qks(∞)=0 (19)可见,Qks(λ)是一个单调函数,如果两个累积分布函数的距离D对应的Qks(λ)大于某一个阈值,即D是足够的小,则认为被比较的累积分布是相似的,用到两个像元的强度的比较上,则认为这两个像元在强度上是同质的;在以像元i为中心的搜索窗口内有j,k,l,m,n,o,p6个像元与它是强度意义下的同质像元,这些像元构成搜索窗口内以像元i相关联的强度同质区域,它是搜索窗口Ω的一个子集Ω’,这个子集是不规则排列的,每一个像元都具有它自己的强度同质子集,并且对堆栈中所有差分干涉图都有效;
步骤七:相位同质的度量
通过前面的步骤,差分干涉图中每一个像元都确定了一个在预先规定的搜索窗口内与该像元强度意义下的同质的子集,斑点噪声在这个子集内被去除了;步骤七是要确定,在这个子集内,在差分相位概率分布的意义下,哪些像元与搜索窗口中心像元是同质的,在强度同质子集内,差分干涉相位服从高斯分布,它的噪声是一个均值为零的高斯噪声;居于此,选用高斯函数来衡量两个像元差分相位的差的分布相似性的度量,用W来表示:
0<Ws,t≤1 (21)
步骤八:差分干涉图的空间域滤波
上一步骤中,像元间差分相位分布相似的程度用(20)式中的Ws,t来度量,亦理解为它们在差分相位意义上同质性的度量,Ws,t还具有(21)式的性质,因此,利用这个差分相位同质性的度量作为差分干涉图的空间域上的滤波参数—权(weight):
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310699869.5A CN103645476B (zh) | 2013-12-18 | 2013-12-18 | 一种合成孔径雷达差分干涉图序列的时空同质滤波方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201310699869.5A CN103645476B (zh) | 2013-12-18 | 2013-12-18 | 一种合成孔径雷达差分干涉图序列的时空同质滤波方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN103645476A true CN103645476A (zh) | 2014-03-19 |
CN103645476B CN103645476B (zh) | 2015-11-04 |
Family
ID=50250725
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201310699869.5A Expired - Fee Related CN103645476B (zh) | 2013-12-18 | 2013-12-18 | 一种合成孔径雷达差分干涉图序列的时空同质滤波方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN103645476B (zh) |
Cited By (12)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104122553A (zh) * | 2014-07-23 | 2014-10-29 | 中国国土资源航空物探遥感中心 | 一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法 |
CN105608674A (zh) * | 2015-12-16 | 2016-05-25 | 同济大学 | 一种基于图像配准、插值和去噪的图像增强方法 |
CN106780361A (zh) * | 2016-11-21 | 2017-05-31 | 上海航天测控通信研究所 | 一种用于城市区域合成孔径雷达图像的噪声抑制方法 |
CN109375222A (zh) * | 2018-12-17 | 2019-02-22 | 中国国土资源航空物探遥感中心 | 一种合成孔径雷达干涉测量电离层相位估计和补偿方法 |
CN110261839A (zh) * | 2019-07-04 | 2019-09-20 | 河海大学 | 一种基于双倍样本的增强谱分集方位向偏移量估计方法 |
CN110646792A (zh) * | 2019-11-04 | 2020-01-03 | 中国人民解放军空军工程大学 | 一种基于观察哨数字望远镜的雷达搜索窗口设置方法 |
WO2021009904A1 (ja) * | 2019-07-18 | 2021-01-21 | 日本電気株式会社 | 画像処理装置および画像処理方法 |
WO2021009905A1 (ja) * | 2019-07-18 | 2021-01-21 | 日本電気株式会社 | 画像処理装置および画像処理方法 |
CN112419198A (zh) * | 2020-11-27 | 2021-02-26 | 中国矿业大学 | 一种用于sar干涉图滤波的非局部均值定权方法 |
CN112711021A (zh) * | 2020-12-08 | 2021-04-27 | 中国自然资源航空物探遥感中心 | 一种多分辨率InSAR交互干涉时序分析方法 |
CN112904337A (zh) * | 2021-05-07 | 2021-06-04 | 北京东方至远科技股份有限公司 | 一种基于Offset Tracking技术的边坡形变时序监测方法 |
CN113272897A (zh) * | 2018-11-19 | 2021-08-17 | 珀金埃尔默健康科学有限公司 | 用于信号处理的降噪滤波器 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US4758838A (en) * | 1984-09-07 | 1988-07-19 | Hitachi, Ltd. | Method of reconstructing images from synthetic aperture radar's data |
US7245250B1 (en) * | 2005-08-16 | 2007-07-17 | Itt Manufacturing Enterprises, Inc. | Synthetic aperture radar image compression |
JP2011208974A (ja) * | 2010-03-29 | 2011-10-20 | Mitsubishi Electric Corp | レーダ画像処理装置 |
CN102565798A (zh) * | 2012-01-09 | 2012-07-11 | 中国民航大学 | 基于相关加权联合子空间投影的InSAR干涉相位估计方法 |
CN103076608A (zh) * | 2013-01-27 | 2013-05-01 | 西安电子科技大学 | 轮廓增强的聚束式合成孔径雷达成像方法 |
-
2013
- 2013-12-18 CN CN201310699869.5A patent/CN103645476B/zh not_active Expired - Fee Related
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US4758838A (en) * | 1984-09-07 | 1988-07-19 | Hitachi, Ltd. | Method of reconstructing images from synthetic aperture radar's data |
US7245250B1 (en) * | 2005-08-16 | 2007-07-17 | Itt Manufacturing Enterprises, Inc. | Synthetic aperture radar image compression |
JP2011208974A (ja) * | 2010-03-29 | 2011-10-20 | Mitsubishi Electric Corp | レーダ画像処理装置 |
CN102565798A (zh) * | 2012-01-09 | 2012-07-11 | 中国民航大学 | 基于相关加权联合子空间投影的InSAR干涉相位估计方法 |
CN103076608A (zh) * | 2013-01-27 | 2013-05-01 | 西安电子科技大学 | 轮廓增强的聚束式合成孔径雷达成像方法 |
Non-Patent Citations (2)
Title |
---|
夏耶: "Bam earthquake: Surface deformation measurement using radar interferometry", 《ACTA SEISMOLOGICA SINICA》 * |
毛建旭等: "合成孔径雷达干涉成像技术及其应用", 《系统工程与电子技术》 * |
Cited By (22)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN104122553A (zh) * | 2014-07-23 | 2014-10-29 | 中国国土资源航空物探遥感中心 | 一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法 |
CN105608674A (zh) * | 2015-12-16 | 2016-05-25 | 同济大学 | 一种基于图像配准、插值和去噪的图像增强方法 |
CN106780361A (zh) * | 2016-11-21 | 2017-05-31 | 上海航天测控通信研究所 | 一种用于城市区域合成孔径雷达图像的噪声抑制方法 |
CN113272897A (zh) * | 2018-11-19 | 2021-08-17 | 珀金埃尔默健康科学有限公司 | 用于信号处理的降噪滤波器 |
CN109375222A (zh) * | 2018-12-17 | 2019-02-22 | 中国国土资源航空物探遥感中心 | 一种合成孔径雷达干涉测量电离层相位估计和补偿方法 |
CN110261839B (zh) * | 2019-07-04 | 2023-02-28 | 河海大学 | 一种基于双倍样本的增强谱分集方位向偏移量估计方法 |
CN110261839A (zh) * | 2019-07-04 | 2019-09-20 | 河海大学 | 一种基于双倍样本的增强谱分集方位向偏移量估计方法 |
WO2021009904A1 (ja) * | 2019-07-18 | 2021-01-21 | 日本電気株式会社 | 画像処理装置および画像処理方法 |
WO2021009905A1 (ja) * | 2019-07-18 | 2021-01-21 | 日本電気株式会社 | 画像処理装置および画像処理方法 |
JPWO2021009905A1 (zh) * | 2019-07-18 | 2021-01-21 | ||
JPWO2021009904A1 (zh) * | 2019-07-18 | 2021-01-21 | ||
JP7188594B2 (ja) | 2019-07-18 | 2022-12-13 | 日本電気株式会社 | 画像処理装置および画像処理方法 |
JP7188595B2 (ja) | 2019-07-18 | 2022-12-13 | 日本電気株式会社 | 画像処理装置および画像処理方法 |
US11846702B2 (en) | 2019-07-18 | 2023-12-19 | Nec Corporation | Image processing device and image processing method |
US20220268922A1 (en) * | 2019-07-18 | 2022-08-25 | Nec Corporation | Image processing device and image processing method |
CN110646792B (zh) * | 2019-11-04 | 2022-04-12 | 中国人民解放军空军工程大学 | 一种基于观察哨数字望远镜的雷达搜索窗口设置方法 |
CN110646792A (zh) * | 2019-11-04 | 2020-01-03 | 中国人民解放军空军工程大学 | 一种基于观察哨数字望远镜的雷达搜索窗口设置方法 |
CN112419198A (zh) * | 2020-11-27 | 2021-02-26 | 中国矿业大学 | 一种用于sar干涉图滤波的非局部均值定权方法 |
CN112419198B (zh) * | 2020-11-27 | 2024-02-02 | 中国矿业大学 | 一种用于sar干涉图滤波的非局部均值定权方法 |
CN112711021A (zh) * | 2020-12-08 | 2021-04-27 | 中国自然资源航空物探遥感中心 | 一种多分辨率InSAR交互干涉时序分析方法 |
CN112904337B (zh) * | 2021-05-07 | 2021-11-05 | 北京东方至远科技股份有限公司 | 一种基于Offset Tracking技术的边坡形变时序监测方法 |
CN112904337A (zh) * | 2021-05-07 | 2021-06-04 | 北京东方至远科技股份有限公司 | 一种基于Offset Tracking技术的边坡形变时序监测方法 |
Also Published As
Publication number | Publication date |
---|---|
CN103645476B (zh) | 2015-11-04 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN103645476B (zh) | 一种合成孔径雷达差分干涉图序列的时空同质滤波方法 | |
Almonacid-Caballer et al. | Evaluation of annual mean shoreline position deduced from Landsat imagery as a mid-term coastal evolution indicator | |
Debella-Gilo et al. | Sub-pixel precision image matching for measuring surface displacements on mass movements using normalized cross-correlation | |
Emmanuel et al. | Temporal and spatial variability of rainfall at the urban hydrological scale | |
CN101727662B (zh) | Sar图像非局部均值去斑方法 | |
EP3639054A1 (en) | Land deformation measurement | |
Shu et al. | Shoreline extraction from RADARSAT-2 intensity imagery using a narrow band level set segmentation approach | |
CN101650439B (zh) | 基于差异边缘和联合概率一致性的遥感图像变化检测方法 | |
CN103822598A (zh) | 地基sar在时间去相关严重区域的形变监测方法 | |
Debella-Gilo et al. | Locally adaptive template sizes for matching repeat images of Earth surface mass movements | |
CN101980293A (zh) | 一种基于刃边图像的高光谱遥感系统mtf检测方法 | |
CN102472815A (zh) | 对从在相同区域上采集的sar图像获得的干涉图进行滤波的方法 | |
Schumann et al. | Downscaling coarse grid hydrodynamic model simulations over large domains | |
CN109613531B (zh) | 一种微变感知预警雷达的多阈值优化形变反演方法及系统 | |
Gabrieli et al. | A low-cost landslide displacement activity assessment from time-lapse photogrammetry and rainfall data: Application to the Tessina landslide site | |
CN106197383A (zh) | 一种海冰体积的遥感估算方法 | |
CN111982822B (zh) | 一种长时间序列高精度植被指数改进算法 | |
Jiang et al. | Application of multitemporal InSAR covariance and information fusion to robust road extraction | |
Jiang et al. | Delineation of built-up land change from SAR stack by analysing the coefficient of variation | |
CN114091274A (zh) | 一种滑坡易发性评价方法及系统 | |
CN112051572A (zh) | 一种融合多源sar数据三维地表形变监测方法 | |
Girón et al. | Nonparametric edge detection in speckled imagery | |
Ju et al. | An improved algorithm for computing local fractal dimension using the triangular prism method | |
Gao et al. | Monitoring terrain elevation of intertidal wetlands by utilising the spatial-temporal fusion of multi-source satellite data: A case study in the Yangtze (Changjiang) Estuary | |
CN101893540B (zh) | 一种基于八邻域算法估算悬沙反演模型尺度误差的方法 |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
C10 | Entry into substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
C14 | Grant of patent or utility model | ||
GR01 | Patent grant | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20151104 Termination date: 20171218 |
|
CF01 | Termination of patent right due to non-payment of annual fee |