CN111696207B - 一种基于引导滤波的多基线dem融合方法 - Google Patents

一种基于引导滤波的多基线dem融合方法 Download PDF

Info

Publication number
CN111696207B
CN111696207B CN202010467199.4A CN202010467199A CN111696207B CN 111696207 B CN111696207 B CN 111696207B CN 202010467199 A CN202010467199 A CN 202010467199A CN 111696207 B CN111696207 B CN 111696207B
Authority
CN
China
Prior art keywords
dem
baseline
initialized
local window
detail
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
CN202010467199.4A
Other languages
English (en)
Other versions
CN111696207A (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.)
University of Electronic Science and Technology of China
Original Assignee
University of Electronic Science and Technology of China
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 University of Electronic Science and Technology of China filed Critical University of Electronic Science and Technology of China
Priority to CN202010467199.4A priority Critical patent/CN111696207B/zh
Publication of CN111696207A publication Critical patent/CN111696207A/zh
Application granted granted Critical
Publication of CN111696207B publication Critical patent/CN111696207B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00—Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/05—Geographic models
    • 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
    • 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/904—SAR modes
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/10—Image enhancement or restoration using non-spatial domain filtering
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T5/00—Image enhancement or restoration
    • G06T5/50—Image enhancement or restoration using two or more images, e.g. averaging or subtraction
    • 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
    • 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/20172—Image enhancement details
    • G06T2207/20192—Edge enhancement; Edge preservation
    • 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/20212—Image combination
    • G06T2207/20221—Image fusion; Image merging

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Physics & Mathematics (AREA)
  • Radar, Positioning & Navigation (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Computer Graphics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明公开了一种基于引导滤波的多基线DEM融合方法，该方法首先以先验DEM为基础将各单基线干涉处理得到的DEM进行两尺度分解，先验DEM作为其中的基准层，分解出来的另一部分作为细节层，然后以幅度图作为引导图像，利用引导滤波对单基线对应的细节层DEM进行加权融合，权重由不同基线对应的干涉测高精度来决定，最后将基准层DEM与引导滤波融合后的细节层DEM相加即为多基线融合后的DEM。该方法通过利用先验DEM辅助引导滤波对多基线DEM进行融合，有效地保留了观测场景的地形细节信息，有效提高了多基线InSAR数字高程模型的精度。

Description

一种基于引导滤波的多基线DEM融合方法
技术领域：
本发明属于雷达技术领域，它特别涉及干涉合成孔径雷达(InSAR)高程测绘中多基线 数据融合技术领域。
技术背景：
作为一种工作在微波波段的有源雷达，合成孔径雷达(Synthetic ApertureRadar,SAR)具 有全天时、全天候的成像能力，即无论是白天或黑夜、晴天还是雷雨风雪天气，都可以随时 随地成像，克服了光学和红外系统不能在晚上和复杂天气条件进行成像的缺点。传统的SAR 成像一般只具有二维成像分辨率，在一些起伏比较大的地方比如陡峭的山峰、峡谷以及城市 中矗立挺拔的高楼时，传统SAR成像存在的失真(阴影遮挡效应、空间模糊、顶底倒置等) 导致空间的一些重要信息(比如高度)丢失，所以获取目标场景的数字高程模型(Digital Elevation Model，DEM)是非常有必要的。
干涉合成孔径雷达(Interferometric SAR,InSAR)测量技术是在合成孔径雷达的基础上， 利用两个或者多个位置不同的天线观测同一个目标场景，根据目标到不同天线的斜距差获得 测量数据的干涉相位，再通过平台与地面观测场景的几何关系反演出地面场景的数字高程模 型的技术。目前InSAR已经成为当前提取大面积地表三维图像和地形高程变化信息的一项重 要遥感技术，在地形测绘、自然灾害监测和自然资源调查等领域发挥越来越大的作用。近年 来，随着空间数据基础设施的建设和“数字地球”战略的实施，更加快了DEM与地理信息 系统、遥感等一体化进程，在地理信息系统进行空间分析和决策中，DEM的高精度获取显得 尤为重要。
在干涉测量中，垂直基线决定了InSAR的测高精度，随着垂直基线的增大，数字高程模 型的估计精度提高，但是干涉图像对之间的相干性降低，干涉相位估计的误差增大。相较于 单基线InSAR系统，多基线InSAR测量利用多基线可以构成多个干涉通道，可通过多个视角 分别成像，获取多个基线对应的干涉相位图，丰富测绘场景的干涉相位信息，可优势互补， 在保持相位解缠的稳健性同时，减小干涉相位误差，提高高程估计精度。
传统的多基线InSAR数据融合方式是在相位解缠时，仅利用基线长度的比值作为权重将 不同基线对应的干涉相位融合到最长基线对应的干涉相位上，并不能很有效地融合多基线InSAR数据。因此，为了更好地将多基线InSAR数据融合，保留观测场景的地形细节信息， 提高多基线InSAR数字高程模型的精度，本发明提出了一种基于引导滤波的多基线DEM融 合方法。
发明内容：
本发明公开了一种基于引导滤波的多基线DEM融合方法，该方法首先以先验DEM为基 础将各单基线干涉处理得到的DEM进行两尺度分解，先验DEM作为其中的基准层，分解出来的另一部分作为细节层，然后以幅度图作为引导图像，利用引导滤波对单基线对应的细节 层DEM进行加权融合，权重由不同基线对应的干涉测高精度来决定，最后将基准层DEM与引导滤波融合后的细节层DEM相加即为多基线融合后的DEM。该方法通过利用先验DEM 辅助引导滤波对多基线DEM进行融合，有效地保留了观测场景的地形细节信息，在很大程 度上提高了多基线InSAR数字高程模型的精度。
为了方便描述本发明的内容，首先作以下术语定义：
定义1、合成孔径雷达(SAR)
合成孔径雷达是将雷达固定于载荷运动平台上，结合平台的运动以合成等效阵列以实现 阵列向的分辨率，再利用雷达波束向回波延时实现距离一维成像，从而实现对观测目标二维 成像的一种合成孔径雷达技术。详见文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科 技大学出版社。
定义2、干涉合成孔径雷达(InSAR)
干涉合成孔径雷达是指利用在同一观测场景不同观测角度获得的两幅或者两幅以上SAR 图像进行干涉成像处理，提取干涉相位信息，然后结合雷达系统参数、雷达平台几何位置参 数和观测地形信息反演地形高度及高程变化信息的合成孔径雷达技术。详见文献“合成孔径 雷达成像原理”，皮亦鸣等编著，电子科技大学出版社。
定义3、垂直基线
基线长度是指合成孔径雷达系统中两天线之间的距离，而垂直基线是指实际基线在航迹 法平面内和雷达视线垂直的分量。在本发明中干涉合成孔径雷达系统垂直基线记为B⊥。详见 文献“合成孔径雷达成像原理”，皮亦鸣等编著，电子科技大学出版社。
定义4、单基线干涉合成孔径雷达和多基线干涉合成孔径雷达
干涉合成孔径雷达系统只有一个基线的情况称为单基线干涉合成孔径雷达，有两个或更 多基线的情况称为多基线干涉合成孔径雷达。
定义5、数字高程模型(DEM)
数字高程模型(Digital Elevation Model)，简称DEM，是通过有限的地形高程数据实现 对地面的数字化模拟(即地形表面形态的数字化表达)，它是用一组有序的数值阵列形式表示 地面高程的一种实体地面模型的一个分支，其他各种地形特征值均可由此派生。详见文献 “DEM基于粗数字高程模型信息的干涉相位图生成方法”，郭交等，电子与信息学报， 2010,32(11):2642-2647。
定义6、干涉相位
干涉相位是指将经过配准后的两幅复图像对应的点共轭相乘得到的干涉图。在本发明中 干涉合成孔径雷达系统干涉相位记为
Figure BDA0002513049440000032
详见文献“星载合成孔径雷达干涉成像”，王超等 编著，科学出版社出版。
定义7、解缠相位：
解缠相位是指干涉处理中从缠绕相位经过相位解缠处理后恢复的真实相位。在本发明中 干涉合成孔径雷达系统解缠相位记为φ。详见文献“星载合成孔径雷达干涉测量”，王超等编 著，科学出版社。
定义8、参考斜距
干涉合成孔径雷达系统参考斜距是指干涉合成孔径雷达系统中主天线在合成孔径长度中 间位置到成像空间参考点的距离，在本发明中干涉合成孔径雷达系统参考斜距记为R。
定义9、引导滤波
引导滤波是一种图像滤波技术，通过一张引导图，对初始图像(输入图像)进行滤波处 理，使得最后的输出图像大体上与初始图像相似，但是纹理部分与引导图相似。详见文献“图 像引导滤波的局部多尺度Retinex算法”，方帅等,中国图象图形学报,2012,17(07):748-755。
本发明提出的基于引导滤波的多基线DEM融合方法，它包括以下步骤：
步骤1、初始化基于引导滤波的多基线DEM融合方法所需的参数
初始化基于引导滤波的多基线DEM融合方法所需的参数，包括：多基线InSAR数据组 数，记为L，第i组DEM数据，记为Hi，i＝1,2,…,L；第i组干涉相位，记为
Figure BDA0002513049440000031
i＝1,2,…,L；第i组解缠相位，记为φi，i＝1,2,…,L；第i组垂直基线长度，记为B⊥i，i＝1,2,…,L；第i 组两幅SAR复图像间相干系数，记为γi，i＝1,2,…,L；参考斜距，记为R；雷达入射角，记 为θ；雷达载频波长，记为λ；先验信息构成的基准层地形，记为B；干涉场景原始SAR图 像幅度信息，记为G；引导滤波的局部窗大小，记为ωk；控制引导滤波程度的正则化参数， 记为ε。
步骤2、计算细节层地形
采用公式Di＝Hi-B，计算得到第i组DEM数据的细节层地形，记为Di，i＝1,2,…,L。其中Hi为步骤1中初始化的多基线InSAR第i组DEM数据，B为步骤1中初始化的先验信 息构成的基准层地形。
步骤3、在局部窗中，计算细节地形层的均值，原始SAR幅度图像的均值和方差
采用公式
Figure BDA0002513049440000041
计算得到第i组细节地形层在局部窗ωk中的均值，记为
Figure BDA0002513049440000042
i＝1,2,…,L。采用公式
Figure BDA0002513049440000043
计算得到原始SAR幅度图像在局部窗ωk中的均值，记 为μk。采用公式
Figure BDA0002513049440000044
计算得到原始SAR幅度图像在局部窗ωk中的方差，记为
Figure BDA0002513049440000045
其中
Figure BDA0002513049440000046
为第i组细节地形层在局部窗ωk中第j个像素点的取值，N表示 包含当前像素点j的局部窗个数，L表示步骤1中初始化的多基线InSAR数据组数，ωk表示 步骤1中初始化的局部窗大小，j表示在ωk局部窗中的第j个像素点，Gj为原始SAR幅度图 像在局部窗ωk中的第j个像素点的取值。
步骤4、计算局部窗对应的线性系数
采用公式
Figure BDA0002513049440000047
计算得到第i组细节地形层在局部窗ωk中对应的线性 乘积系数，记为
Figure BDA0002513049440000048
i＝1,2,…,L。采用公式
Figure BDA0002513049440000049
计算得到第i组细节地形层在局 部窗ωk中对应的线性相加系数，记为
Figure BDA00025130494400000410
i＝1,2,…,L。其中Gj为原始SAR幅度图像在局部 窗ωk中的第j个像素点的取值，
Figure BDA00025130494400000411
为第i组细节地形层在局部窗ωk中第j个像素点的取值，
Figure BDA00025130494400000412
为步骤3中的第i组细节地形层在局部窗ωk中的均值，μk为步骤3中的原始SAR幅度图 像在局部窗ωk中的均值，
Figure BDA00025130494400000413
为步骤3中的原始SAR幅度图像在局部窗ωk中的方差，N表示 包含当前像素点j的局部窗个数，L表示步骤1中初始化的多基线InSAR数据组数，ε表示 步骤1中初始化的正则化参数，ωk表示步骤1中初始化的局部窗大小，j表示图像的第j个像 素点。
步骤5、计算所有局部窗对应的线性系数平均值
采用公式
Figure BDA0002513049440000051
计算得到第i组细节地形层在包含了像素点j的所有局部窗对应 的线性乘积系数平均值，记为
Figure BDA0002513049440000052
i＝1,2,…,L。采用公式
Figure BDA0002513049440000053
计算得到第i组细 节地形层在包含了像素点j的所有局部窗对应的线性相加系数平均值，记为
Figure BDA0002513049440000054
i＝1,2,…,L。 其中
Figure BDA0002513049440000055
为第i组细节地形层在包含了像素点j的所有局部窗中第n个局部窗对应的线性乘积 系数，
Figure BDA0002513049440000056
为第i组细节地形层在包含了像素点j的所有局部窗中第n个局部窗对应的线性相加 系数，N表示包含当前像素点j的局部窗个数，L表示步骤1中初始化的多基线InSAR数据 组数，ω(j)表示包含像素点j的所有局部窗集合，n表示在包含了像素点j的所有局部窗中第 n个局部窗。
步骤6、对各单基线细节层DEM进行引导滤波处理
采用公式
Figure BDA0002513049440000057
计算得到第i组经引导滤波处理后的细节层DEM，记为GDi， i＝1,2,…,L。其中
Figure BDA0002513049440000058
为步骤5中的第i组细节地形层在包含了像素点j的所有局部窗对应的 线性乘积系数平均值，
Figure BDA0002513049440000059
为步骤5中第i组细节地形层在包含了像素点j的所有局部窗对应 的线性相加系数平均值，Gj为原始SAR幅度图像在局部窗ωk中的第j个像素点的取值，j表 示图像的第j个像素点，X表示图像的所有像素点集合。
步骤7、计算多基线InSAR融合权重
采用公式
Figure RE-GDA00025811664700000510
计算 得到第i组干涉相位的概率密度函数，记为
Figure RE-GDA00025811664700000511
i＝1,2,…,L。采用公式
Figure RE-GDA00025811664700000512
计算得到第i组干涉相位的估计方差，记为
Figure RE-GDA00025811664700000513
i＝1,2,…,L。 采用公式
Figure RE-GDA00025811664700000514
计算得到第i组高程模糊度，记为Hamb,i，i＝1,2,…,L。采用公式
Figure RE-GDA00025811664700000515
计算得到第i组高度估计标准差，记为σh,i，i＝1,2,…,L。采用公式
Figure RE-GDA00025811664700000516
计算得到第i组DEM的融合权重，记为Wi，i＝1,2,…,L。其中
Figure RE-GDA00025811664700000517
为步骤1 中初始化的第i组干涉相位，φi为步骤1中初始化的第i组解缠相位，B⊥i为步骤1中初始化 的第i组垂直基线长度，γi为步骤1中初始化的第i组两幅SAR复图像间相干系数，R为步 骤1中初始化的雷达参考斜距，θ为步骤1中初始化的入射角，λ为步骤1中初始化的雷达 载频波长，L为步骤1中初始化的多基线InSAR数据组数。
步骤8、重构出基于引导滤波的多基线InSAR融合DEM
采用公式
Figure BDA0002513049440000063
计算得到加权融合后的细节层地形，记为
Figure BDA0002513049440000064
采用公式
Figure BDA0002513049440000065
计算得到多基线InSAR融合DEM，记为DEMfused。其中GDi为步骤6 中的第i组经引导滤波处理后的细节层DEM，Wi表示步骤7中的第i组DEM的融合权重，L 表示步骤1中初始化的多基线InSAR数据组数，B为步骤1中初始化的先验信息构成的基准 层地形。
经过以上步骤，得到基于引导滤波的多基线DEM融合方法处理的数字高程模型DEMfused。
本发明的创新点在于提出了一种基于引导滤波的多基线DEM融合方法，本发明用先验 地形来提取细节层DEM，利用引导滤波对各单基线对应的细节层DEM进行加权融合，能有 效保留地形的细节信息，有效提高了多基线InSAR数字高程模型的测高精度。
本发明的优点是利用先验地形信息，在一定程度上保证地形的整体变化趋势是正确的； 采用基于局部线性模型的引导滤波方法，考虑地形在局部窗内的缓变特性，能有效抑制脉冲 噪声干扰；采用原始SAR图像的幅度信息作为引导，充分利用干涉场景的空间幅度信息，能 够相互补偿低相干区域、阴影叠掩、高误差区域；此外通过对高度值进行加权平均来减少原 始DEM中的随机误差，有效保留地形边缘的细节信息，有效提高多基线InSAR数字高程模 型精度。
附图说明
图1为本发明所提供方法的流程示意框图；
图2为本发明所提供方法的多基线DEM融合方法仿真参数。
具体实施方式
本发明主要采用仿真实验的方法进行验证，所有步骤、结论都在MATLAB R2017b软件 上验证正确。具体实施步骤如下：
步骤1、初始化基于引导滤波的多基线DEM融合方法所需的参数
初始化基于引导滤波的多基线DEM融合方法所需的参数，包括：多基线InSAR数据组 数，记为L＝4，第i组DEM数据，记为Hi，i＝1,2,…,4；第i组干涉相位，记为
Figure BDA0002513049440000071
i＝1,2,…,4； 第i组解缠相位，记为φi，i＝1,2,…,4；垂直基线长度分别为B⊥1＝101m，B⊥2＝64m，B⊥3＝94m， B⊥4＝151m；两幅SAR复图像间相干系数分别为γ1＝0.746，γ2＝0.771，γ3＝0.757，γ4＝0.752； 参考斜距，记为R＝618km；雷达入射角，记为θ＝35.25°；雷达载频波长，记为λ＝31mm；先 验信息构成的基准层地形，记为B；干涉场景原始SAR图像幅度信息，记为G；引导滤波的 局部窗大小，记为ωk＝49；控制引导滤波程度的正则化参数，记为ε＝0.1。
步骤2、计算细节层地形
采用公式Di＝Hi-B，计算得到第i组DEM数据的细节层地形，记为Di，i＝1,2,…,L。其中Hi为步骤1中初始化的多基线InSAR第i组DEM数据，B为步骤1中初始化的先验信 息构成的基准层地形。
步骤3、在局部窗中，计算细节地形层的均值，原始SAR幅度图像的均值和方差
采用公式
Figure BDA0002513049440000072
计算得到第i组细节地形层在局部窗ωk中的均值，记为
Figure BDA0002513049440000073
i＝1,2,…,L。采用公式
Figure BDA0002513049440000074
计算得到原始SAR幅度图像在局部窗ωk中的均值，记 为μk。采用公式
Figure BDA0002513049440000075
计算得到原始SAR幅度图像在局部窗ωk中的方差，记为
Figure BDA0002513049440000076
其中
Figure BDA0002513049440000077
为第i组细节地形层在局部窗ωk中第j个像素点的取值，N＝49 表示包含当前像素点j的局部窗个数，L＝4表示步骤1中初始化的多基线InSAR数据组数， ωk＝49表示步骤1中初始化的局部窗大小，j表示在ωk局部窗中的第j个像素点，Gj为原始 SAR幅度图像在局部窗ωk中的第j个像素点的取值。
步骤4、计算局部窗对应的线性系数
采用公式
Figure BDA0002513049440000078
计算得到第i组细节地形层在局部窗ωk中对应的线性 乘积系数，记为
Figure BDA0002513049440000081
i＝1,2,…,L。采用公式
Figure BDA0002513049440000082
计算得到第i组细节地形层在局 部窗ωk中对应的线性相加系数，记为
Figure BDA0002513049440000083
i＝1,2,…,L。其中Gj为原始SAR幅度图像在局部 窗ωk中的第j个像素点的取值，
Figure BDA0002513049440000084
为第i组细节地形层在局部窗ωk中第j个像素点的取值，
Figure BDA0002513049440000085
为步骤3中的第i组细节地形层在局部窗ωk中的均值，μk为步骤3中的原始SAR幅度图 像在局部窗ωk中的均值，
Figure BDA0002513049440000086
为步骤3中的原始SAR幅度图像在局部窗ωk中的方差，N＝49 表示包含当前像素点j的局部窗个数，L＝4表示步骤1中初始化的多基线InSAR数据组数， ε＝0.1表示步骤1中初始化的正则化参数，ωk＝49表示步骤1中初始化的局部窗大小，j表示 图像的第j个像素点。
步骤5、计算包含了像素点j的所有局部窗对应的线性系数平均值
采用公式
Figure BDA0002513049440000087
计算得到第i组细节地形层在包含了像素点j的所有局部窗对应 的线性乘积系数平均值，记为
Figure BDA0002513049440000088
i＝1,2,…,L。采用公式
Figure BDA0002513049440000089
计算得到第i组细 节地形层在包含了像素点j的所有局部窗对应的线性相加系数平均值，记为
Figure BDA00025130494400000810
i＝1,2,…,L。 其中
Figure BDA00025130494400000811
为第i组细节地形层在包含了像素点j的所有局部窗中第n个局部窗对应的线性乘积 系数，
Figure BDA00025130494400000812
为第i组细节地形层在包含了像素点j的所有局部窗中第n个局部窗对应的线性相加 系数，N＝49表示包含当前像素点j的局部窗个数，L＝4表示步骤1中初始化的多基线InSAR 数据组数，ω(j)表示包含像素点j的所有局部窗集合，n表示在包含了像素点j的所有局部窗 中第n个局部窗。
步骤6、对各单基线细节层DEM进行引导滤波处理
采用公式
Figure BDA00025130494400000813
计算得到第i组经引导滤波处理后的细节层DEM，记为GDi， i＝1,2,…,L。其中
Figure BDA00025130494400000814
为步骤5中的第i组细节地形层在包含了像素点j的所有局部窗对应的 线性乘积系数平均值，
Figure BDA00025130494400000815
为步骤5中第i组细节地形层在包含了像素点j的所有局部窗对应 的线性相加系数平均值，Gj为原始SAR幅度图像在局部窗ωk中的第j个像素点的取值，j表 示图像的第j个像素点，X表示图像的所有像素点集合。
步骤7、计算多基线InSAR融合权重
采用公式
Figure RE-GDA0002581166470000091
计算 得到第i组干涉相位的概率密度函数，记为
Figure RE-GDA0002581166470000092
i＝1,2,…,L。采用公式
Figure RE-GDA0002581166470000093
计算得到第i组干涉相位的估计方差，记为
Figure RE-GDA0002581166470000094
i＝1,2,…,L。 采用公式
Figure RE-GDA0002581166470000095
计算得到第i组高程模糊度，记为Hamb,i，i＝1,2,…,L。采用公式
Figure RE-GDA0002581166470000096
计算得到第i组高度估计标准差，记为σh,i，i＝1,2,…,L。采用公式
Figure RE-GDA0002581166470000097
计算得到第i组DEM的融合权重，记为Wi，i＝1,2,…,L。其中
Figure RE-GDA0002581166470000098
为步骤1中初始化的第i组干涉相位，φi为步骤1中初始化的第i组解缠相位，B⊥i为步骤1中初始化的第i组垂直基线长度，γi为步骤1中初始化的第i组两幅SAR复图像间相干系数，R＝618km为步骤1中初始化的雷达参考斜距，θ＝35.25°为步骤1中初始化的入射角，λ＝31mm为步骤1中初始化的雷达载频波长，L＝4为步骤1中初始化的多基线InSAR数据组数。
步骤8、重构出基于引导滤波的多基线InSAR融合DEM
采用公式
Figure BDA0002513049440000099
计算得到加权融合后的细节层地形，记为
Figure BDA00025130494400000910
采用公式
Figure BDA00025130494400000911
计算得到多基线InSAR融合DEM，记为DEMfused。其中GDi为步骤6 中的第i组经引导滤波处理后的细节层DEM，Wi表示步骤7中的第i组DEM的融合权重， L＝4表示步骤1中初始化的多基线InSAR数据组数，B为步骤1中初始化的先验信息构成的 基准层地形。
经过以上步骤，可得到基于引导滤波的多基线DEM融合方法处理的数字高程模型DEMfused。
经过计算机仿真及实测数据结果证明，本发明通过用先验地形来提取细节层DEM，利用 引导滤波对各单基线对应的细节层DEM进行加权融合，与传统的多基线InSAR数据在相位 解缠时进行融合的方法相比，本发明所提方法极大地减少了原始DEM中的随机误差，有效 保留了地形边缘的细节信息，有效提高了多基线InSAR数字高程模型精度。

Claims (1)

1.一种基于引导滤波的多基线DEM融合方法，其特征是它包括以下步骤：
步骤1、初始化基于引导滤波的多基线DEM融合方法所需的参数
初始化基于引导滤波的多基线DEM融合方法所需的参数，包括：多基线InSAR数据组数，记为L，第i组DEM数据，记为Hi，i＝1,2,…,L；第i组干涉相位，记为
Figure RE-FDA0002581166460000011
第i组解缠相位，记为φi，i＝1,2,…,L；第i组垂直基线长度，记为B⊥i，i＝1,2,…,L；第i组两幅SAR复图像间相干系数，记为γi，i＝1,2,…,L；参考斜距，记为R；雷达入射角，记为θ；雷达载频波长，记为λ；先验信息构成的基准层地形，记为B；干涉场景原始SAR图像幅度信息，记为G；引导滤波的局部窗大小，记为ωk；控制引导滤波程度的正则化参数，记为ε；
步骤2、计算细节层地形
采用公式Di＝Hi-B，计算得到第i组DEM数据的细节层地形，记为Di，i＝1,2,…,L；其中Hi为步骤1中初始化的多基线InSAR第i组DEM数据，B为步骤1中初始化的先验信息构成的基准层地形；
步骤3、在局部窗中，计算细节地形层的均值，原始SAR幅度图像的均值和方差
采用公式
Figure RE-FDA0002581166460000012
计算得到第i组细节地形层在局部窗ωk中的均值，记为
Figure RE-FDA0002581166460000013
i＝1,2,…,L；采用公式
Figure RE-FDA0002581166460000014
计算得到原始SAR幅度图像在局部窗ωk中的均值，记为μk；采用公式
Figure RE-FDA0002581166460000015
计算得到原始SAR幅度图像在局部窗ωk中的方差，记为
Figure RE-FDA0002581166460000016
其中
Figure RE-FDA0002581166460000017
为第i组细节地形层在局部窗ωk中第j个像素点的取值，N表示包含当前像素点j的局部窗个数，L表示步骤1中初始化的多基线InSAR数据组数，ωk表示步骤1中初始化的局部窗大小，j表示在ωk局部窗中的第j个像素点，Gj为原始SAR幅度图像在局部窗ωk中的第j个像素点的取值；
步骤4、计算局部窗对应的线性系数
采用公式
Figure RE-FDA0002581166460000021
计算得到第i组细节地形层在局部窗ωk中对应的线性乘积系数，记为
Figure RE-FDA0002581166460000022
采用公式
Figure RE-FDA0002581166460000023
计算得到第i组细节地形层在局部窗ωk中对应的线性相加系数，记为
Figure RE-FDA0002581166460000024
其中Gj为原始SAR幅度图像在局部窗ωk中的第j个像素点的取值，
Figure RE-FDA0002581166460000025
为第i组细节地形层在局部窗ωk中第j个像素点的取值，
Figure RE-FDA0002581166460000026
为步骤3中的第i组细节地形层在局部窗ωk中的均值，μk为步骤3中的原始SAR幅度图像在局部窗ωk中的均值，
Figure RE-FDA0002581166460000027
为步骤3中的原始SAR幅度图像在局部窗ωk中的方差，N表示包含当前像素点j的局部窗个数，L表示步骤1中初始化的多基线InSAR数据组数，ε表示步骤1中初始化的正则化参数，ωk表示步骤1中初始化的局部窗大小，j表示图像的第j个像素点；
步骤5、计算所有局部窗对应的线性系数平均值
采用公式
Figure RE-FDA0002581166460000028
计算得到第i组细节地形层在包含了像素点j的所有局部窗对应的线性乘积系数平均值，记为
Figure RE-FDA0002581166460000029
采用公式
Figure RE-FDA00025811664600000210
计算得到第i组细节地形层在包含了像素点j的所有局部窗对应的线性相加系数平均值，记为
Figure RE-FDA00025811664600000211
其中
Figure RE-FDA00025811664600000212
为第i组细节地形层在包含了像素点j的所有局部窗中第n个局部窗对应的线性乘积系数，
Figure RE-FDA00025811664600000213
为第i组细节地形层在包含了像素点j的所有局部窗中第n个局部窗对应的线性相加系数，N表示包含当前像素点j的局部窗个数，L表示步骤1中初始化的多基线InSAR数据组数，ω(j)表示包含像素点j的所有局部窗集合，n表示在包含了像素点j的所有局部窗中第n个局部窗；
步骤6、对各单基线细节层DEM进行引导滤波处理
采用公式
Figure RE-FDA00025811664600000214
计算得到第i组经引导滤波处理后的细节层DEM，记为GDi，i＝1,2,…,L；其中
Figure RE-FDA00025811664600000215
为步骤5中的第i组细节地形层在包含了像素点j的所有局部窗对应的线性乘积系数平均值，
Figure RE-FDA00025811664600000216
为步骤5中第i组细节地形层在包含了像素点j的所有局部窗对应的线性相加系数平均值，Gj为原始SAR幅度图像在局部窗ωk中的第j个像素点的取值，j表示图像的第j个像素点，X表示图像的所有像素点集合；
步骤7、计算多基线InSAR融合权重
采用公式
Figure RE-FDA0002581166460000031
计算得到第i组干涉相位的概率密度函数，记为
Figure RE-FDA0002581166460000032
采用公式
Figure RE-FDA0002581166460000033
计算得到第i组干涉相位的估计方差，记为
Figure RE-FDA0002581166460000034
i＝1,2,…,L；采用公式
Figure RE-FDA0002581166460000035
计算得到第i组高程模糊度，记为Hamb,i，i＝1,2,…,L；采用公式
Figure RE-FDA0002581166460000036
计算得到第i组高度估计标准差，记为σh,i，i＝1,2,…,L；采用公式
Figure RE-FDA0002581166460000037
计算得到第i组DEM的融合权重，记为Wi，i＝1,2,…,L；其中
Figure RE-FDA0002581166460000038
为步骤1中初始化的第i组干涉相位，φi为步骤1中初始化的第i组解缠相位，B⊥i为步骤1中初始化的第i组垂直基线长度，γi为步骤1中初始化的第i组两幅SAR复图像间相干系数，R为步骤1中初始化的雷达参考斜距，θ为步骤1中初始化的入射角，λ为步骤1中初始化的雷达载频波长，L为步骤1中初始化的多基线InSAR数据组数；
步骤8、重构出基于引导滤波的多基线InSAR融合DEM
采用公式
Figure RE-FDA0002581166460000039
计算得到加权融合后的细节层地形，记为
Figure RE-FDA00025811664600000310
采用公式
Figure RE-FDA00025811664600000311
计算得到多基线InSAR融合DEM，记为DEMfused；其中GDi为步骤6中的第i组经引导滤波处理后的细节层DEM，Wi表示步骤7中的第i组DEM的融合权重，L表示步骤1中初始化的多基线InSAR数据组数，B为步骤1中初始化的先验信息构成的基准层地形；
经过以上步骤，得到基于引导滤波的多基线DEM融合方法处理的数字高程模型DEMfused。
CN202010467199.4A 2020-05-28 2020-05-28 一种基于引导滤波的多基线dem融合方法 Active CN111696207B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202010467199.4A CN111696207B (zh) 2020-05-28 2020-05-28 一种基于引导滤波的多基线dem融合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202010467199.4A CN111696207B (zh) 2020-05-28 2020-05-28 一种基于引导滤波的多基线dem融合方法

Publications (2)

Publication Number Publication Date
CN111696207A CN111696207A (zh) 2020-09-22
CN111696207B true CN111696207B (zh) 2022-10-11

Family

ID=72478728

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202010467199.4A Active CN111696207B (zh) 2020-05-28 2020-05-28 一种基于引导滤波的多基线dem融合方法

Country Status (1)

Country Link
CN (1) CN111696207B (zh)

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
ITPI20070012A1 (it) * 2007-02-14 2008-08-15 Univ Pisa Metodo di elaborazione di dati radar multipassaggio per il rilevamento e l'analisi di componenti multiple di retrodiffusori non stazionari
CN101551455A (zh) * 2009-05-13 2009-10-07 西安电子科技大学 干涉合成孔径雷达三维地形成像系统及其高程测绘方法
CN101893710A (zh) * 2009-05-20 2010-11-24 中国科学院电子学研究所 一种非均匀分布的多基线合成孔径雷达三维成像方法
CN103424744A (zh) * 2012-05-16 2013-12-04 中国科学院电子学研究所 干涉sar叠掩区域数字高程模型重建的方法
CN107102333A (zh) * 2017-06-27 2017-08-29 北京航空航天大学 一种星载InSAR长短基线融合解缠方法
CN109242872A (zh) * 2018-08-27 2019-01-18 西安电子科技大学 基于srtm dem的干涉基线估计方法
KR20190043478A (ko) * 2017-10-18 2019-04-26 서울대학교산학협력단 고해상도 수치 표고 모델 생성 시스템 및 방법
CN109993717A (zh) * 2018-11-14 2019-07-09 重庆邮电大学 一种结合引导滤波和ihs变换的遥感图像融合方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
ITCS20070012A1 (it) * 2007-02-14 2008-08-15 Develpack Srl Foglio per cucinare in forno a microonde e in forno tradizionale, per conservare in frigo, per la consumazione o per l asporto di prodotti vari

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
ITPI20070012A1 (it) * 2007-02-14 2008-08-15 Univ Pisa Metodo di elaborazione di dati radar multipassaggio per il rilevamento e l'analisi di componenti multiple di retrodiffusori non stazionari
CN101551455A (zh) * 2009-05-13 2009-10-07 西安电子科技大学 干涉合成孔径雷达三维地形成像系统及其高程测绘方法
CN101893710A (zh) * 2009-05-20 2010-11-24 中国科学院电子学研究所 一种非均匀分布的多基线合成孔径雷达三维成像方法
CN103424744A (zh) * 2012-05-16 2013-12-04 中国科学院电子学研究所 干涉sar叠掩区域数字高程模型重建的方法
CN107102333A (zh) * 2017-06-27 2017-08-29 北京航空航天大学 一种星载InSAR长短基线融合解缠方法
KR20190043478A (ko) * 2017-10-18 2019-04-26 서울대학교산학협력단 고해상도 수치 표고 모델 생성 시스템 및 방법
CN109242872A (zh) * 2018-08-27 2019-01-18 西安电子科技大学 基于srtm dem的干涉基线估计方法
CN109993717A (zh) * 2018-11-14 2019-07-09 重庆邮电大学 一种结合引导滤波和ihs变换的遥感图像融合方法

Also Published As

Publication number Publication date
CN111696207A (zh) 2020-09-22

Similar Documents

Publication Publication Date Title
Massonnet et al. Radar interferometry: limits and potential
CN109633648B (zh) 一种基于似然估计的多基线相位估计装置及方法
WO2016086699A1 (zh) 一种结合局部频率估计的小波域InSAR干涉相位滤波方法
CN109212522B (zh) 一种基于双基星载sar的高精度dem反演方法及装置
CN108663678B (zh) 基于混合整数优化模型的多基线InSAR相位解缠算法
CN108132468B (zh) 一种多基线极化干涉sar建筑物高度提取方法
CN109597074B (zh) 一种sar影像几何定位参数校正方法及系统
CN111856459B (zh) 一种改进的DEM最大似然约束多基线InSAR相位解缠方法
CN112882030B (zh) InSAR成像干涉一体化处理方法
CN109239710B (zh) 雷达高程信息的获取方法及装置、计算机可读存储介质
CN108008383A (zh) 一种迭代多基线高精度四次fft相位解缠方法
CN111696207B (zh) 一种基于引导滤波的多基线dem融合方法
CN110618409B (zh) 顾及叠掩及阴影的多通道InSAR干涉图仿真方法及系统
CN109946682B (zh) 基于ICESat/GLAS的GF3数据基线估计方法
Deo et al. Evaluation of interferometric SAR DEMs generated using TanDEM-X data
Agrawal et al. Accuracy assessment of digital elevation model generated by SAR stereoscopic technique using COSMO-SkyMed data
CN115546264A (zh) 一种星载InSAR图像精细配准与立体测量方法
CN115540908A (zh) 基于小波变换的InSAR干涉条纹匹配方法
Usman et al. Comparative analysis on digital surface model of urban area from Sentinel-1 SAR interferometry and aerial photogrammetry for disaster mitigation plan
CN112946647A (zh) 大气误差改正InSAR干涉图堆叠地质灾害普查方法和装置
Saied et al. Digital Elevation Model Enhancement using CNN-Based Despeckled SAR Images
Nascetti et al. Radargrammetric digital surface models generation from high resolution satellite SAR imagery: Methodology and case studies
Zhu et al. Towards global 3d/4d urban modeling using tandem-x data
Shiping DEM generation using ERS-1/2 interferometric SAR data
Lombardi et al. Accuracy of high resolution CSK interferometric Digital Elevation Models

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