CN105425216B - 基于图像分割的重复航过极化InSAR图像配准方法 - Google Patents

基于图像分割的重复航过极化InSAR图像配准方法 Download PDF

Info

Publication number
CN105425216B
CN105425216B CN201510825264.5A CN201510825264A CN105425216B CN 105425216 B CN105425216 B CN 105425216B CN 201510825264 A CN201510825264 A CN 201510825264A CN 105425216 B CN105425216 B CN 105425216B
Authority
CN
China
Prior art keywords
image
subgraph
piecemeal
piece
offset
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
Application number
CN201510825264.5A
Other languages
English (en)
Other versions
CN105425216A (zh
Inventor
李真芳
解金卫
张芳
王跃锟
张海赢
郭媛
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Xidian University
Original Assignee
Xidian University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Xidian University filed Critical Xidian University
Priority to CN201510825264.5A priority Critical patent/CN105425216B/zh
Publication of CN105425216A publication Critical patent/CN105425216A/zh
Application granted granted Critical
Publication of CN105425216B publication Critical patent/CN105425216B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021SAR image post-processing techniques
    • G01S13/9023SAR image post-processing techniques combined with interferometric techniques
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Details of systems according to groups G01S13/00, G01S15/00, G01S17/00
    • G01S7/02Details of systems according to groups G01S13/00, G01S15/00, G01S17/00 of systems according to group G01S13/00
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01SRADIO 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/00Systems 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/88Radar or analogous systems specially adapted for specific applications
    • G01S13/89Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/904SAR modes
    • G01S13/9076Polarimetric features in SAR

Abstract

本发明属于极化干涉合成孔径雷达探测技术领域,公开了一种基于图像分割的重复航过极化InSAR图像配准方法,包括:获取极化干涉合成孔径雷达的参考图像和预配准图像;计算参考图像和预配准图像的幅度互相关矩阵,并获取参考图像和预变换后的图像的公共位置部分分别作为主图像和辅图像;对主图像和辅图像进行分级配准,求得辅图像中各级子图像相对于主图像中各级子图像的偏移量;根据辅图像中各级子图像相对于主图像中各级子图像的偏移量,计算辅图像相对于主图像的整体偏移量;根据辅图像相对于主图像的整体偏移量,计算辅图像中各个像素点的值,得到配准后的辅图像。

Description

基于图像分割的重复航过极化InSAR图像配准方法
技术领域
本发明涉及极化干涉合成孔径雷达探测技术领域,尤其涉及一种基于图像分割的重复航过极化InSAR图像配准方法,可用于机载雷达对同一场景中重复航过速度或基线差异较大的情况下,图像配准时偏移量存在严重方位向或距离向空变的极化InSAR图像的精确配准。
背景技术
极化干涉合成孔径雷达技术是现代雷达遥感测量中一项前沿技术。它能够充分利用雷达回波所携带的幅度信息、相位信息和极化散射信息,其不仅具有干涉合成孔径雷达对空间分布和高程的敏感性而且具有极化合成孔径雷达对散射体形状和方向的敏感性的特点,使得其在地表沉降监测、地表参数反演、地表地物变化检测中得到了广泛的应用。
无论是地表沉降检测、参数反演还是地物变化检测,图像配准是极化干涉合成孔径雷达图像预处理中的关键步骤,配准精度直接影响干涉相位图的质量。重复航过中同场景的极化InSAR复图像由于各小场景成像时采用的载机速度、多普勒中心和下视角不同等原因造成的不同区域之间旋转、平移和局部形变等变化程度不同,从而使得图像整体配准精度下降,而不精确的配准及图像失配会导致干涉相干系数下降、干涉相位条纹不清晰甚至消失,从而给后续的测量中引入较大误差。现有的极化InSAR图像配准方法主要处理步骤都是以一幅图像为主图像,在主图像中随机选取若干区域,这些区域作为局部匹配窗,然后在辅图像中选定一个较大的搜索窗用匹配窗按一定质量引导规则进行搜索,得到若干偏移量后用多项式拟合出辅图像中所有像素点位置,最后利用插值得到配准后的辅图像。
齐海宁在文章“一种结合最优相干运算的极化干涉SAR相干配准方法”(遥感技术与应用,2004,19(6):512-516)中先利用全局最优相干过程将两组全极化数据窗口转化到最优极化状态下,然后对所得窗口数据以相干系数作匹配质量指标做配准。该方法能够从复杂的散射机制中抽取出占主导地位的散射机制,从而使散射中心得以确定,这样可以有效消除散射中心高度差去相干的影响,然而该方法的不足之处是:只有在待配准图像没有太大形变的情况下,该配准方法比较精确。当两航过前后速度差异较大、基线不平行或者基线较长时,成像所得前后两航过SAR图像形变差异太大,该方法配准精度不理想,难以达到要求。
熊涛在文章“极化干涉合成孔径雷达应用的关键技术研究”(清华大学工学博士学位论文,2009)将每一航过的极化散射矩阵写成散射矢量形式,利用散射矢量构造了一种相似性参数,用该参数作为配准的质量引导因子对极化InSAR图像对进行了配准研究,该方法有效利用了极化InSAR中的极化信息,相对于传统的配准方法配准精度有了较大的提高,然而上述方法的不足之处是:在进行配准时,假设图像的极化通道间是精确配准的,而且空间通道上每一个搜索窗内有唯一的相似性参数峰值;然而通常情况下,极化通道间因为系统的极化通道不完全正交等因素,定标后的极化通道间仍然存在误差;空间通道上,当地物类型单一时,同一个搜索窗内相似性参数往往可能出现多个相同的峰值或者搜索窗内的峰值并不能反映真实的偏移量。
发明内容
本发明针对上述现有技术存在的在两次航过中载机速度或基线差异较大情景下,因为图像旋转伸缩现象严重使得极化InSAR图像配准不精确的缺陷,提出了一种基于图像分割的重复航过极化InSAR图像配准方法,通过分级分块进行配准处理,保证了在方位向和距离向均存在较大偏移量空变性的情况下,可以有效消除图像的中的伸缩旋转问题,使得两幅图像的共同部分精确配准从而获得质量较好的配准结果。
为达到上述目的,本发明的实施例采用如下技术方案予以实现。
一种基于图像分割的重复航过极化InSAR图像配准方法,所述方法包括如下步骤:
步骤1,获取极化干涉合成孔径雷达的参考图像和预配准图像,所述参考图像和预配准图像为所述极化干涉合成孔径雷达两次重复航过中对同一区域的成像结果;
步骤2,计算所述参考图像和所述预配准图像的幅度互相关矩阵,根据所述幅度互相关矩阵对所述预配准图像进行预变换,得到预变换后的图像,并获取所述参考图像中预设像素区域的部分作为主图像,获取所述预变换后的图像中与所述参考图像初步配准的部分作为辅图像;
步骤3,对所述主图像和所述辅图像进行分级分块并配准,求得所述辅图像中各级子图像相对于所述主图像中各级子图像的偏移量;
步骤4,根据所述辅图像中各级子图像相对于所述主图像中各级子图像的偏移量,计算所述辅图像相对于所述主图像的整体偏移量;
步骤5,根据所述辅图像相对于所述主图像的整体偏移量,计算所述辅图像中各个像素点的值,得到配准后的辅图像。
本发明与现有的技术相比具有以下优点:第一,本发明对重复航过极化干涉合成孔径雷达(InSAR)数据的图像配准采用基于图像分割的多级快速配准技术,在成像所得前后两航过极化InSAR图像形变差异较大的情况下有效克服了主辅图像间偏移量空变严重无法精确配准的难题,得到了清晰的干涉条纹图。第二,本发明对不同航过间极化InSAR数据进行配准时,采用了相关系数和复相干系数结合做配准的质量引导指标,既利用了图像的幅度信息也利用了图像的相位信息,当地物类型单一时,在图像前几级分块中因为图像包含的纹理信息较多使用幅度相关法可以得到准确的偏移量,在最后一级分块中因为图像包含像素较少纹理信息也减少,使用复相干系数法加入图像间的相位信息,也能准确搜索到反映图像偏移位置的质量引导指标峰值的位置。因此配准具有很好的鲁棒性。
附图说明
为了更清楚地说明本发明实施例或现有技术中的技术方案,下面将对实施例或现有技术描述中所需要使用的附图作简单地介绍,显而易见地,下面描述中的附图仅仅是本发明的一些实施例,对于本领域普通技术人员来讲,在不付出创造性劳动的前提下,还可以根据这些附图获得其他的附图。
图1为本发明实施例提供的基于图像分割的重复航过极化InSAR图像配准方法的流程示意图;
图2为两航过HH极化通道实测数据的参考图像和预配准图像;
图3为实测数据传统粗配准方式主辅图像生成的相干系数图;
图4为实测数据传统粗配准后进行亚像素级精配准方式主辅图像生成的相干系数图;
图5为实测数据本发明多级配准方法主辅图像生成的相干系数图;
图6为实测数据本传统粗配准、亚像素精配准和本发明分级配准相干系数统计图对比;
图7为实测数据传统粗配准加精配准后主辅图像生成的干涉相位图
图8为实测数据本发明配准完成后主辅图像生成的干涉相位图。
具体实施方式
下面将结合本发明实施例中的附图,对本发明实施例中的技术方案进行清楚、完整地描述,显然,所描述的实施例仅仅是本发明一部分实施例,而不是全部的实施例。基于本发明中的实施例,本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例,都属于本发明保护的范围。
参照附图1,本发明的具体实施步骤如下:
步骤1,获取极化干涉合成孔径雷达的参考图像和预配准图像。
所述参考图像和预配准图像为所述极化干涉合成孔径雷达两次重复航过中对同一区域的成像结果。
将双航过极化干涉合成孔径雷达分别得到的HH极化通道图像输入到系统中,第一航过得到的HH极化通道图像作为参考图像数据,第二航过得到的HH极化通道图像作为预配准图像数据。输入的参考图像数据和预配准图像数据是对同一地区的成像结果,以此来保证参考图像与预配准图像之间通过本发明配准后具有一定相干性。
步骤2,计算所述参考图像和所述预配准图像的幅度互相关矩阵,根据所述幅度互相关矩阵对所述预配准图像进行预变换,得到预变换后的图像,并获取所述参考图像中预设像素区域的部分作为主图像,获取所述预变换后的图像中与所述参考图像初步配准的部分作为辅图像。
利用幅度互相关系数法对两航过同极化通道(HH极化通道)的参考图像和预配准图像进行预配准处理,然后取出预配准后的两幅图像的公共位置部分,在图像预配准的处理过程中配准精度只要求达到像素级。
幅度互相关系数法具体步骤如下:
第一步,利用下式计算参考图像和预配准图像的幅度互相关矩阵:
L=fftshift(ifft2(fft2(abs(M))×conj(fft2(abs(S)))))
其中,L表示参考图像和预配准图像的幅度互相关矩阵,fft2(·)表示对图像数据做二维傅里叶变换,ifft2(·)表示对图像数据做二维逆傅里叶变换,M表示参考图像数据,S表示预配准图像数据,conj(·)表示取共轭操作,abs(·)表示取矩阵元素的幅度值,fftshift(·)将数据先上半部分和下半部分互换然后左半部分和右半部分互换。
第二步,找出幅度互相关矩阵L所有元素中的最大值,得到互相关矩阵L中最大值元素对应的坐标位置。
第三步,用互相关矩阵L的中心坐标位置的值减去幅度互相关矩阵L中元素最大值的坐标位置,得到预配准后的偏移量。
第四步,根据所述预配准所需的偏移量对预配准图像的数据整体做平移,得到平移后的预配准图像,然后截取参考图像中预设像素区域的部分作为主图像,获取预配准图像中与所述参考图像初步配准的部分作为辅图像。
步骤3,对所述主图像和所述辅图像进行分级分块并配准,求得所述辅图像中各级子图像相对于所述主图像中各级子图像的偏移量。
步骤3具体包括如下子步骤:
(3a)对主图像和辅图像的上下左右各预留预设数量的像素;
(3b)对主图像沿距离向和方位向分别进行第一级分块,并对辅图像分别沿距离向和方位向分别进行第一级分块,求解第一级分块后辅图像中每块子图像对应主图像中每块子图像的偏移量,并记录主图像第一级分块中每块子图像的中心点位置坐标;
(3c)对主图像进行第一级分块后的每块子图像再进行下一级分块,按照第一级分块中辅图像的每块子图像对应主图像的每块子图像的偏移量,对辅图像的每块子图像加上偏移量并进行下一级分块,求解下一级分块后辅图像中每块子图像对应主图像中每块子图像的偏移量,并记录主图像下一级分块中每块子图像的中心点位置坐标;
(3d)依次类推,对上一级分块后的子图像再进行下一级分块,直到下一级分块后得到的子图像的像素个数小于设定的像素阈值时停止分级;然后,预设偏移量阈值,得到倒数第二级分块后计算得到的辅图像中子图像的偏移量大于预设偏移量阈值的子图像,并令该子图像的偏移量为零;
(3e)对主图像和辅图像最后一级分块得到的每块子图像去坡度;
所述对主图像和辅图像最后一级分块得到的每块子图像去坡度,具体包括:
(1)利用下式求取主图像和辅图像最后一级分块得到的每块子图像的干涉图,
subms=subm.*conj(subs)
其中,subm表示主图像最后一级分块得到的任一子图像,subs表示辅图像最后一级分块得到与所述主图像的任一子图像对应的子图像,conj(·)表示取共轭操作,subms表示主图像和辅图像最后一级分块得到的子图像的干涉图;
(2)对干涉图沿每一个方位向做N1点傅里叶变换,N1大于等于方位向点数,然后分别沿每一个距离向复数相加,得到M1个复数相加后的值,M1为距离向点数,得到M1个复数相加后的值中的最大值A;
对干涉图沿每一个距离向做N2点傅里叶变换,N2大于等于距离向点数,然后分别沿每一个方位向复数相加,得到M2个复数相加后的值,M2为方位向点数,得到M2个复数相加后的值中的最大值B,当2π*A/N或2π*B/N的值大于π时要将其减去2π,保证其值在[-π,π]范围内;
(3)利用下式得到主图像和辅图像最后一级分块得到的每块子图像的二维坡度,
slop=exp(-j*2π*A/N*P1).'*exp(-j*2π*B/N*P2)
其中,j表示虚数单位(·).'表示转置,exp(·)表示取相位值,slop表示子块的二维坡度,P1表示子图像在方位向上的点数,P2表示子图像在距离向上的点数;
(4)将主图像和辅图像中每块子图像的二维坡度与每块子图像的干涉图做内积去除子块的坡度。
(3f)对主图像和辅图像最后一级分块得到的每块子图像进行配准,求解最后一级分块后辅图像中每块子图像对应主图像中每块子图像的偏移量。
对于各区域偏移量存在严重空变的SAR图像,经预配准无法使得所有像素达到一个像素的配准精度,需要进一步进行多级配准,此处假设经预配准后图像大小为4120*4120。
多级配准的具体步骤为:
(a)将取出的主辅图像(即主图像和辅图像)的公共位置部分上下左右各空出12个像素量,取出方位向4096点距离向4096点的中间部分作为待分级的主辅图像块;
(b)将待分级主辅图像进行第一级分块,第一级分块可将主辅图像沿距离向和方位向平均分成2*2块,每块大小为2048*2048;对应主辅图像块按预配准方式分别做整体幅度相关求偏移量,求出第一级每块对应距离向及方位向偏移量数值并分别记录在X1、Y1向量中,主图像第一级每一块的中心点位置坐标分别记录在x1、y1向量中;
(c)将主图像第一级分得的子块分别按照类似于第一级的分块方法进行第二级分块,即主图像第一级每一子块再分成2*2个小块,每块大小为1024*1024,辅图像要取第一级得到的子块整体加上该子块在第一级得到的偏移量所对应区域进行第二级分块,每块大小为1024*1024,同样利用幅度相关系数法求出第二级每块对应距离向及方位向相对于第一级的偏移量数值并分别记录在X2、Y2向量中,主图像第二级每一块的中心点位置坐标分别记录在x2、y2向量中;
(d)依次进行第三级、第四级等分块配准,第三级主图像分块在第二级主图像子块上进行,第三级辅图像分块在第二级辅图像子块整体加第二级偏移量的位置上进行,第四级主图像分块在第三级主图像子块上进行,第四级辅图像分块在第三级辅图像子块整体加第三级偏移量的位置上进行,每一级分块大小为2*2,每一级对应距离向及方位向相对于上一级的偏移量数值并分别记录在Xn、Yn,(n为分块的级数)向量中,主图像每一级所分的块的中心点位置坐标分别记录在xn、yn,(n为分块的级数)向量中。
依此类推,直到所分块大小小于或等于设定阈值像素数T时(这里T设为16*16,即块大小小于16*16时认为块内信息太少不适合继续分级)停止分块,最后一次分块前要根据预先设置的偏移量门限剔除相关性低于门限P(这里设置门限为5,即最后一级主辅图像子块的偏移量大于5个像素时求得的偏移量不准确)的偏移量,即将该级偏移量置零,为了避免分级的级数太多造成数据量庞大造成配准效率下降,在第四级分块时由于前三级分块配准已经基本将图像配准的一个像素内可以直接将第四级分为32*32块,每块大小16*16;
(e)将主辅图像最后一级也就是第四级得到的子块进行去坡度,这里指的主辅图像干涉相位坡度;
去坡度的目的是避免干涉相位条纹太细而造成的利用图像相位进行配准时误差太大,去坡度的具体步骤为:
(1)利用下式求取第四级子块主辅图像的干涉图;
subms=subm.*conj(subs)
其中,subm表示第四级子块主图像,subs表示第四级子块辅图像,conj(·)表示取共轭操作,subms表示第四级子块的干涉图;
(2)对干涉图的每一个方位向做N(N大于等于方位向点数,这里设N为512)点FFT,然后沿距离向复数相加,最后得到复数相加后的最大值A,对干涉图的每一个距离向做N(N大于等于方位向点数,这里设N为512)点FFT,然后沿方位向复数相加,最后得到复数相加后的最大值B,当2π*A/N或2π*B/N的值大于π时要将其减去2π,保证其值在[-π,π]范围内;
(3)利用下式得到子块的二维坡度;
slop=exp(-j*2π*A/N*16).'*exp(-j*2π*B/N*16)
其中,j表示虚数单位(·).'表示转置,exp(·)表示取相位值,slop表示子块的二维坡度;
(4)将得到的二维坡度与子块干涉图进行做内积就去除了子块的坡度。
(f)将去坡度后的第四级子块主图像利用双线性内插法进行十倍插值,辅图像相应子块以块中心点位置取15*15大小的窗口,插值后的主图像从方位向和距离向第一点开始每隔10个像素进行采样也取出15*15大小的窗,计算主辅图向两个窗采用质量引导因子为带相位信息的复相干系数,其中复相干系数按如下公式计算:
其中,s1为主图像,s2表示幅图像,E(·)表示取空间平均,辅图像在主图像中按行列搜索求取最大复相干系数的位置,从而得到正负一个像素内偏移量并将其作为最后一级偏移量,求出第四级辅图像每块对应距离向及方位向相对于第三级的偏移量数值并分别记录在X4、Y4向量中,主图像第四级每一块的中心点位置坐标分别记录在x4、y4向量中。
步骤4,根据所述辅图像中各级子图像相对于所述主图像中各级子图像的偏移量,计算所述辅图像相对于所述主图像的整体偏移量。
步骤4具体包括如下子步骤:
(4a)将辅图像中各级分块得到子图像的偏移量与其所属的上一级子图像的偏移量进行累加,从而得到最后一级分块得到的每块子图像的最终偏移量;
(4b)将最后一级分块得到的每块子图像的最终偏移量按照主图像最后一级子图像的中心位置坐标重新排列成矩阵形式,形成偏移量矩阵,对偏移量矩阵做中值滤波;
(4c)对中值滤波后的偏移量矩阵再进行双线性插值,得到辅图像每个像素点对应的偏移量。
具体的,以步骤3中的示例为基础,步骤4具体包括:
(a)将每一级分块得到的距离向和方位向偏移量按下式累加:
X=(X1+X2+X3+X4)
Y=(Y1+Y2+Y3+Y4)
其中,X为距离向总的偏移量,由于每一级分块存储的偏移量向量长度不同,偏移量累加时需要将将前面三级的偏移量向量长度扩展,使得其长度和第四级偏移量向量长度相同,以此保证相加后向量内元素是每一级对应子块的偏移量之和。
(b)将累加后的距离向和方位向偏移量按照主图像最后一级子块中心的位置坐标重新排列成矩阵形式XX和YY,对偏移量矩阵做二维中值滤波,滤波窗口大小为5*5,使偏移量更加平滑;
(c)对滤波后的偏移量进行双线性插值,得到辅图像每一像素点对应的距离向和方位向偏移量,两个偏移量矩阵大小均为4096*4096。
步骤5,根据所述辅图像相对于所述主图像的整体偏移量,计算所述辅图像中各个像素点的值,得到配准后的辅图像。
按照步骤4中计算所得偏移量按下式求辅图像最后一级子块对应整个主图像偏移后位置:
Xpos=x4-XX;
Ypos=y4-YY;
根据求得的辅图像像素偏移后位置对整个辅图像进行双线性插值重采样,得到偏移后像素点对应值,从而得到配准后的辅图像。
下面结合实测数据实验对本发明的效果做进一步的说明。
处理实测数据时,使用的数据是中电十四所重复航过同一场景的机载X波段HH极化通道的极化InSAR数据。
图2(a)表示配准前的第一航过的HH极化通道极化InSAR数据图像,为主图像,图2(b)表示配准前的第二航过的HH极化通道极化InSAR数据图像,为辅图像,两图的横坐标轴表示距离向,单位为像素,纵坐标轴表示方位向,单位为像素。对比图2(a)和图2(b)能够明显看到图像中的一条主干道呈现旋转形变现象,其他建筑物呈现较大位置平移现象;图3表示传统对主辅图像做整体幅度互相关求去偏移量配准后主辅图像的相干系数图,横坐标轴表示距离向,单位为像素,纵坐标轴表示方位向,单位为像素,图像中的像素值在图像右侧色度条中的位置反应了相干系数的大小,也反应出了粗配准后主辅图像的相关性;图4为实测数据传统粗配准后进行亚像素级精配准方式主辅图像生成的相干系数图,横坐标轴表示距离向,单位为像素,纵坐标轴表示方位向,单位为像素;图5为实测数据本发明多级配准方法主辅图像生成的相干系数图,横坐标轴表示距离向,单位为像素,纵坐标轴表示方位向,单位为像素。图3、图4和图5对比可以看出,传统对主辅图像做整体幅度互相关求取偏移量配准后主辅图像的相干系数图相干性很低只有中间一部分相干系数值较高,高于0.5,大部分像素值低于0.5。传统粗配准后进行亚像素级精配准方式主辅图像生成的相干系数图比对主辅图像做整体幅度互相关求取偏移量配准后主辅图像的相干系数值略有提高,但是整体还是偏低,大部分低于0.5。而本发明多级配准方法主辅图像生成的相干系数图整体较前两种方法有所提高,0.5以上的像素数明显增多;图6为实测数据本传统粗配准、粗配准后的亚像素精配准和本发明分级配准相干系数统计图对比,横坐标为相干系数的值,纵坐标表示该区域各个值的相干系数像素数占该区域总像素数的百分比,单位为1。其中,十字标注的曲线代表传统粗配准方法得到主辅图像的相干系数图中各个值的像素数占该区域总像素数的百分比,圆圈标注的曲线代表传统粗配准后的亚像素精配准方法得到主辅图像的相干系数图中各个值的像素数占该区域总像素数的百分比,五角星标注的曲线代本发明分级配准方法得到主辅图像的相干系数图中各个值的像素数占该区域总像素的百分比,由相干系数统计图可以看出传统粗配准方法和粗配准后的亚像素精配准方法得到的主辅图像的相干系数值大部分分布在0.15左右,而本发明得到的主辅图像的相干系数值大部分在0.3左右,本发明得到的主辅图像的相干系数较高的像素数占总像素百分比明显高于前两种方法占总像素的百分比;图7为实测数据传统粗配准加亚像素级精配准后主辅图像生成的干涉相位图,横坐标轴表示距离向,单位为像素,纵坐标轴表示方位向,单位为像素;图8为实测数据本发明配准完成后主辅图像生成的干涉相位图,横坐标轴表示距离向,单位为像素,纵坐标轴表示方位向,单位为像素。对比图7和图8,可以看出传统粗配准加亚像素级精配准得到的主辅图像求取的干涉相位图基本看不到干涉条纹,而本发明分级配准方法配准后的主辅图像求取的干涉相位图有干涉条纹,且干涉条纹比较清晰。
由上述分析可知,本发明处理偏移量空变量较大的极化InSAR实测数据时,仍然可以得到比现有处理方法质量更好的干涉相位图,进一步验证了本发明对配准误差有很好的鲁棒性,并且说明了本发明的实用性和有效性。
以上所述,仅为本发明的具体实施方式,但本发明的保护范围并不局限于此,任何熟悉本技术领域的技术人员在本发明揭露的技术范围内,可轻易想到变化或替换,都应涵盖在本发明的保护范围之内。因此,本发明的保护范围应以所述权利要求的保护范围为准。

Claims (5)

1.一种基于图像分割的重复航过极化InSAR图像配准方法,其特征在于,所述方法包括如下步骤:
步骤1,获取极化干涉合成孔径雷达的参考图像和预配准图像,所述参考图像和预配准图像为所述极化干涉合成孔径雷达两次重复航过中HH极化通道对同一区域的成像结果;
其中,所述参考图像为第一航过得到的HH极化通道图像,所述预配准图像为第二航过得到的HH极化通道图像;
步骤2,计算所述参考图像和所述预配准图像的幅度互相关矩阵,根据所述幅度互相关矩阵对所述预配准图像进行预变换,得到预变换后的图像,并获取所述参考图像中预设像素区域的部分作为主图像,获取所述预变换后的图像中与所述参考图像初步配准的部分作为辅图像;
步骤3,对所述主图像和所述辅图像进行分级分块并配准,求得所述辅图像中各级子图像相对于所述主图像中各级子图像的偏移量;
步骤4,根据所述辅图像中各级子图像相对于所述主图像中各级子图像的偏移量,计算所述辅图像相对于所述主图像的整体偏移量;
步骤5,根据所述辅图像相对于所述主图像的整体偏移量,计算所述辅图像中各个像素点的值,得到配准后的辅图像。
2.根据权利要求1所述的一种基于图像分割的重复航过极化InSAR图像配准方法,其特征在于,步骤2具体包括如下子步骤:
(2a)利用下式计算参考图像和预配准图像的幅度互相关矩阵:
L=fftshift(ifft2(fft2(abs(M))×conj(fft2(abs(S)))))
其中,L表示参考图像和预配准图像的幅度互相关矩阵,fft2(·)表示对图像数据做二维傅里叶变换,ifft2(·)表示对图像数据做二维逆傅里叶变换,M表示参考图像数据,S表示预配准图像数据,conj(·)表示取共轭操作,abs(·)表示取矩阵元素的幅度值;fftshift(·)表示将矩阵数据先上半部分和下半部分互换,然后再左半部分和右半部分互换;
(2b)获取所述幅度互相关矩阵L中最大值元素,得到幅度互相关矩阵L中最大值元素对应的坐标位置;
(2c)用幅度互相关矩阵L的中心坐标位置的值减去幅度互相关矩阵L中最大值元素的坐标位置,得到预配准所需的偏移量;
(2d)根据所述预配准所需的偏移量对预配准图像的数据整体做平移,得到平移后的预配准图像,然后截取参考图像中预设像素区域的部分作为主图像,获取平移后的预配准图像中与所述参考图像初步配准的部分作为辅图像。
3.根据权利要求1所述的一种基于图像分割的重复航过极化InSAR图像配准方法,其特征在于,步骤3具体包括如下子步骤:
(3a)对主图像和辅图像的上下左右各预留预设数量的像素;
(3b)对主图像沿距离向和方位向分别进行第一级分块,并对辅图像分别沿距离向和方位向分别进行第一级分块,求解第一级分块后辅图像中每块子图像对应主图像中每块子图像的偏移量,并记录主图像第一级分块中每块子图像的中心点位置坐标;
(3c)对主图像进行第一级分块后的每块子图像再进行下一级分块,按照第一级分块中辅图像的每块子图像对应主图像的每块子图像的偏移量,对辅图像的每块子图像加上偏移量并进行下一级分块,求解下一级分块后辅图像中每块子图像对应主图像中每块子图像的偏移量,并记录主图像下一级分块中每块子图像的中心点位置坐标;
(3d)依此类推,对上一级分块后的子图像再进行下一级分块,直到下一级分块后得到的子图像的像素个数小于设定的像素阈值时停止分级;然后,预设偏移量阈值,得到倒数第二级分块后计算得到的辅图像中子图像的偏移量大于预设偏移量阈值的子图像,并令该子图像的偏移量为零;
(3e)对主图像和辅图像最后一级分块得到的每块子图像去坡度;
(3f)对去坡度后的主图像和去坡度后的辅图像最后一级分块得到的每块子图像进行配准,求解最后一级分块后辅图像中每块子图像对应主图像中每块子图像的偏移量。
4.根据权利要求3所述的一种基于图像分割的重复航过极化InSAR图像配准方法,其特征在于,所述对主图像和辅图像最后一级分块得到的每块子图像去坡度,具体包括:
(1)利用下式求取主图像和辅图像最后一级分块得到的每块子图像的干涉图,
subms=subm.*conj(subs)
其中,subm表示主图像最后一级分块得到的任一子图像,subs表示辅图像最后一级分块得到与所述主图像的任一子图像对应的子图像,conj(·)表示取共轭操作,subms表示主图像和辅图像最后一级分块得到的子图像的干涉图;
(2)对干涉图沿每一个方位向做N1点傅里叶变换,N1大于等于方位向点数,然后分别沿每一个距离向复数相加,得到M1个复数相加后的值,M1为距离向点数,得到M1个复数相加后的值中的最大值A;
对干涉图沿每一个距离向做N2点傅里叶变换,N2大于等于距离向点数,然后分别沿每一个方位向复数相加,得到M2个复数相加后的值,M2为方位向点数,得到M2个复数相加后的值中的最大值B,当2π*A/N或2π*B/N的值大于π时要将其减去2π,保证其值在[-π,π]范围内;
(3)利用下式得到主图像和辅图像最后一级分块得到的每块子图像的二维坡度,
slop=exp(-j*2π*A/N*P1).′*exp(-j*2π*B/N*P2)
其中,j表示虚数单位(·).′表示转置,exp(·)表示取相位值,slop表示子块的二维坡度,P1表示子图像在方位向上的点数,P2表示子图像在距离向上的点数;
(4)将主图像和辅图像中每块子图像的二维坡度与每块子图像的干涉图做内积去除子块的坡度。
5.根据权利要求1所述的一种基于图像分割的重复航过极化InSAR图像配准方法,其特征在于,步骤4具体包括如下子步骤:
(4a)将辅图像中各级分块得到子图像的偏移量与其所属的上一级子图像的偏移量进行累加,从而得到最后一级分块得到的每块子图像的最终偏移量;
(4b)将最后一级分块得到的每块子图像的最终偏移量按照主图像最后一级子图像的中心位置坐标重新排列成矩阵形式,形成偏移量矩阵,对偏移量矩阵做中值滤波;
(4c)对中值滤波后的偏移量矩阵再进行双线性插值,得到辅图像每个像素点对应的偏移量。
CN201510825264.5A 2015-11-24 2015-11-24 基于图像分割的重复航过极化InSAR图像配准方法 Active CN105425216B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201510825264.5A CN105425216B (zh) 2015-11-24 2015-11-24 基于图像分割的重复航过极化InSAR图像配准方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201510825264.5A CN105425216B (zh) 2015-11-24 2015-11-24 基于图像分割的重复航过极化InSAR图像配准方法

Publications (2)

Publication Number Publication Date
CN105425216A CN105425216A (zh) 2016-03-23
CN105425216B true CN105425216B (zh) 2018-02-02

Family

ID=55503544

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201510825264.5A Active CN105425216B (zh) 2015-11-24 2015-11-24 基于图像分割的重复航过极化InSAR图像配准方法

Country Status (1)

Country Link
CN (1) CN105425216B (zh)

Families Citing this family (12)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6351419B2 (ja) * 2014-07-17 2018-07-04 株式会社ミツトヨ 球形状測定方法及び装置
CN106530334B (zh) * 2016-10-21 2019-10-08 北京无线电测量研究所 一种机载干涉合成孔径雷达复图像配准方法及复图像配准系统
CN107085208B (zh) * 2017-04-20 2022-12-13 中国人民解放军海军工程大学 基于分段曲面拟合的干涉合成孔径声纳复图像配准方法
CN106990412B (zh) * 2017-05-11 2019-06-11 哈尔滨工业大学 一种干涉合成孔径激光雷达系统的图像配准方法
CN107945216B (zh) * 2017-11-10 2019-10-11 西安电子科技大学 基于最小二乘估计的多图像联合配准方法
CN108333562A (zh) * 2018-01-30 2018-07-27 西安电子科技大学 一种地形高程自适应的降维图像配准方法
CN108682025B (zh) * 2018-05-23 2022-03-15 东软医疗系统股份有限公司 一种图像配准方法及装置
CN109633639B (zh) * 2018-11-20 2020-07-31 上海无线电设备研究所 Topsar干涉数据的高精度快速配准方法
CN113763295B (zh) * 2020-06-01 2023-08-25 杭州海康威视数字技术股份有限公司 图像融合方法、确定图像偏移量的方法及装置
CN111640147B (zh) * 2020-06-04 2023-10-31 北京无线电测量研究所 用于步进频频带拼接的sar子图像配准方法、装置、计算机设备
CN117616455A (zh) * 2022-06-20 2024-02-27 北京小米移动软件有限公司 多帧图像对齐方法、多帧图像对齐装置及存储介质
CN114820739B (zh) * 2022-07-01 2022-10-11 浙江工商大学 一种面向多光谱相机的图像快速配准方法及装置

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102955157B (zh) * 2011-08-26 2014-04-02 中国科学院空间科学与应用研究中心 一种用于干涉合成孔径雷达图像精配准的快速相干系数法
CN102520407B (zh) * 2011-12-29 2013-06-26 北京理工大学 基于相干系数线性度最优的极化干涉sar图像配准方法

Also Published As

Publication number Publication date
CN105425216A (zh) 2016-03-23

Similar Documents

Publication Publication Date Title
CN105425216B (zh) 基于图像分割的重复航过极化InSAR图像配准方法
Gallup et al. Variable baseline/resolution stereo
CN105300316B (zh) 基于灰度重心法的光条中心快速提取方法
CN102621549B (zh) 多基线/多频段干涉相位解缠频域快速算法
CN102122359B (zh) 一种图像配准方法及装置
CN106683173A (zh) 一种基于邻域块匹配提高三维重建点云稠密程度的方法
CN106093939B (zh) 一种基于相位差统计模型的InSAR图像相位解缠方法
CN103267496B (zh) 一种基于小波变换的改进窗口傅里叶三维测量法
CN104361590A (zh) 一种控制点自适应分布的高分辨率遥感影像配准方法
CN102693216B (zh) 基于分数阶微分的点特征跟踪方法
CN103454636B (zh) 基于多像素协方差矩阵的差分干涉相位估计方法
CN104504740A (zh) 压缩感知框架下的图像融合方法
CN103871039A (zh) 一种sar图像变化检测差异图生成方法
CN103809180B (zh) 用于InSAR地形测量的方位向预滤波处理方法
CN108022259A (zh) 干涉sar复图像配准方法和系统
CN109239710A (zh) 雷达高程信息的获取方法及装置、计算机可读存储介质
CN108573280A (zh) 一种无人船自主通过桥梁的方法
CN102445690A (zh) 合成孔径雷达三维成像的qr分解方法
CN101996415A (zh) 眼球的三维建模方法
CN106408600B (zh) 一种用于太阳高分辨率图像中图像配准的方法
CN105116410B (zh) 基于线性模型匹配的干涉相位图自适应滤波算法
CN103226194A (zh) 一种基于经验模式分解的InSAR干涉相位滤波方法
CN110297242A (zh) 基于压缩感知的合成孔径雷达层析三维成像方法及装置
CN104700401B (zh) 一种基于K‑Means聚类法的图像仿射变换控制点选取方法
CN106056551A (zh) 基于局部相似样例学习的稀疏去噪方法

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant