CN104122553A - 一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法 - Google Patents
一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法 Download PDFInfo
- Publication number
- CN104122553A CN104122553A CN201410353927.3A CN201410353927A CN104122553A CN 104122553 A CN104122553 A CN 104122553A CN 201410353927 A CN201410353927 A CN 201410353927A CN 104122553 A CN104122553 A CN 104122553A
- Authority
- CN
- China
- Prior art keywords
- frame
- phase
- deformation
- track
- coherent
- 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
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
- G01C—MEASURING DISTANCES, LEVELS OR BEARINGS; SURVEYING; NAVIGATION; GYROSCOPIC INSTRUMENTS; PHOTOGRAMMETRY OR VIDEOGRAMMETRY
- G01C5/00—Measuring height; Measuring distances transverse to line of sight; Levelling between separated points; Surveyors' levels
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Electromagnetism (AREA)
- Computer Networks & Wireless Communication (AREA)
- Radar Systems Or Details Thereof (AREA)
Abstract
一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法,它有七大步骤:步骤一:长条带雷达数据分块与配准;步骤二:差分干涉纹图集生成;步骤三:轨道残余相位去除;步骤四:CTInSAR时序分析;步骤五:长条带相邻分块检验与校正;步骤六:多轨道坐标系统统一;步骤七:多轨道信息检验与参考基准统一。本发明解决了InSAR在大区域形变调查与监测中遇到的两个主要问题,一是如何减弱形变信息提取过程中基线误差、大气波动以及相位解缠噪声等在空间域随范围增大而表现明显的低频特性的误差影响;二是如何将相邻轨道的形变监测结果在空间位置上进行精确配准,并统一不同轨道间形变测量的参考基准。
Description
技术领域
本发明涉及一种集成多轨道、长条带CTInSAR(Coherent Target SAR Interferometry,相干目标合成孔径雷达干涉测量)的区域性地面沉降监测方法,属于合成孔径雷达干涉测量(InSAR)技术应用领域。它克服了大区域CTInSAR数据处理过程中误差传播与基准偏差的问题,达到了区域性地面沉降的高精度调查与监测。
背景技术
当前我国存在大范围的地面沉降,华北平原、长江三角洲、汾渭谷地三大主要沉降区的影响范围超过数万平方公里,且有进一步扩展连片分布的特点。其他地区如江汉平原、松嫩平原、下辽河以及珠江三角洲等也有不同程度的地面沉降发生。传统的主要地面沉降监测手段(水准测量、基岩标、分层标和GPS监测网等)均面临空间覆盖有限、监测频率低、监测点密度低、观测时间长且难以同步等问题。星载InSAR技术是根据星载雷达图像的相位数据来获取地面形变信息的一门技术,能满足调查监测精度要求的同时,具有大覆盖范围、快速重访、高监测点位密度等特点,且易于组织,人为干扰少。
大区域地面沉降场调查与监测要求雷达数据具备二维方向上的完整覆盖,而单个图幅重复轨InSAR主要解决局部地区(几百至几千平方公里)地表形变场的高精度连续监测,以常用的ENVISAT ASAR扫描模式的数据来说,单个图幅覆盖范围在100km×100km,涉及几万平方公里大区域的地面沉降监测必然需要多轨道、长条带的InSAR成果集成。多轨道集成包含相邻平行轨道集成(如图1所示)以及升降轨模式集成,其中相邻平行轨道集成是解决大范围覆盖的主要技术途径。
多轨道、长条带集成方法基于单图幅CTInSAR时序分析技术,实现大区域地面沉降监测。但随着空间覆盖范围的增大,InSAR数据处理过程中误差因素也随之增大。这些误差主要包括基线误差、大气波动以及相位解缠噪声(相位不连续)等在空间域表现为低频特性的因素。随着覆盖范围的增大,低频信号的累加效应表现明显,影响了沉降场监测结果的准确性,往往使得沉降场观测结果中附加某种趋势信号。
另一个问题是不同轨道(入射角)下InSAR成果具有不同的坐标系和参考基准。无论采用何种时序分析方法,每个轨道下的形变监测结果都存在于独立的坐标系中(雷达坐标系),其形变场参考基准也不相同。不同轨道下SAR数据均与各自独立的主图像进行精确配准,因而,其观测结果位于不同主图像的坐标系内。斜距成像的特点决定了不同入射的雷达图像存在不同的变形(扭曲、旋转等),特别是在地形起伏较大的地区。即使同一景雷达图像内,其象元的分辨率也存在差别。此外,InSAR技术所获取的形变量均是相对于某一个特定的参考基准而言的,其参考基准是在时序分析过程中选择的某一特定相干目标,其活动状况往往未知。这使得不同轨道下的形变场存在参考基准偏差,涉及到整体偏差补偿问题。
因此,为解决区域性地面沉降InSAR监测中存在的问题,本专利通过将长条带雷达数据分块,各分块独立完成CTInSAR时序分析,获取各分块高精度地面沉降监测成果,有效地降低了长条带合成孔径雷达数据的干涉处理误差,并通过对同一条带不同分块和多轨道不同条带下形变监测结果的检验与校正,最终达到大区域地面沉降的全覆盖、同时段调查与监测。
发明内容
1、目的
本发明的目的是为了获取全覆盖、同时段的高精度区域性地面沉降信息,解决了InSAR在大区域形变调查与监测中遇到的两个主要问题,一是如何减弱形变信息提取过程中基线误差、大气波动以及相位解缠噪声等在空间域随范围增大而表现明显的低频特性的误差影响;二是如何将相邻轨道的形变监测结果在空间位置上进行精确配准,并统一不同轨道间形变测量的参考基准。
2、技术方案
本发明是一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法,建立了满足低相干条件下区域性地面沉降监测的CTInSAR技术,方法具体步骤如下:
步骤一:长条带雷达数据分块与配准
符合InSAR测量要求的星载SAR数据通常以长条带方式获取,数据通常为单视复数数据(single look complex digital image,SLC),记录了目标在影像上的散射强度以及雷达回波的相位信息。
假设同一长条带获得m个重轨数据,在对长条带数据进行处理时,首先选择一幅主影像M,对其进行方位向上的分块裁切,得到n个分块分别为Mframe1、Mframe2、……、Mframen,每个分块称为一个Frame,每个Frame在雷达影像的方位向尽量与距离向保持相同地面距离范围,相邻分块之间保留20%的重叠区域。
将其他SAR数据做为辅影像分别以Mframe1、Mframe2、……、Mframen为主影像进行精密配准,得到n个frame配准后的SLC数据集[M,S1,S2,……,Sm-1]。针对每个frame的SLC数据集执行步骤二至步骤四。
步骤二:差分干涉纹图集生成
考虑空间基线与时间基线对SAR图像对造成的失相干影响,选择具有较短空间基线与时间基线的图像对两两进行组合,生成差分干涉纹图集。
每个图像对主副影像共轭相乘得到干涉纹图,干涉纹图的相位组成包括:
φint=δφflat+δφtopo+δφmov+δφatm+δφnoise (1)
式中:δφflat为斜距变化产生的平地相位,δφtopo为地形相位,δφmov是两次获取SAR图像时地表移动引起的沿雷达视线向形变相位,δφatm为大气波动引起的相位变化,δφnoise为噪声相位。
根据已知参数模拟平地相位与地形相位:
其中:λ为雷达波长,r为雷达到目标之间的斜距,b为垂直基线长,Δr和Δh分别为像元间的斜距和高程增量。
去除平地、地形相位得到差分干涉相位:
φdiff=δφmov+δφtoperror+δφatm+δφnoise (4)
式(4)中δφtoperror为差分干涉纹图中由于外部DEM误差带来的相位。
步骤三:轨道残余相位去除
地形相位与平地相位的计算均与基线长度相关,而轨道的精度直接影响干涉相对空间基线的计算。由于轨道的不准确性形成的残余相位,即Orbit ramp或Orbit residual。应用精密轨道参数修正初始轨道参数可以有效降低这种因素的影响,但当干涉图范围扩大时,二维平面上轨道残余相位的积累将使得形变图上附加一种趋势,这种趋势随着数据覆盖范围的增大而变化。图2所示为沉降场附加轨道残余相位的剖面。
本发明中基于最小二乘进行趋势面拟合进而去除轨道残余相位。首先依据先验知识挑选出差分干涉纹图集中有明显残余相位的干涉纹图,针对有误差的干涉纹图首先进行相位解缠,通过提取解缠相位图中一定采样间隔的解缠相位,利用二次或者高次函数进行相位回归,进而拟合趋势面。这里选择二次函数(式5)作为拟合函数进行相位回归。
式中:为采样点解缠相位;(x,y)为采样点图像坐标;ai为多项式系数。基于最小二乘原则,计算出ai,从差分相位φdiff中减去拟合的轨道残余相位。
步骤四:CTInSAR时序分析
当SLC数据集中数据量较少时(m<20),依据相干目标频谱变化小这一特性进行相干目标的选择。点目标的实际几何尺寸往往小于雷达分辨单元,与分布目标具有不同的信号特征,不具备斑点噪声的特性,相位稳定,因而可将其视为相干目标。由于点目标稳定的后向散射特性,在获取原始信号过程中各个子视图像上表现出相似的后向散射特征,不同距离向和方位向的子带宽内的散射强度大致相同。因而可将SLC图像作反变换,生成若干个子视图(sub-looks),将子视图作谱相关分析,根据其相关性的大小来识别点目标。假设各子视内的某像元的信号xi(n)=aig(n),1≤i≤L为视数,g(n)为该像元信号,按照能量比标准衡量L个子视信号的整体相似性为:
当a1=a2=…=aL时,E=1。能量比越大,整体相似性越好,每幅SLC中频谱特征统计值E高于某一给定阈值的像元都被识别为点目标候选点,在候选点中将雷达目标的散射强度值高于某一阈值的点作为相干点。
利用相干点将去除了轨道残余相位的差分干涉纹图集的相位进行提取,得到差分干涉相位的点集。相位组成如式(4)所示,其中形变相位δφmov包含线性形变δφlinear和非线性形变δφnonlinear两部分。差分干涉相位可以表达为
式(7)中,λ为雷达波长,v为线性形变速率,T为生成差分干涉纹图的图像对时间基线,r为雷达到目标之间的斜距,b为垂直基线长,ε为高程误差,θ为雷达入射角。
由于在整个干涉纹图集中每个相干点的线性速率和DEM误差均为常量,是模型中的线性分量,利用以下二维模型可以很好地实现形变速率和高程改正的估计:
Ti为第i个干涉像对的时间基线。对于相干目标,由于差分干涉纹图中存在不同的相位补偿,单个像元的相位值在处理过程中并没有实际意义。在确定相位解缠参考点后连接相邻两个相干目标(xm,ym)和(xn,yn),通过式(9)求时序模型函数的γmodel最大值可以得到式(8)中v和ε的估算结果vestimated和εestimated。
式(9)中N为干涉纹图集中包含的干涉纹图个数。
去除线性分量,模型残余相位可用下式表示:
式中:δφdiff为初始差分相位,δφlinearestmodel为线性模型将线性形变和DEM误差合成得到的相位。模型残余相位中包含大气相位影响δφatmos和非线性移动δφnonliear(x,y,Ti)以及失相干等因素造成的噪声δφnoise(x,y,Ti)。
利用1km×1km二维移动平均窗口进行空间低通滤波,滤波完成后,残余相位的噪声δφnoise(x,y,Ti)被滤除,最后,利用时域高通滤波估计大气贡献δφatm,区分出非线性形变相位δφnonliear(x,y,Ti)。
利用前面估计出的线性分量与非线性分量获得完整的形变演变量。
依据上述相干目标时序分析方法,通过相干目标时序分析模型迭代补偿DEM误差、通过滤波减小大气相位以及噪声的影响,获取到各相干目标点的地面沉降量。
步骤五:长条带相邻分块检验与校正
对长条带中的每个frame重复步骤二至步骤四,得到条带下各分块地面沉降信息。长条带的每个frame获取的形变信息图的参考基准位置不同,需要对连续frame的参考基准进行整体偏差补偿。整体偏差估算依赖于重叠部分相干目标形变参数的统计分析,依据同名点按照式(12)统计连续frame之间形变速率的整体偏差,参与统计的相干目标均为满足相关模型的点。在此过程中选择连续frame中的一个作为基准,其余均相对于该frame进行基准偏差补偿。一般情况下,在每个frame时序分析过程中,根据时间基线的变化关系,将在不同时间间隔下稳定区域内的相干目标作为参考基准,以此降低frame之间形变参数的整体偏差。以两块连续frame之间的检验、校正为例,选择其中一个frame作为基准(以先验知识为基础进行主frame的选择),依据frame之间重叠区域同名点的统计信息进行基准偏差的整体补偿
式(12)和(13)中:Δρoff为frame之间基准偏差;ρmi为基准frame重叠区中相干目标i的形变量;ρsi为需调整frame重叠区中相干目标i的形变量;N表示两个frame包含的同名点个数;为调整的frame中所有相干目标改正后的形变向量。
对于重叠区域内的各相干点目标形变量ρi的取值采用距离加权平均方法求解:
式中Pmi和Psi分别为基准frame和调整frame第i个相干目标的权,ρmi和分别为重叠区内第i个相干目标在基准frame和调整frame中的形变量,D为重叠区长度,Dmi和Dsi分别为基准frame和调整frame中相干目标距离重叠区内各自影像方位向边缘的距离(图3a、b)。
步骤六:多轨道坐标系统统一
通过对不同轨道下的多个长条带InSAR地面沉降信息集成进而获取大覆盖范围的监测成果,首先需要对不同的轨道数据进行坐标系统一,具有相同的空间坐标基准。
首先需要获取研究区的DEM数据,根据外部DEM将雷达图像与监测成果进行正射校正,即地理编码,使得不同轨道下的雷达图像位于相同的地面坐标系下,消除因地形起伏以及成像几何关系引起的畸变等因素的影响。高程精度优于10m的DEM数据保证了平坦地区配准精度优于一个像元。依据卫星参数将DEM模拟至SAR坐标系,建立初始坐标变换表,再将DEM模拟的影像与SAR图像进行精密配准改进坐标变换表。通过反向变换,将雷达坐标系下各轨道地面形变信息地理编码至DEM所在的地面坐标系内。
步骤七:多轨道信息检验与参考基准统一
位于相同地面坐标系的不同轨道形变信息集成过程中,受制于雷达成像入射角的影响,不能像长条带连续frame重叠区内的相干目标在相邻frame中基本一致,相邻轨道下相干目标是完全不同的(图4)。
选择其中一条轨道为主轨道,相邻另一条为辅轨道,首先通过反距离加权插值方法将辅轨道上相干目标分析结果进行空间域插值为面,设辅轨道重叠区内共有n个相干目标点,其插值后平面各点坐标(x,y)的沉降速率值通过下式计算得出:
式中ρi为第i个相干目标点的沉降量,di为点(x,y)到第i个相干目标点的距离。再利用主轨道上重叠区内相干点目标位置提取其在辅轨道上插值面上的形变值ρs。
根据相同相干点目标在主轨道的沉降值ρm和在辅轨道的沉降值ρs进行检验与校正,检验与校正方法同长条带相邻分块的检验与校正。
3、优点及功效:本发明通过对长条带、多轨道CTInSAR分析处理,达到了大区域地面沉降的高精度调查与监测。
(1)长条带数据的分块处理,减小了雷达数据常规干涉测量中轨道、大气、相位解缠等误差的传播影响。针对每个分块利用相干目标时序分析方法获取高精度形变信息。
(2)利用长条带连续frame之间重叠区内的地面沉降信息,完成多frame之间的检验与校正,保证了长条带地面沉降信息的完整性与正确性。
(3)通过对长条带数据进行地理编码达到了不同轨道数据的空间基准统一。
(4)利用相邻轨道长条带间重叠区内的相干目标点,通过插值的方法,完成了多轨道的形变信息检验与校正,最终实现了对大区域地面沉降的高精度监测。本方法对我国经济快速发展过程中区域性地面沉降调查与监测,具有普适性。
附图说明
图1:大区域地面沉降InSAR多轨道测量。多轨道长条带数据覆盖大区域范围,通过对这些数据的处理达到大区域地面沉降调查与监测的需要。
图2:残余相位对形变场的影响。沿误差传播方向做剖面,可以看到残余相位为形变量求解带来的趋势性影响,当干涉图范围扩大时,二维平面上轨道残余相位的积累将使得形变图上附加一种明显的趋势,这种趋势随着数据覆盖范围的增大而增大。
图3a:为长条带相邻frame重叠区内相干目标确定权重的距离示意图,Dm为相干目标距离基准frame重叠区边缘的距离,Ds为相干目标距离调整frame重叠区边缘的距离。
图3b:为相邻轨道长条带重叠区内相干目标确定权重的距离示意图,Dm为相干目标距离主轨道重叠区边缘的距离,Ds为相干目标距离辅轨道重叠区边缘的距离。
图4:多轨道相邻长条带重叠区内相干目标位置差异示意图。对于相邻的条带,由于雷达成像入射角的影响,相邻轨道重叠区内相干目标位置是完全不同的,图中十字符号为主轨道相干目标,圆圈符号为辅轨道相干目标。
图5:轨道残余相位去除示意图
图6:Track-218与Track-447重叠区相干目标分布示意图。
图7:Track-218与Track-447重叠区相干目标沉降量互差统计直方图。
图8:多轨道、长条带CTInSAR监测成果示意图。图为京津沉降区2008年度InSAR监测地面沉降速率图,图中散点(相干目标)由浅到深表明沉降愈加严重。
图9:本发明数据处理流程图。
具体实施方式
以京津主要沉降区为例,说明本发明在实施中的具体步骤。
实施过程中利用完整覆盖京津主要沉降区的ENVISAT雷达数据(IS2,strip mode),涉及三个轨道,Track-218、Track-447和Track-175,三个轨道的相邻关系如图1所示,每个条带的处理长度均为200km,宽度为100km,条带间横向重叠度约为40%,获取的雷达数据时间从2007年12月截止2009年1月。所用的雷达数据获取时间如表1:
表1雷达数据列表
Track-218 | Track-447 | Track-175 | |
1 | 20071219 | 20080104 | 20080120 |
2 | 20080123 | 20080208 | 20080330 |
3 | 20080227 | 20080314 | 20080504 |
4 | 20080402 | 20080418 | 20080713 |
5 | 20080507 | 20080523 | 20080817 |
6 | 20080611 | 20080627 | 20080921 |
7 | 20080716 | 20080801 | 20081026 |
8 | 20080820 | 20080905 | 20090104 |
9 | 20080924 | 20081010 | |
10 | 20081029 | 20081219 | |
11 | 20090107 | 20090123 |
实施过程利用的外部DEM为SRTM DEM,轨道参数为doris精密轨道。
见图9,本发明是一种集成多轨道、长条带CTInSAR,的区域性地面沉降监测方法,该方法具体步骤如下:
步骤一:长条带雷达数据分块与配准
长条带雷达数据配准与分块以Track-218为例,根据获取数据的实际长度,以Track-218获取的第一景影像20071219为主影像,将其裁切为2个frame Mpatch1和Mpatch2,两个Frame在雷达影像的方位向与距离向分别保留25000和5000像元,多视后方位向与距离向的宽度均为5000像元,距离约为100km,两个frame之间保留20%的重叠区域。
将Track218下的其他SAR数据分别以Mpatch1、Mpatch2为主影像进行精密配准,得到2个分块配准后的SLC数据集。
针对三个条带下的每个分块执行步骤二至步骤四。
步骤二:差分干涉纹图集生成
首先利用DORIS精密轨道进行SLC数据轨道修正。选择空间基线小于350m,时间基线大于60天小于400天的图像对进行组合,Track-218、Track-447和Track-175三个轨道的各frame分别形成35个、34个和21个干涉图像对组合,每个图像对主副影像共轭相乘得到干涉纹图集。根据已知参数及式(2)计算平地相位,从干涉纹图中去除平地相位。依据卫星参数将STRM DEM模拟至SAR影像坐标系,建立初始坐标变换表,再将DEM模拟SAR影像与SAR图像通过地形特征进行精密配准,进而改进坐标变换表,将STRM DEM通过坐标变换表变换到SAR影像坐标系,利用式(3)计算地形相位。从干涉纹图减去平地相位与地形相位得到差分干涉纹图。
步骤三:轨道残余相位去除
依据先验知识判断轨道残余相位的影响,确定需要进行残余相位去除的差分干涉纹图。Track-218下的干涉纹图集中有24个干涉纹图进行了残余相位的去除。对于每个需要进行残余相位去除的差分干涉纹图进行相位解缠,通过提取解缠相位图中一定采样间隔的解缠相位,依据式(3)利用二次函数进行相位回归,拟合趋势面。从差分相位图中减去拟合的趋势面,得到修正了轨道残余相位差分干涉纹图,利用去除了残余的差分干涉纹图进行下一步处理将避免产生附加的形变趋势,如图2。如图5为轨道残余相位去除效果示意。
步骤四:CTInSAR时序分析
依据点目标识别方法,将E阈值取0.4,强度值阈值取图像平均强度的10倍,得到点集。用点位提取差分干涉相位得到差分干涉相位点集pdiff,Track-218各分块的pdiff就是包含35个干涉图像对相位的点集,针对pdiff进行相干目标时序分析,获取到相干目标点的地面沉降量(包含线性与非线性沉降)。
对3个轨道下每2个分块执行步骤二至步骤四,得到3个条带中每个分块的所覆盖区域的地面沉降信息。
步骤五:长条带相邻分块检验与校正
对3个轨道下的每2个frame的解算参考基准进行整体偏差补偿。以Track-218为例,其中第一个frame包含北京市区,依照先验知识,北京市中心稳定,在相干目标时序分析时选择市中心一点做为解算起点,将这个frame做为基准frame。参考基准整体偏差估算依据重叠部分相干目标形变参数的统计分析,依据同名点按照式(12)统计两个frame之间形变速率的整体偏差为5.26365mm,将位于下方的frame所有相干目标的形变量调整-5.26365mm,对于重叠区内的各相干点目标利用公式(14)进行计算,各相干点权重确定如图3a,得到Track-218的长条带地面沉降信息。
步骤六:多轨道坐标系统统一
通过坐标变换表将3个轨道下的长条带SAR图像及地面沉降信息反向变换到STRMDEM所在的WGS-84坐标系内。
步骤七:多轨道信息检验与参考基准统一
以Track-218为基准调整Track-447,Track-218的北京地区和Track-447的天津地区北部的重叠部分,覆盖范围约为40km×100km,包含相干目标167851个,如图6,放大到局部,两相邻轨道的相干目标差异如图4所示。依据式(16)对Track-447的相干目标插值生成面数据,再由Track-218位于重叠区的相干目标进行沉降量提取,得到重叠区相干目标沉降量差值的统计分布如图7所示,显示差值均值为-0.09408mm,差值范围为[-13.0496mm,12.760325mm],位于[-5mm,5mm]以内的相干目标占总体样本的96.63%。将Track-447下所有相干点目标的沉降量值减去0.09408mm,对于Track-218与Track-447重叠区内的相干点目标取值依式(14)进行计算,各相干点权重确定如图3b。
同样的方法以调整后的Track-447为基准调整Track-175轨道下的沉降量信息,最终得到京津主要沉降区2007年12月到2009年01月的全覆盖、相同时间段的沉降量信息,如图8为京津主要沉降多轨道长条带InSAR监测地面沉降结果图。具体实施过程中获取了天津地区和北京地区的地面沉降水准测量值,视水准观测为沉降量的最或然值时,通过对两组观测值的检验,得到InSAR测量获取的形变量的点位中误差为±5mm。
由此,本发明的方法满足区域性地面沉降调查与监测的测量要求,能够为我国区域性地面沉降的防治工作提供可靠的基础数据。
Claims (1)
1.一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法,其特征在于:该方法具体步骤如下:
步骤一:长条带雷达数据分块与配准
符合InSAR测量要求的星载SAR数据通常以长条带方式获取,数据通常为单视复数数据即SLC,记录了目标在影像上的散射强度以及雷达回波的相位信息;
设同一长条带获得m个重轨数据,在对长条带数据进行处理时,首先选择一幅主影像M,对其进行方位向上的分块裁切,得到n个分块分别为Mframe1、Mframe2、……、Mframen,每个分块称为一个Frame,每个Frame在雷达影像的方位向尽量与距离向保持相同地面距离范围,相邻分块之间保留20%的重叠区域;
将其他SAR数据做为辅影像分别以Mframe1、Mframe2、……、Mframen为主影像进行精密配准,得到n个frame配准后的SLC数据集[M,S1,S2,……,Sm-1];针对每个frame的SLC数据集执行步骤二至步骤四;
步骤二:差分干涉纹图集生成
考虑空间基线与时间基线对SAR图像对造成的失相干影响,选择具有较短空间基线与时间基线的图像对两两进行组合,生成差分干涉纹图集;
每个图像对主副影像共轭相乘得到干涉纹图,干涉纹图的相位组成包括:
φint=δφflat+δφtopo+δφmov+δφatm+δφnoise (1)
式中:δφflat为斜距变化产生的平地相位,δφtopo为地形相位,δφmov是两次获取SAR图像时地表移动引起的沿雷达视线向形变相位,δφatm为大气波动引起的相位变化,δφnoise为噪声相位;
根据已知参数模拟平地相位与地形相位:
其中:λ为雷达波长,r为雷达到目标之间的斜距,b为垂直基线长,Δr和Δh分别为像元间的斜距和高程增量;
去除平地、地形相位得到差分干涉相位:
φdiff=δφmov+δφtoperror+δφatm+δφnoise (4)
式(4)中δφtoperror为差分干涉纹图中由于外部DEM误差带来的相位;
步骤三:轨道残余相位去除
地形相位与平地相位的计算均与基线长度相关,而轨道的精度直接影响干涉相对空间基线的计算;由于轨道的不准确性形成的残余相位,即Orbit ramp或Orbit residual;应用精密轨道参数修正初始轨道参数能有效降低这种因素的影响,但当干涉图范围扩大时,二维平面上轨道残余相位的积累将使得形变图上附加一种趋势,这种趋势随着数据覆盖范围的增大而变化;基于最小二乘进行趋势面拟合进而去除轨道残余相位,首先依据先验知识挑选出差分干涉纹图集中有明显残余相位的干涉纹图,针对有误差的干涉纹图首先进行相位解缠,通过提取解缠相位图中一定采样间隔的解缠相位,利用二次或者高次函数进行相位回归,进而拟合趋势面;这里选择二次函数式(5)作为拟合函数进行相位回归;
式中:为采样点解缠相位;(x,y)为采样点图像坐标;ai为多项式系数;基于最小二乘原则,计算出ai,从差分相位φdiff中减去拟合的轨道残余相位;
步骤四:CTInSAR时序分析
当SLC数据集中数据量较少时即m<20,依据相干目标频谱变化小这一特性进行相干目标的选择;点目标的实际几何尺寸往往小于雷达分辨单元,与分布目标具有不同的信号特征,不具备斑点噪声的特性,相位稳定,因而将其视为相干目标;由于点目标稳定的后向散射特性,在获取原始信号过程中各个子视图像上表现出相似的后向散射特征,不同距离向和方位向的子带宽内的散射强度大致相同,因而将SLC图像作反变换,生成复数个子视图,将子视图作谱相关分析,根据其相关性的大小来识别点目标;设各子视图内的某像元的信号xi(n)=aig(n),1≤i≤L为视数,g(n)为该像元信号,按照能量比标准衡量L个子视信号的整体相似性为:
当a1=a2=…=aL时,E=1;能量比越大,整体相似性越好,每幅SLC中频谱特征统计值E高于某一给定阈值的像元都被识别为点目标候选点,在候选点中将雷达目标的散射强度值高于某一阈值的点作为相干点;
利用相干点将去除了轨道残余相位的差分干涉纹图集的相位进行提取,得到差分干涉相位的点集;相位组成如式(4)所示,其中形变相位δφmov包含线性形变δφlinear和非线性形变δφnonlinear两部分;差分干涉相位表达为
式(7)中,λ为雷达波长,v为线性形变速率,T为生成差分干涉纹图的图像对时间基线,r为雷达到目标之间的斜距,b为垂直基线长,ε为高程误差,θ为雷达入射角;
由于在整个干涉纹图集中每个相干点的线性速率和DEM误差均为常量,是模型中的线性分量,利用以下二维模型能实现形变速率和高程改正的估计:
Ti为第i个干涉像对的时间基线;对于相干目标,由于差分干涉纹图中存在不同的相位补偿,单个像元的相位值在处理过程中并没有实际意义;在确定相位解缠参考点后连接相邻两个相干目标(xm,ym)和(xn,yn),通过式(9)求时序模型函数的γmodel最大值得到式(8)中v和ε的估算结果vestimated和εestimated;
式(9)中N为干涉纹图集中包含的干涉纹图个数;
去除线性分量,模型残余相位用下式表示:
式中:δφdiff为初始差分相位,δφlinearestmodel为线性模型将线性形变和DEM误差合成得到的相位;模型残余相位中包含大气相位影响δφatmos和非线性移动δφnonliear(x,y,Ti)以及失相干因素造成的噪声δφnoise(x,y,Ti);
利用1km×1km二维移动平均窗口进行空间低通滤波,滤波完成后,残余相位的噪声δφnoise(x,y,Ti)被滤除,最后,利用时域高通滤波估计大气贡献δφatm,区分出非线性形变相位δφnonliear(x,y,Ti);
利用前面估计出的线性分量与非线性分量获得完整的形变演变量;
依据上述相干目标时序分析方法,通过相干目标时序分析模型迭代补偿DEM误差、通过滤波减小大气相位以及噪声的影响,获取到各相干目标点的地面沉降量;
步骤五:长条带相邻分块检验与校正
对长条带中的每个frame重复步骤二至步骤四,得到条带下各分块地面沉降信息;长条带的每个frame获取的形变信息图的参考基准位置不同,需要对连续frame的参考基准进行整体偏差补偿;整体偏差估算依赖于重叠部分相干目标形变参数的统计分析,依据同名点按照式(12)统计连续frame之间形变速率的整体偏差,参与统计的相干目标均为满足相关模型的点;在此过程中选择连续frame中的一个作为基准,其余均相对于该frame进行基准偏差补偿,一般情况下,在每个frame时序分析过程中,根据时间基线的变化关系,将在不同时间间隔下稳定区域内的相干目标作为参考基准,以此降低frame之间形变参数的整体偏差;
以两块连续frame之间的检验、校正,选择其中一个frame作为基准,依据frame之间重叠区域同名点的统计信息进行基准偏差的整体补偿
式(12)和(13)中:Δρoff为frame之间基准偏差;ρmi为基准frame重叠区中相干目标i的形变量;ρsi为需调整frame重叠区中相干目标i的形变量;N表示两个frame包含的同名点个数;为调整的frame中所有相干目标改正后的形变向量;
对于重叠区域内的各相干点目标形变量ρi的取值采用距离加权平均方法求解:
式中Pmi和Psi分别为基准frame和调整frame第i个相干目标的权,ρmi和分别为重叠区内第i个相干目标在基准frame和调整frame中的形变量,D为重叠区长度,Dmi和Dsi分别为基准frame和调整frame中相干目标距离重叠区内各自影像方位向边缘的距离;
步骤六:多轨道坐标系统统一
通过对不同轨道下的多个长条带InSAR地面沉降信息集成进而获取大覆盖范围的监测成果,首先需要对不同的轨道数据进行坐标系统一,具有相同的空间坐标基准;
首先需要获取研究区的DEM数据,根据外部DEM将雷达图像与监测成果进行正射校正,即地理编码,使得不同轨道下的雷达图像位于相同的地面坐标系下,消除因地形起伏以及成像几何关系引起的畸变因素的影响;高程精度优于10m的DEM数据保证了平坦地区配准精度优于一个像元,依据卫星参数将DEM模拟至SAR坐标系,建立初始坐标变换表,再将DEM模拟的影像与SAR图像进行精密配准改进坐标变换表,通过反向变换,将雷达坐标系下各轨道地面形变信息地理编码至DEM所在的地面坐标系内;
步骤七:多轨道信息检验与参考基准统一
位于相同地面坐标系的不同轨道形变信息集成过程中,受制于雷达成像入射角的影响,不能像长条带连续frame重叠区内的相干目标在相邻frame中基本一致,相邻轨道下相干目标是完全不同的,选择其中一条轨道为主轨道,相邻另一条为辅轨道,首先通过反距离加权插值方法将辅轨道上相干目标分析结果进行空间域插值为面,设辅轨道重叠区内共有n个相干目标点,其插值后平面各点坐标(x,y)的沉降速率值通过下式计算得出:
式中ρi为第i个相干目标点的沉降量,di为点(x,y)到第i个相干目标点的距离;再利用主轨道上重叠区内相干点目标位置提取其在辅轨道上插值面上的形变值ρs;根据相同相干点目标在主轨道的沉降值ρm和在辅轨道的沉降值ρs进行检验与校正,检验与校正方法同长条带相邻分块的检验与校正。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410353927.3A CN104122553B (zh) | 2014-07-23 | 2014-07-23 | 一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201410353927.3A CN104122553B (zh) | 2014-07-23 | 2014-07-23 | 一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN104122553A true CN104122553A (zh) | 2014-10-29 |
CN104122553B CN104122553B (zh) | 2017-01-25 |
Family
ID=51768042
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201410353927.3A Expired - Fee Related CN104122553B (zh) | 2014-07-23 | 2014-07-23 | 一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN104122553B (zh) |
Cited By (17)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105204079A (zh) * | 2015-08-31 | 2015-12-30 | 中国科学院测量与地球物理研究所 | 一种利用TanDEM-X双站InSAR提取地震滑坡体积的方法 |
CN105467390A (zh) * | 2016-01-12 | 2016-04-06 | 重庆大学 | 一种基于地基InSAR的桥梁形变近距离监测方法 |
CN106093938A (zh) * | 2016-05-17 | 2016-11-09 | 长安大学 | 一种基于人工角反射器偏移量的矿区形变监测方法 |
CN108050964A (zh) * | 2018-01-30 | 2018-05-18 | 长沙深之瞳信息科技有限公司 | 一种基于微波干涉的二维面形变监测方法及系统 |
CN108088358A (zh) * | 2017-12-18 | 2018-05-29 | 电子科技大学 | 一种基于多基线雷达轨道形变检测方法 |
CN108303735A (zh) * | 2018-01-30 | 2018-07-20 | 单新建 | 基于最优参数设置的子带干涉测量的地震形变获取方法 |
CN108802727A (zh) * | 2018-04-13 | 2018-11-13 | 长沙理工大学 | 一种顾及流变参数的时序InSAR公路形变监测模型及解算方法 |
CN109059849A (zh) * | 2018-09-28 | 2018-12-21 | 中国科学院测量与地球物理研究所 | 一种基于遥感中InSAR技术的地面沉降预测方法 |
CN109633648A (zh) * | 2019-01-22 | 2019-04-16 | 北京航空航天大学 | 一种基于似然估计的多基线相位估计装置及方法 |
CN109696152A (zh) * | 2019-02-13 | 2019-04-30 | 太原理工大学 | 一种低相干性区域地面沉降量估算方法 |
CN109886910A (zh) * | 2019-02-25 | 2019-06-14 | 榆林市国土资源局 | 外部数字高程模型dem修正方法及装置 |
CN110018476A (zh) * | 2019-05-20 | 2019-07-16 | 月明星(北京)科技有限公司 | 时间差分基线集时序干涉sar处理方法 |
CN110058236A (zh) * | 2019-05-21 | 2019-07-26 | 中南大学 | 一种面向三维地表形变估计的InSAR和GNSS定权方法 |
CN110132229A (zh) * | 2019-05-10 | 2019-08-16 | 西南交通大学 | 一种铁路轨道控制网三角高程测量与数据处理的方法 |
CN111273293A (zh) * | 2020-03-03 | 2020-06-12 | 中南大学 | 一种顾及地形起伏的InSAR残余运动误差估计方法及装置 |
CN112485790A (zh) * | 2020-11-23 | 2021-03-12 | 湖南中大检测技术集团有限公司 | 基于k波段雷达的轨道非接触式变形高精度测量方法 |
CN117630932A (zh) * | 2023-11-29 | 2024-03-01 | 中国科学院空天信息创新研究院 | 一种基于分组稀疏的机载阵列干涉sar层析成像方法 |
Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20090237297A1 (en) * | 2008-02-06 | 2009-09-24 | Halliburton Energy Services, Inc. | Geodesy Via GPS and INSAR Integration |
WO2012072976A2 (en) * | 2010-11-30 | 2012-06-07 | Halliburton Energy Services Inc. | Evaluating surface data |
CN103645476A (zh) * | 2013-12-18 | 2014-03-19 | 中国国土资源航空物探遥感中心 | 一种合成孔径雷达差分干涉图序列的时空同质滤波方法 |
CN103675790A (zh) * | 2013-12-23 | 2014-03-26 | 中国国土资源航空物探遥感中心 | 一种基于高精度DEM提高InSAR技术监测地表形变精度的方法 |
-
2014
- 2014-07-23 CN CN201410353927.3A patent/CN104122553B/zh not_active Expired - Fee Related
Patent Citations (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
US20090237297A1 (en) * | 2008-02-06 | 2009-09-24 | Halliburton Energy Services, Inc. | Geodesy Via GPS and INSAR Integration |
WO2012072976A2 (en) * | 2010-11-30 | 2012-06-07 | Halliburton Energy Services Inc. | Evaluating surface data |
CN103645476A (zh) * | 2013-12-18 | 2014-03-19 | 中国国土资源航空物探遥感中心 | 一种合成孔径雷达差分干涉图序列的时空同质滤波方法 |
CN103675790A (zh) * | 2013-12-23 | 2014-03-26 | 中国国土资源航空物探遥感中心 | 一种基于高精度DEM提高InSAR技术监测地表形变精度的方法 |
Non-Patent Citations (1)
Title |
---|
葛大庆: ""区域性地面沉降InSAR监测关键技术研究"", 《中国博士学位论文全文数据库 基础科学辑》 * |
Cited By (24)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN105204079A (zh) * | 2015-08-31 | 2015-12-30 | 中国科学院测量与地球物理研究所 | 一种利用TanDEM-X双站InSAR提取地震滑坡体积的方法 |
CN105467390A (zh) * | 2016-01-12 | 2016-04-06 | 重庆大学 | 一种基于地基InSAR的桥梁形变近距离监测方法 |
CN106093938A (zh) * | 2016-05-17 | 2016-11-09 | 长安大学 | 一种基于人工角反射器偏移量的矿区形变监测方法 |
CN108088358A (zh) * | 2017-12-18 | 2018-05-29 | 电子科技大学 | 一种基于多基线雷达轨道形变检测方法 |
CN108088358B (zh) * | 2017-12-18 | 2019-08-20 | 电子科技大学 | 一种基于多基线雷达轨道形变检测方法 |
CN108303735A (zh) * | 2018-01-30 | 2018-07-20 | 单新建 | 基于最优参数设置的子带干涉测量的地震形变获取方法 |
CN108050964B (zh) * | 2018-01-30 | 2023-04-18 | 长沙深之瞳信息科技有限公司 | 一种基于微波干涉的二维面形变监测方法及系统 |
CN108050964A (zh) * | 2018-01-30 | 2018-05-18 | 长沙深之瞳信息科技有限公司 | 一种基于微波干涉的二维面形变监测方法及系统 |
CN108802727A (zh) * | 2018-04-13 | 2018-11-13 | 长沙理工大学 | 一种顾及流变参数的时序InSAR公路形变监测模型及解算方法 |
CN109059849A (zh) * | 2018-09-28 | 2018-12-21 | 中国科学院测量与地球物理研究所 | 一种基于遥感中InSAR技术的地面沉降预测方法 |
CN109633648A (zh) * | 2019-01-22 | 2019-04-16 | 北京航空航天大学 | 一种基于似然估计的多基线相位估计装置及方法 |
CN109633648B (zh) * | 2019-01-22 | 2022-04-29 | 北京航空航天大学 | 一种基于似然估计的多基线相位估计装置及方法 |
CN109696152A (zh) * | 2019-02-13 | 2019-04-30 | 太原理工大学 | 一种低相干性区域地面沉降量估算方法 |
CN109886910A (zh) * | 2019-02-25 | 2019-06-14 | 榆林市国土资源局 | 外部数字高程模型dem修正方法及装置 |
CN109886910B (zh) * | 2019-02-25 | 2021-07-13 | 榆林市国土资源局 | 外部数字高程模型dem修正方法及装置 |
CN110132229A (zh) * | 2019-05-10 | 2019-08-16 | 西南交通大学 | 一种铁路轨道控制网三角高程测量与数据处理的方法 |
CN110018476A (zh) * | 2019-05-20 | 2019-07-16 | 月明星(北京)科技有限公司 | 时间差分基线集时序干涉sar处理方法 |
CN110058236A (zh) * | 2019-05-21 | 2019-07-26 | 中南大学 | 一种面向三维地表形变估计的InSAR和GNSS定权方法 |
CN110058236B (zh) * | 2019-05-21 | 2023-04-07 | 中南大学 | 一种面向三维地表形变估计的InSAR和GNSS定权方法 |
CN111273293B (zh) * | 2020-03-03 | 2021-11-23 | 中南大学 | 一种顾及地形起伏的InSAR残余运动误差估计方法及装置 |
CN111273293A (zh) * | 2020-03-03 | 2020-06-12 | 中南大学 | 一种顾及地形起伏的InSAR残余运动误差估计方法及装置 |
CN112485790A (zh) * | 2020-11-23 | 2021-03-12 | 湖南中大检测技术集团有限公司 | 基于k波段雷达的轨道非接触式变形高精度测量方法 |
CN112485790B (zh) * | 2020-11-23 | 2023-11-24 | 中大智能科技股份有限公司 | 基于k波段雷达的轨道非接触式变形高精度测量方法 |
CN117630932A (zh) * | 2023-11-29 | 2024-03-01 | 中国科学院空天信息创新研究院 | 一种基于分组稀疏的机载阵列干涉sar层析成像方法 |
Also Published As
Publication number | Publication date |
---|---|
CN104122553B (zh) | 2017-01-25 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN104122553B (zh) | 一种集成多轨道、长条带CTInSAR的区域性地面沉降监测方法 | |
CN104111456B (zh) | 一种高速铁路沿线地表形变高分辨率InSAR监测方法 | |
CN104111457B (zh) | 一种升降轨PSInSAR地面沉降监测结果的互检验与时序融合方法 | |
Ventura et al. | Tracking and evolution of complex active landslides by multi-temporal airborne LiDAR data: The Montaguto landslide (Southern Italy) | |
Ding et al. | Ground subsidence monitoring in Hong Kong with satellite SAR interferometry | |
CN104123464A (zh) | 一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法 | |
CN109388887B (zh) | 一种地面沉降影响因素定量分析方法及系统 | |
Notti et al. | Models to predict Persistent Scatterers data distribution and their capacity to register movement along the slope | |
Favalli et al. | LIDAR strip adjustment: Application to volcanic areas | |
CN109031301A (zh) | 基于PSInSAR技术的山区地形形变提取方法 | |
CN102176016A (zh) | 一种大斜视滑动聚束sar成像处理方法 | |
Liao et al. | Ionospheric correction of InSAR data for accurate ice velocity measurement at polar regions | |
CN114791273B (zh) | 一种针对滑坡的InSAR形变监测结果解释方法 | |
CN113281749B (zh) | 一种顾及同质性的时序InSAR高相干点选取方法 | |
Suganthi et al. | Microwave D-InSAR technique for assessment of land subsidence in Kolkata city, India | |
Rigo et al. | Monitoring of Guadalentín valley (southern Spain) through a fast SAR Interferometry method | |
Ni et al. | Co-registration of two DEMs: impacts on forest height estimation from SRTM and NED at mountainous areas | |
González et al. | Relative height accuracy estimation method for InSAR-based DEMs | |
CN108872989A (zh) | 一种基于最大周期图的PS-InSAR精确搜索方法 | |
Haque et al. | Time series analysis of subsidence in Dhaka City, Bangladesh using Insar | |
Cuong et al. | Ground subsidence monitoring in Vietnam by multi-temporal InSAR technique | |
Zhu et al. | Monitoring of Surface Subsidence of the Mining Area Based on SBAS. | |
Zhu et al. | Monitoring of surface subsidence of the mining area based on SBAS | |
Takami | Monitoring Artificial Islands Subsidence in North Jakarta Using Persistent and Distributed Scatterers InSAR Analysis | |
Schumann et al. | Hydraulic and event knowledge to reduce the positional uncertainty in SAR flood images for improved flood model calibration and development |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
C06 | 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 | ||
CF01 | Termination of patent right due to non-payment of annual fee |
Granted publication date: 20170125 Termination date: 20180723 |