SAR干涉测量主要内容基本概念SAR干涉测量原理SAR干涉测量基本算法SAR差分干涉测量理论SAR干涉测量应用1本章主要内容基本概念SAR干涉测量原理SAR干涉测量基本算法SAR差分干涉测量理论SAR干涉测量应用相关词汇Synthetic Aperture Radar InterferometryInterferometricSARInSARIFSARInterferogram(干涉纹图)Phase Unwrapping(相位解缠)Permanent Scatterers(PS)2What is InSAR?InterferometricSAR(InSAR) exploits the phase differencesof at least two complex-valued SAR images acquired from different orbit position and/ or at different times. The information derived from these interferometricdata sets can be used to measure several geophysical quantities, such as topography, deformations(volcanoes, earthquakes, ice fields), glacier flows, ocean currents, vegetation properties, etc.-R.Bamler,DLRWhy InSAR?如何利用SAR数据提取目标的高程信息？B传统的方法HR1R2立体测量P3SAR立体测量提取DEM的精度精度（m）卫星地形类型同侧异侧SIR-A起伏大100-SIR-B中等变化起伏大25-6036ERS-1/2中等变化起伏大202045-JERS-1起伏大75-ALMAZ-1起伏大30~50-Radarsat-1中等变化较平坦8~1020起伏大15~204025~30-但雷达图像立体测图技术仍存在一定的问题，主要是侧视雷达图像的某些特性的不定性所造成的，如在山地，透视收缩与阴影都会造成图像形变，影响其测绘精度。WHO DOES THE FIRST??1946年，Ryle和Vonberg构造了类似Michelson-Morley干涉仪产生的无线电波，并能对一些新的宇宙电波进行定位。1969年Rogers和Ingalls将无线电波干涉测量技术用于观测金星表面。美国军方首次将机载SAR干涉测量技术应用InSAR研究的蓬勃发展于地形测绘，利用相位差图像获取高程信息，并于1971年申请了相关专利。1972年Zisk采用同样的方法来测量月球的地形。1974年Graham利用机载合成孔径雷达数据获取了能满足1:25万地形图要求的高程数据，开创了InSAR技术在对地观测中获取三维信息的先河。4InSAR的发展1969年Rogers 和Ingalls首次将InSAR技术用于观测金星表面。1974年Graham利用机载合成孔径雷达数据获取了能满足1：25万地形图要求的高程数据，开创了InSAR技术在对地观测中获取三维信息的先河。1986年Zebker和Goldstein首次将Seasat数据应用于地形测绘。1989年Gabriel等将Seasat数据应用于微小高程变化的制图工作中。1991年欧洲空间局（ESA）ERS-1卫星发射，InSAR技术应用及其数据处理引起人们极大的关注。1994年航天飞机成像雷达SIR-C/X-SAR任务将干涉测量向前迈了一大步，首次实现了多频多极化干涉测量。1995年开始了InSAR研究的一个新分支——极化InSAR的研究。2000年2月11日美国宇航局（NASA）和国家影像与测绘局（NIMA）联合进行了为期11天的航天飞机雷达地形测绘任务（SRTM）。2002年3月发射了ENVISAT-1，ASAR具有同极化和交叉极化两种极化模式以及可以多角度，多模式成像。2006年1月日本计划发射的ALOS/PALSAR具有获取多极化、极化、极化干涉和多种模式数据的能力。2007年6月德国TerraSAR-X发射，COSMO-SkyMed2于同年12月发射。2007底12月加拿大在RADARSAT-1基础上发射的RADARSAT-2也具有全极化和干涉测量的功能。19691974198619891991199419952000200220062007本章主要内容基本概念SAR干涉测量原理SAR干涉测量基本算法SAR差分干涉测量理论PSINSAR技术SAR干涉测量应用5HOW DOES IT WORK?A2B天线A1Return could befrom anywhere 天线1on this circle天线2Return comes fromintersectionSARInSARSatellite orbit 2PerpendicularBaselineSatellite orbit 1InterferometerBaselineAzimuthAntenStrip -mapna footprintGround rangeSlant range6Interferometriczsystem: monostatic(ERS) Su(t)exp(−jω0t)S12sPr)∝u(t−τ)exp(−jωr1(10(t−τ1))exp(jϕs(P))+n1(tτ=2rPr)11()/cs(Pr)∝u(t−τr22)exp(−jω0(t−τ2))exp(jϕs(P))+n2(t)θτ(Pr2=2r2)/cHPy (ground range)rInterferometriczsystem: bistatic(SRTM 2000)SS1u(t)exp(−jω20t)s(Pr)∝u(t−τt−τr11)exp(−jω0(1))exp(jϕs(P))+n1(τrrPrt)1=(1(P)+r1())/cs(Pr)∝u(t−τr22)exp(−jω0(t−τ2))exp(jϕs(P))+n2(Hθτrrt)2=(r1(P)+r2(P))/cPy (ground range)r7IMPLEMENTATION ADVANTAGES DISADVANTAGES Simultaneous Baseline • Known baseline • Difficult to get adequate baseline in space • No temporal decorrelation • High data rate from two radars • Typically better performance • Typically higher cost Repeat Track • Lower data rate from one • Temporal decorrelation radar • Lower cost • Baseline not well known and may be changing • Depending on orbit, any baseline can be realized The SAR interferogram is generated by:Cross -multiplying pixel by pixel the first SARimage times the second one complex conjugated.Thus, the interferogram amplitude is the amplitude of the first image times that of the second one.Whereas its phase (called interferometric phase) is the phase difference between the two images.由于入射角的差异使得两幅SAR图像不是完全重合，对它们进行配准处理后，配准后的图像对进行复共轭相乘就得到了复干涉纹图（interferogram）8InteferometricphasesInSAR的几何原理s1(R)=u1(R)exp(iφ(R))s2(R+∆R)=u2(R+∆R)exp(iφ(R+∆R))zA2AB1αϕ2π1=2λR+arg{u1}θϕ2=22πλ(R+∆R)+arg{u2}R+∆RhRs*1(R)s2(R+∆R)=s**⎛⎜4π⎞1s2expi(ϕ1−ϕ2)=s1s2exp⎝−iλ∆R⎟⎠Z(y)φ=−4πλ∆R+2πNN=0,±1,±2,LyInSAR成像几何示意图（如果地表无形变）9=(R+∆R)2−R2sin(θ−α)−B2B||=Bsin(θ−α)2RBzA2B⊥=Bcos(θ−α)AB1αz=H−RcosθθR+∆R∆R≈Bsin(θ−α)+B22RhRθ=α−arcsin⎡⎢λφ⎤⎣4πB⎥⎦Z(y)如果知道天线位置参数和雷达y成像系统参数等，就可以从相位中计算出地表的高程值。InSAR成像几何示意图（如果地表无形变）SAR interferometric phase: geometric SAR interferometric phase: geometric contributioncontributionPhase≈4πrBλ∆nBr2re1gnar tnalSGround rangeAzimuthAzimuth10斜距变化和高度变化两个方面研究干涉相位的变化情况。4ABA2π1αφ′=−θR−α)θλBsin(θ0+∆0∆φR=φ′−φ=−4πλ[Bsin(θ0+∆θR−α)−Bsin(θ0−α)]=−4πλBcos(θ0−α)∆θRΔRPP'R∆θR≈Rsin∆θR=∆R/tanθ0高度不变，斜距发生变化∆φ4πBcos(θ0−α)∆R4πB⊥∆R无高程变化的平坦地表也会产生线性R=−变化的干涉相位，称之为平地效应。λRtanθ=−0λRtanθ0由于平地效应的影响，常常会造成干涉条纹过密，影响后面的相位解缠11Interferogramflattening Mt. Vesuvius, baseline 250 m.12Interferogramflattening BA2A4π1αφ′=−sin(θ0+∆θz−α)θλB0∆φ4πz=φ′−φ=−[Bsin(θ0+∆θz−α)−Bsin(θ0−α)]RRλP'=−4πλBcos(θ0−α)∆θzΔzPR∆θz≈Rsin∆θz=∆z/sinθ0高度变化，斜距不变∆φ4πBcos(θ0−α)∆z4πB⊥∆zz=−λRsinθ=−0λRsinθ0比较∆φπBcos(θ0−α)∆RR=−4λRtanθ=−4πB⊥∆R0λRtanθ013重要的指标！干涉相位的高度敏感度∆φz∆z=−4πB⊥λRsinθ0高度模糊数——引起一个2π相位变化所对应的高度变化，以此来表征干涉测量对高度变化的敏感程度∆zλRsinθ02π=−2B⊥Elevation An elevation ∆qover the earth surface, gives, to θthe flattened interferogram, the differential phase: ∆φ=−∆qBn4πsinθ⋅R⋅0λThe altitude of ambiguity(∆φ=2π) is then:∆q×Bλn=sinθR0O2P×97∆q’≈97ye.g.: Bn=120 Æ∆q=80 mP14Interferogramflattening eganr tanlSAzimuthMt. Vesuvius, baseline 50 m.The role of the normal baselineThe role of the normal baselineBn=50Bn=250egnar tnalSAzimuthThe The normal (or perpendicular) baselinenormal (or perpendicular) baselineis a keyis a keyparameter in SAR interferometryparameter in SAR interferometry15Interferogramand DEM∆zRsinθ02π=−λ2B⊥16本章主要内容基本概念SAR干涉测量原理SAR干涉测量基本算法SAR差分干涉测量理论SAR干涉测量应用主图像SLC辅图像SLC算法高精度配准流干涉纹图程去平地效应简图相位解缠HOT地理编码DEM17FIRST：数据的选择数据选择的原则：基线距与时间相干（1）SLC数据（实部＋虚部）（2）提取DEM的干涉图像对的基线距在200-500米之间。（3）ERS-1/2 tandem模式时间去相干较小。临界基线距BWRtan(θ−ϑ)⊥C=λ2ERS-1/2 SAR系统临界基线距大约为1100m 18SECOND：高精度配准XβRkrRMasterSlavekxX通常，图像的配准误差必须在1/8个像元以下才对干涉纹图的质量没有明显的影响SECOND：高精度配准距离向配⎧sin(πµr)准误差γcoreg,r=⎪⎨sinc(µr)=πµ,0≤µr≤1⎪r⎩0,µr>1通常，图像的配准误差必须在1/8个像元以下才对干涉纹图的质量没有明显的影响19距离向21fDC2fDC1方位向在重复轨道的两次成像过程中，主辅图像在方位向因为斜视角不同，导致回波信号多普勒中心频率之间存在偏移，两幅图像方位谱也不匹配。在距离向，因基线距而导致影像所反映的地距频率范围也不相同，这个偏移反映到距离频谱上，使得主、辅图像的频谱也存在非公共距离频率偏移。（a）主像滤波前；（b）辅图像滤波前；方位滤波前后方位频谱对比（c）主图像滤波后；（d）辅图像滤波后THIRD：生成干涉纹图配准（interferogram）后的图像对进行复共轭相乘就得到了复干涉纹图s*⎛4π1(R)s2(R+∆R)=s*=s*⎞1s2expi(ϕ1−ϕ2)1s2exp⎜⎝−iλ∆R⎟⎠虽然对单幅图像来说相位是随机的，但干涉纹图的相位是确定的，它取决于信号的路径差ΔR。20Complex multilookA “cleaned”interferogram is achieved by averaging in areas of uniform phase.SNR improves ∝the number of looks. Usually the averaging window is adaptive. Noise sources in interferometry SAR acquisitionFocusing SceneHIRF1HFOC1reflectivityMulti-lookv(xFocusing filtering0,r0)SAR acquisitionHIRF2HFOC2Scene decorrelation noiseSystem decorrelation noiseInterferogram(temporal, volume, geometry)(thermal, quantization, alias)Phase noise sourcesProcessor decorrelation noise(misalignment, phase error ...)For evaluating noise performances, the interferometric phase difference (that is usually topographic-dependent) is assumed to be zero.21重要的指标！干涉相干γ=s1s*2s**1s1s2s2量度干涉条纹的质量提供散射体的重要信息，可用于地物分类∑∑MNs∗jφ(m,n)1(m,n)⋅s2(m,n)⋅e−MLE相干估算器γˆ=m=0n=0∑∑MNMNs1(m,n)2n)2m=0n=0m∑∑s2(m,=0n=0去相干源主要可分为六大类：①基线或几何去相干，由视角差异造成的；②多普勒质心去相干，主要是两幅图像的多普勒质心存在差异；③体散射去相干；④热噪声去相干，主要受系统特征影响，包括增益和天线特征；⑤时间去相干，由获取图像对期间地表散射特性发生变化造成的；⑥数据处理去相干，主要取决于数据处理过程所采用的算法。γtot=γgeom×γDC×γvol×γthermal×γtemporal×γprocessing22Correlation mapsInterferogramsAll images from Space Shuttle (SIR-C) span Apr-Oct From: Rosen et al., 199623干涉纹图干涉相干图Coherence mapsAmplitudePhaseCoherence24Atmospheric artifactsAtmospheric inhomogeneities, due to variations in Pressure, Temperature, Humidity affects light speed velocity. The delay variation in the repeat intervalresults in a “phase screen ”, hence a noise.“Atmospheric”noise is fractal (f-a) power spectrum (a~2/3), hence it is correlated in space (~ hundred of meters). Coherence maps cannot measure this noise.Atmospheric artifacts can be up to two fringes. This is converted in elevation error, depending on the baseline. The error cannot be estimated or recovered.25Atmospheric artifactsAtmospheric artifacts can be detected in a multi-baseline environment, since atmosphere is incorrelatedin time.EtnaParis (+ reflectivity)FOURTH：去平地效应无高程变化的平坦地表也会产生线性变化的干涉相位，称之为平地效应。A4π2φ′=−ABλBsin(θ0+∆θR−α)1αθ0∆φR=φ′−φ=−4πλ[Bsin(θ0+∆θR−α)−Bsin(θ0−α)]π∆θ=−4λBcos(θ0−α)∆θRrΔRPP'πBcos(θ0−α)∆R4πB⊥∆R平地效应示意图∆φR=−4λRtanθ=−0λRtanθ026目前常见的平地效应消除方法主要有：基于轨道参数和成像区域中心点的大地经纬度计算平地效应；基于已有的DEM数据计算平地效应。通过测量距离向和方位向的占优势的干涉条纹频率来计算的平地相位，然后进行平地效应的消除。以苏州地区为例27FIFTH：相位解缠一切将相位由主值或相位差值恢复为真实值的过程统称为相位解缠(phase unwrapping)。除了合成孔径雷达干涉测量要进行相位解缠，在信号与图像处理（如同态解卷积）、合成孔径声纳、自适应光学、声音成像、光学与微波干涉、补偿式成像、核磁共振（MR）、地震处理等方面相位解缠都有重要的应用。Phase unwrapping28Phasez8πPhase⇒⇒6π4π2π2π00yyyACTUAL ELEVATION PROFILEWRAPPED PHASEUNWRAPPED PHASE解缠前后相位图29φi,j=ψi,j+2πki,jki,j为整数,−π≤ψi,j<πi∈[0,M−1],j∈[0,N−1]目前所有的相位解缠可分为两个步骤：①基于缠绕相位计算解缠相位的相位梯度估算值；②积分。根据所采用的积分方法，相位解缠方法主要分为两大类：路径跟踪法和最小范数法。现有的相位解缠算法都是基于这样一个假设：有可能推导出缠绕相位的离散偏导数，即邻近像元的相位差，并且这些相位差的绝对值都小于π。通过这些离散的偏导数，可以重建解缠相位。0.5 0.75 0.0 0.25 0.5 0.75 0.0 0.25(a)-0.5 -0.25 0.0 0.25 0.5 0.75 1.0 1.25(b)-0.5 -0.25 1.0 1.25 1.5 1.75 2.0 2.25(c)一维相位解缠的实例：(a)缠绕相位（每个值乘以2π）；(b) 满足解缠假设的解缠结果；(c)另一种可能的相位解缠结果，但箭头所指处的相位差大于2π。30−π≤∆[φ]<π二维相位场具有一致性xi,j−π≤∆y[φi,j]<π相位场是连续的（consistent）或是“无旋”的：①相位值处处有定义。如果某个像元点的干涉强度为0，则该像元就没有任何相位信息。“一致性”要求干涉纹图的强度必须不为0；②相邻像元的真实相位差在[−π,π)之间。但真实干涉纹图很多像元的相位信息都不满足一致性条件（有旋场），主要由以下二类因素造成的：①低相干；②地形因素，如顶底位移、伪相位信号以及雷达阴影等。∫∇φ(rr)⋅drr+∫rrrCC[∇×A(r)]⋅dr无旋分量有旋分量m0.1Δ0.24residuenΔ1Δ3-0.1Δ2-0.2∆1=W[ψm,n+1−ψm,n]∆2=W[ψm+1,n+1−ψm,n+1]∆3=W[ψm+1,n−ψm+1,n+1]∆4=W[ψm,n−ψm+1,n]4q=∑∆⎧=0,连续节点i⎨i=1⎩≠0，非连续节点31m0.2Δ40.3nΔ1Δ3-0.2Δ2-0.4有否残差点？•路径跟踪的相位解缠方法主要包括枝切法（切割线法）、掩模割线法（Mask-cut法）、像元排序或像元扩散法、最小跨越树法、加权的最小不连续法和区域生长法等。输入：缠绕相位（相位质量信息）步骤1：相位连续性/不连续性检测：识别残差点，生成枝切线步骤2：计算/建立相位质量图步骤3：相位积分：在枝切线周围或在质量图的指导下如果缠绕相位包含多个独立的区域，选择起始的图像块重复（Repeat）选择起始积分点确立将要进行相位解缠的候选数据列重复（Repeat）更新候选数据列对候选数据列进行相位解缠（依据对应的数据质量）直到（Until）无候选数据直到（Until）所有的缠绕数据块被处理完毕32+_+_+_残差点之间连接枝切线的示意图：(a)正残差点与负残差点；(b)连接残差点的枝切线，用深灰色的像元表示；(c)枝切线不包括左上方像元。+++___Goldstein枝切法的基本操作过程：(a)正残差点与负残差点；(c)搜索窗扩大为5×5，枝切线用深灰色像元表示。(b)3×3的搜索窗；+++_+_+_+表示在有第三个残差点的情况下如何连接枝切线的示意图：(a)在图像的右下方有一新的正残差点；(b)连接新残差点的枝切线；(c)将残差点与边界相连。利用枝切法设置的枝切线33干涉纹图与利用Goldstein算法得到的解缠结果A real example33221100-1-1-2-2-3-3Original phaseIrrotationalphase3502401300-120-210-30Rotational phaseUnwrapped irr. phase34Unwrapping 1Identification of disconnected zones35Unwrapping 2Interpolation36最小范数法Min{J}=Min⎧M−2N−1M−1N−2⎨⎩∑∑φpp⎫i+1,j−φi,j−∆xi,j+∑∑φi,j+1−φi,j−∆yi,j⎬i=0j=0i=0j=0⎭(φi+1,j−φ)(−φxi,j−∆xi,jUi,j+φi,j+1i,j−∆yi,j)Vi,j−(φi,j−φi−1,j−∆i−1,j)Ui−1,j+(φyi,j−φi,j−1−∆i,j−1)Vi,j−1=0⎧⎪φ−φp−2Ui=0,1,2,i,j=⎨i+1,ji,j−∆xi,jKM−2;j=0,1,2,KN−1⎪⎩0otherwisep−2V⎧⎪i=0,1,2,,j=⎨φi,j+1−φi,j−∆yii,jKM−1;j=0,1,2,KN−2⎪⎩0otherwise当p=2时，最小范数的相位解缠问题等效于求解牛曼边界的泊松方程。(φi+1,j−2φi,j+φi−1,j)+(φi,j+1−2φi,j+φi,j−1)=ρi,jρi,j=(∆xi,j−∆xi−1,j)+(∆yi,j−∆yi,j−1)⎧⎪x⎨∆−1,j=0⎪⎩∆xM−1,j=0⎧⎨∆yi,−1=0⎩∆yi,N−1=037最小二乘相位解缠算法9基于FFT的无权重最小二乘算法9基于DCT的无权重最小二乘算法9基于误差方程的无权重最小二乘算法9PCG算法9加权重的多网格算法9Lp-Norm算法（a) 513×513的干涉纹图；(b)基于FFT的无权重的最小二乘算法相位解缠结果；(c)相干图；(d) Goldstein枝切法相位解缠结果38LAST：高程转换及地理编码∆ϕ=4πλRBn∆q0sinθ&phase to height conversiongeocoding实验区1:50000地形图地理编码后的数字高程模型39新疆卡拉喀什地区三维透视图三维透视图404142SRTMOn February 11, 2000, the Shuttle Radar Topography Mission (SRTM) payload onboard Space Shuttle Endeavour launched into space. With its radars sweeping most of the Earth's surfaces, SRTM acquired enough data during its ten days of operation to obtain the most complete near-global high-resolution database of the Earth's topography. To acquire topographic (elevation) data, the SRTM payload was outfitted with two radar antennas. One antenna was located in the shuttle's payloadbay, the other on the end of a 60-meter (200-foot) mast that extended from the payload pay once the Shuttle was in space.43SRTM DEM ProductRaster size 1\"x1\" Lon& LatHeight levels 1mDatum (horizontal)WGS84Datum (vertical)WGS84 or MSL (optionally)Data format16-bit Signed Integerhorizontal accuracy (absolute)±20m 90% circular errorhorizontal accuracy (relative)±15m 90% circular errorvertical accuracy (absolute)±16m 90% vertical errorvertical accuracy (relative)±6m 90% vertical errorSalt Lake City, Utah44SAR差分干涉测量45464ππA3φ1≈−4πλB||−λ∆Rd=−4πλBsin(θ1−α1)−4λ∆RdAB'2B||B−αAαα2−4π1φ2≈2)1λB||′=−4πλB′sin(θ2RR32ZR1Pφd=φ1−B||4πB′φ2=−λ∆Rd||三轨法用“去平地”相位重新表示得φd=φ−B⊥f1B′φ4πf2=−⊥λ∆Rd如果在差，造成位移图的变形。“去平地”处理的过程中，如果所用的基线距的值并不是真实的基线距值，将会引入误差分相位对地表形变的敏感度φd∆R=−4πdλ47基准辅SAR复影像基准主SAR复影像观测SAR复影像复图像粗配准复图像粗配准复图像精配准复图像精配准辅影像重采样观测影像重采样辅影像滤波主影像滤波观测影像滤波复相干图1单视干涉纹图1单视干涉纹图2复相干图2平地效应消除1平地效应消除2相位解缠1相位解缠2差分干涉测量区域性地壳垂直形变484950Major earthquakes studied using radar interferometryzzzzzzzzzzzzzzGigiearthquake, CHinaZhangbeiParkfieldearthquake,ChinaCreep San Andres fault/ ItalyColfiorito, Umbria-Marche, GreeceGrevrenaearthquake, Elat(Aqaba)Nuweibaearthquake, Gulf of earthquake,JapanKagoshima-kenhokuseibuCaliforniaNorthridge earthquake, Eureka valley,CaliforniaIzimitManyiearthquake,TurkeyearthquakeHector earthquake,Californiaearthquake,CaliforniaLanders Kobe,Japan张北-尚义地震About 80% houses in VIII intensity region were destroyed due to their poor quality of constructions, which were mostly made of rubles or adobe walls.51NQijiaTaolizhuangShangyidispalcementcenterFriendship XiaosuangouXinhereservoirHuaianDianzichun12KMSAR image of studyERS-1/2 SAR Data of Study areamissionorbitframeacquisitionBtime time//(m)B⊥(m)intervalERS-2 SLCI12683278123/9/1997000ERS-1SLCI32356278122/9/1997-231-269-1ERS-2SLCI16190278126/5/199849-8824552Coherence map (the tandem image pair)combination of preseismicSAR interferogramdata pairformed from Before earthquake53Interferogrambefore earthquake overlaid with SAR intensity map The LOS surface displacement map overlaid with SAR intensity image54Earthquake deformation from the Landers earthquake, CaliforniaLand subsidence in Las Vegas, Nevada55Glacial movement of the PetermannGlacier, Greenland This has led to the recent discovery thatSAR interferometrycan be is much more complex than previouslyice flow within the Antarctic interiorused to calculate ice flowthought surface velocities.56Volcano inflation (simulated) of Mt. Kilauea on the big island of Hawaii 57