小基线集技术监测北京沉降.docx

land subsidence monitoring in beijing based on sbas methodsupervisor: xxxxxxxx universityjune xxxx第一章 绪论1.1问题提出及研究意义 地面沉降广泛存在于全世界的城市中，这是全球在城市化进程中出现的不可忽视的环境地质问题，这种连续且渐进的地质灾害一旦形成很难恢复，不仅给城市带来巨大危害，同时还影响到人类生命财产安危和区域的可持续发展。目前，我国已经有96个大中型城市发生了地面沉降，主要分布在长江三角洲地区，松辽、黄淮海、东南沿海和内陆河谷等平原地区，沉降总面积超48655 km2；其中全国累计沉降量超过200 mm的地区达到7.9万km2，比较严重的地区有上海、天津、台北、西安、宁波、苏州等地，地面累计沉降量达4602780 mm，地面年沉降速率达1056 mm/a，已经造成上千亿元的经济损失1。当下长江流域、黄河流域、珠三角、松辽平原和环渤海地区及东南沿海城市的地面沉降由于过量开采地下流体资源已经形成巨型漏斗地面沉降在区域上呈扩张趋势，并且对地区基础设施建设造成极坏的影像，其形成和发展速度令世界震惊2。地面沉降现象日益恶化和严重，并在部分地区伴生地裂缝、岩溶塌陷等环境地质问题，其破坏性较大，给人民生命财产与社会经济发展带来严重威胁，已经成为制约区域环境、社会经济可持续发展的主要地质灾害之一。严重的地面沉降可产生许多危害，特别是一些大中城市可能存在着较大的潜在的危险，主要表现在以下几个方面： 地面标高损失明显，地面低洼积水难排，河道泄洪能力下降，洪涝灾害加重； 地面沉降、地裂缝使建筑物基础下沉、工程设施损坏，长时期严重的地面沉降会造成地下排水管道的破坏，严重影响城市生活污水和雨后积水的排放； 地面沉降导致地下水环境恶化，影响居民生活用水。为了达到有效预防地面沉降以保证城市可持续发展和为城市发展规划提供可靠地质资料数据的目的，必须建立长期有效的地面沉降调查与监测机制并制订合理措施对沉降灾害进行预防与治理。常规的地面沉降监测一般包括高等级水准测量、地下精细观测方法（基岩标、分层标等）和gps网络技术监测等方法，这些方法虽然能够达到较高的地表形变观测精度，但这些方法往往耗费大量的人力物力，观测周期长，实测费用高，同时这些技术仅能获取离散点的形变量，而且在量测大范围的沉降灾害监测时这些技术劣势明显3。下表是几种地面沉降监测方法的比较。测量方法形变分量测量方式测量精度（mm）测量成本测量密度（每次）精密水准全站仪gpsinsar垂直水平/垂直水平/垂直视线向线线网面0.1-102/55/2010高高较高低101001010010100105107地面沉降给社会经济和人民的生命安全等方面造成巨大的危害，直接影响区域经济的可持续发展，而实时监测地面沉降可以提前做出风险预警、灾害预报与评估，避免重大的灾难事故的发生；可以为合理开采利用地下水资源、煤、油、气等矿产资源做出指导；可以作为城市的合理规划、发展、立法的依据。因此，实时监测地面的沉降对于维护人民的生命安全和经济的可持续发展具有很高的意义。 及时开展 insar 技术在地面沉降监测中的研究，可以实现地面沉降的大面积、低成本、高效率监测，尽快推动 insar 技术在我国地表形变监测领域的实际应用。因此，研究 insar 差分干涉测量的理论与实用化技术、研究利用 insar 技术监测地面沉降具有很高的研究价值和研究意义。1.2 insar技术及其发展近二十多年来，空间对地遥感观测技术发展迅速，合成孔径雷达干涉测量（interferometric synthetic aperture radar，sar）作为一种新兴的卫星微波遥感技术，尤其是在此基础上发展起来的雷达差分干涉测量技术（differential insar，d-insar）为实现地面沉降的准确、高精度、大范围测量提供了全新的手段。与传统地表形变监测方法相比，d-insar技术在理论上可以实现mm精度（高程）、同时达到高分辨率（水平方向米级）、大范围（100*100 km2以上）的探测雷达视线方向的地表形变。1.2.1 insar技术国内外现状1969年rogers等首先利用雷达复数数据提取得到地表三位信息，1972年zisk报告了用雷达干涉技术用于月球表面测绘的可行性。1974年，graham首次利用机载双天线合成孔径雷达得到地形高度信息。随后，1986年，推进器实验室的zebker和goldstein研究利用机载sar干涉获取地形数据，引起了学者的关注。1988年goldstein利用机载sar干涉技术处理了seasat sar数据，并且获得了理想的结果。li及goldstein在1990年研究了不同基线对地形制图产生的影响，因此提出了最优基线的概念。1991年，欧空局发射了c波段的ers-1，方便人们利用相隔一段时间的图像进行干涉处理。1995年，ers-2大大缩短了原有时间基线，提高了处理精度。2000年初，美国施行航天飞机测图计划（srtm），在短短的11天内，利用奋进号航天飞机成功获取了覆盖全球80%的地形数据，从而实现了全球数字高程模型的构建。在地表形变监测方面，1989年，gabriel等首次论证了利用合成孔径雷达差分干涉测量(d-insar)技术监测微小地表形变的可能性，并利用seasat l波段的sar测量了美国加里福尼亚东南部的imperial valley灌溉区的地表形变。1993年，massonnet等成功地用d-insar技术获取了1992年美国加州landers地震(m=73)的同震形变场，与测量数据及弹性形变模型获得的结果吻合很好，这一结果发表在nature上，引起了国际地学界的极大关注。d-insar早期主要是研究比较明显的地震、火山等活动的形变，随着理论和技术的发展，d-insar技术研究范围涵盖了地面沉降、山体滑坡等微小持续的地面形变。欧洲也进行了d-insar应用在地面沉降监测的试验，选取区域在波兰的煤田，该煤田由于开采，地面普遍存在下降现象。1998年，试验者通过采用50张该区域的ers雷达影像，采用d-insar的方法得到了煤田区域的形变量。在美国的加利福尼亚地区、法国的gardanne地区，分别选择2颗ers卫星的数据，利用insar技术监测了belridge、lost hills油田和gardanne附近的煤矿开采区的地面沉降，分别提取了不同时间周期内的沉降量变化值，与对应时间周期内的水准测量结果的对比表明，insar处理结果十分令人满意，证实采用insar技术进行沉降监测是可行的。1999年，wegmuller等利用1992年8月到1996年5月期间的ers数据监测意大利bologna城沉降情况，并与水准和gps数据进行分析比较，得到较一致的形变场和形变速率。为了抑制失相关和大气影响对形变信息提取的负面影响，sandwell和price率先提出了基于多干涉图堆积（stacking）平滑的方法。为了削弱大气延迟的影响，一些学者提出了利用地面气象数据或gps所获得的天顶延迟量对d-insar结果进行大气改正，然而，这些数据分辨率过低只能在一定程度上降低大气延迟的影响，无法改正大气相位，因此从d-insar结果中完全扣除大气影响较为困难。意大利人ferretti等于2000年率先提出了永久散射体干涉（persistent scatterer interferometry，ps-insar）的技术途径，不仅可以抑制失相关的影响，而且还可以分理出大气延迟与形变信息，从而提高形变观测的精度和可靠性。1999年，usai提出了利用多幅干涉图反演地表形变时间序列的方法，该方法将形变问题转化成最小二乘问题。2002年，bernardino等人bernardino和lanari等人对模型进行改进，提出了小基线集（small baseline subsets，sbas）方法，并将其用于研究大尺度的形变。sbas方法将sar数据根据垂直基线大小组合成若干个集合，有效的减弱了空间去相干影响，而且增加了时间分辨率，具有高精度的地表形变监测能力，同时该方法有多余观测，增加了形变监测的可靠性，在此基础上研究了意大利南部的campi flegrei火山和naples市区在空间低分辨率下(约100m*100m) 的时间序列形变，利用奇异值分解法的原则获得了很好的监测结果。2003年，mora等人结合了ps和sbas方法的特点，进一步提出了进行形变分析的方法，2005年，lauknes等在对挪威首都oslo的地表形变监测中，比较了sbas和ps两种不同的干涉方法，认为二者得到的相干性分布和形变模式是一致的。2005年，casu等利用小基线集(sbas)方法测量意大利naples湾和美国los angles的地表形变，并与水准测量及gps数据进行比较。2006年，casu认为参考像元的选择对标准偏差的计算值有0.05mm/km的影响。在利用 d-insar 技术进行地表形变监测应用研究方面，国内还处在探索研究阶段，但也开展了许多研究，主要的研究有： 2000 年，王超、张红等人利用差分干涉测量技术完成了1998 年张北地震同震形变场的监测，提取了地震同震形变场，并对形变特征和震源构造进行了分析；2002 年，单新建等利用西藏那曲地区玛尼乡的六景ers-1/2 数据成功地获取了玛尼地区同震形变场。2002 年，路旭和li tao等利用ers数据经过差分干涉测量处理提取了天津地区不同时间周期的地面沉降变化值，并与同时间周期的水准测量结果相比较，取得了满意的效果。王超等利用1993年至2000年ers-1 sar数据对苏州地面沉降进行了监测，通过与水准数据进行比较，所得结果与水准观测近似度达0.943，测量精度达到5mm。刘国祥等利用ers-2 数据成功获取了香港赤腊角机场在近乎一年内的非均匀沉降场，与水准测量结果相比吻合较好。张红利用雷达差分干涉测量技术成功获取了苏州地区的地面沉降场。2004年，李陶等利用ers-1 sar数据对天津市1992年至1997年地面沉降进行了研究，发现四个沉降漏斗，测得其中一个漏斗沉降速率与实测结果一致。2006年，gong等利用jers-1 sar数据对沧州地区由于地下水开采进行了研究，结果表明实测结果与所得结果近似一致。这些实例都证明了雷达干涉测量技术在我国地面沉降监测广泛应用，也预示了广阔的发展前景。1.2.2 insar影像处理软件软件名称软件开发者运行平台备注gamma瑞士gamma remote sensingwindows，linux商业化软件earth view加拿大 atlantis公司sun，windows商业化软件sarscape瑞士sarmap公司windows，linux商业化软件isar意大利sundoris荷兰代尔夫特理工windows，linux，sun源代码开放roi_pac推进器实验室（jpl）sun，linux，mac os源代码开放idiot柏林理工大学sun，linux，mac os，windows源代码开放roi斯坦福大学unix源代码开放rat柏林理工大学windows+idl源代码开放polsarpro欧洲空间局windows免费diapason法国cneswindows，linux商业化软件sarinfors中国林科院windows国内软件grace中科院遥感所windows国产软件ifsar美国erdas公司windows+erdas不能做差分1.2.3 几种典型的星载sar卫星系统国际上主要卫星sar系统及相关参数卫星sar系统发射年份轨道高度（km）波段/波长（cm）侧视角（度）重复周期（天）地面分辨率（m）影像幅宽（km）almaz-11987300s/1021655715303045ers-1/21：19922：1995790c/5.6233525100jers-11992568l/23.5384425800radarsat-11995790c/5.623652483050500envisat2002800c/5.615453525100100405alos2006700l/23.5860461010020350radarsat-22007798c/5.6104924310025500terrasar-x2007514x/3.120451111610100tandem-x2010514x/3.120451111610100cosmo-skymed（四星座）20072010620x/3.12060416110010200hj-1c2012499s/10254731310040100第二章 insar基础2.1 雷达基本原理雷达（radar）意即无线电探测与测距，最早由军方研制使用，用来探测硬目标（一般为金属点目标）及测距，这些雷达系统并不产生图像。后期发展的雷达遥感技术把地形地貌作为主要探测目标。真实孔径雷达（rar）是最早的成像雷达系统，其方位向分辨率受天线尺寸的限制，随着理论研究的深入，天线设计、信号处理及计算机软硬件的发展，合成孔径雷达（sar）开始出现，并逐步取代了真实孔径雷达。sar 系统既可以装载在飞机上，也可以装载在航天飞机或卫星上。20 世纪 90 年代机载、星载雷达遥感经历了蓬勃的发展，达到了雷达遥感发展的高潮。雷达遥感可分为主动和被动两种方式。目前多数星载雷达采用主动方式，即由遥感平台发射电磁波，然后接收辐射和散射回波信号。主动方式雷达在极短的时间内发出高频雷达波束，然后接受地面反射信号，这样便得到了地表雷达图像。雷达成像不依赖于光照和天气条件，具有区别于光学传感器的独特优势，一定条件下，微波还具有一定程度的穿透特性。真实孔径雷达（real aperture radar，rar），其雷达天线长度就是实际长度，雷达波的发射和接收均是以其自身实际有效长度的效率直接体现在显示记录中。真实孔径雷达图像可以分为距离向和方位向分辨率，统称为空间分辨率。距离向分辨率是指在脉冲发射方向上能分辨地面两个目标的最小距离；方位向分辨率是指在与辐射波束垂直方向上相邻的两束脉冲之间，能分辨出两个目标的最小距离。空间分辨率主要取决于雷达孔径的大小、脉冲持续时间、波长以及天线到地面目标的距离等。相对于真实孔径雷达（rar），为了得到方位分辨率较高的雷达图像，科学家们根据多普勒效应原理，于 20 世纪 50 年代初期由美国科学家最先提出了“合成孔径”的概念。早期主要是用来满足军事侦察对雷达高分辨率的需求，合成孔径雷达的成像方法是一种相干的主动微波成像。由于平台系统提供了精确位置以及多普勒，且本身也具有完整的获取数据的设计，所以合成孔径的处理可以顺利的完成。雷达波束覆盖时间段内，地面物体的回波幅度以及相位信息被遥感平台中的相应设备记录在回波库中。不需要物理长天线及短波长，只要对回波信号进行普勒频移处理，就可得到一个很窄天线效果。2.2 insar 的基本原理2.2.1 insar干涉测量工作方式在insar的数据采集模式中，依据接收天线位置的几何关系的不同，可以将sar干涉测量的工作模式分为：单轨双天线横向模式（xti） 、单轨双天线纵向模式(ati)、重复轨道单天线模式（rti）。下面将详细介绍每种工作方式的成像几何关系。（1）单轨双天线横向（crosstrack interferometry，xti）模式图31是sar交轨干涉测量的示意图。交轨干涉测量是利用一副天线发射雷达波，由两幅天线同时接受回波。由于两幅天线接收的回波具有一定的相干性，两幅天线与地面目标物之间的路径差为影像的相位差。相位差是由两幅天线与地面点之间的路径差造成的，路径差是与地面点的高程相联系的。因此，获取干涉测量系统的几何参数，就可以将相位信息转换为高程信息。即除了图像点的方位向和距离向的坐标外，相位差提供了高程信息，这样就能够重建地面目标点的三维特征。由于进行sar交轨干涉测量需要同时安装两幅sar天线，因此通常应用于机载系统或航天飞机任务。（2）单轨双天线纵向（along-track interferometry，ati）模式sar顺轨干涉测量技术是由goldstein和zebker在1987年提出的。实现方法也是在一般的机载sar系统上加装另一副天线，但与xti不同的是，这两副天线是顺着平台飞行轨道来安装，前后基线之间的长度通常为2-20m。假设天线之间的间距是20m，飞机的飞行速度是200m/s，则两副天线所观测到同一地物的时问间隔是0.1s，在这个时间间隔内时钟漂移和传播延迟的变化都可以忽略不计的，因此测量到的相位差反映的是地面点的高程变化信息。这项技术被称为“sar顺轨干涉测量”(ati)。安装的同一架飞机上的两幅天线是沿着飞行方向的。ati常用于水流制图、动目标检测以及定向波谱的测量。上面两种方式都属于单轨双天线干涉，这种方式的优点有：由于可以同时获取两幅sar影像，在两幅影像中地物具有相同的的后向散射系数，获取数据的相干性比较好；这种方式是一个天线发送，由两个天线同时接收，即两幅影像的回波是来自同一区域的同一次扫描，地面上的各点在两幅影像中位置能很好的对应起来，这样大大减少了配准的难度，从而可以避免配准产生的误差；由于天线的位置是固定的，两天线的距离即基线也是固定的，这样也减少了基线求解引入的误差。但是这种方式也有缺点，例如需要安装两幅天线和两套接收通道，使得雷达硬件系统变得比较复杂；对导航等硬件的技术水平和对飞机飞行状况都要求较高；而且空间基线b的选择余地较小，受到飞行平台的尺寸限制。（3）重复轨道单天线（repeat-track interferometry，rti）模式重复轨道干涉测量(rti)只要求安装一副天线，通过在不同的时间对同一块地区进行成像，同一区域的成像时间仍然保持一定的相关性，进而实现干涉测量。重复轨道的sar干涉测量(rti)既可以测量地形高程，也常用来监测地表运动。目前星载sar系统可以通过轨道进行重访，即对同一地物进行两次重复观测。这种模式用时间基线来描述时间间隔，理想状态下间隔越短越好，否则会出现严重的去相关现象。与机载系统相比，星载sar系统由于受大气的影响较小，具有准确并且稳定的运行轨道，因此星载系统比机载系统具有更大的优势。但是，大多数的重复轨道干涉测量轨道并不是完全重合的，所以干涉相位信号还包括视线向位移信息和地形信息。差分干涉测量(d-insar)即是去除了地形信息得到的目标运动或形变量的技术。差分干涉测量可监测陆地地表的微小形变，监测时间可从几天到几年，可以获取大范围的、高精度的地面高程变化信息。从平台载体来看insar分为机载和星载两类。机载insar的机动性好，测量精度高，但是地面覆盖能力有限。星载insar具有快速覆盖全球的能力，有利于进行大范围测绘和动态过程的长期监测，特别适合危险地区或人类无法进入地区的研究工作，并且由于其受大气的影响较小从而具有更准确、更稳定的运行轨道，因此具有更大的优势，得到了广泛的应用。在星载insar中大多要求使用相当大的天线(可达10m以上)，如果要在一个航天器上使用两副这样的天线并且具有可观的基线长度，这在空间能力和成本上可能都是不允许的。因此需要两副天线的单轨双天线横向模式和单轨双天线纵向模式在星载insar中的应用受到了限制，所以最常采用的是星载insar的重复轨道单天线模式。2.2.2 insar基本原理合成孔径雷达干涉测量是利用卫星或飞机搭载的合成孔径雷达系统，通过两副天线同时观测（单轨模式），或两次近平行的观测(重复轨道模式)，获取地面同一景观的复影像对。由于目标与两天线位置的几何关系，在复图像上产生了相位差，形成干涉条纹图。干涉条纹图中包含了斜距向上的点与两天线位置之差的精确信息。根据复雷达图像的相位差信息，利用传感器高度、雷达波长、波束视向及天线基线距之间的几何关系，通过成像处理、干涉数据处理和几何转换等来提取地面目标地形的三维信息。本文采用的为重复轨道于涉测量得到的sar数据，因此下面将以重复轨道干涉测量(rti)为例，介绍insar的基本原理。重复轨道干涉测量的卫星轨道与目标的相对几何关系如图32所示。a1、a2是在不同时间同一天线对同一地区进行观测，r1和r2分别为雷达两次成像时天线中心到目标点p的斜距，b是两天线之间的基线，为雷达天线a1对目标点p成像时的侧视角，为基线与水平方向的夹角，卫星传感器高度为h，地面点p高程为z，此外，令r2=r1+r当两副天线所形成的复数影像精确配准后，对应象素值共轭相乘就形成干涉条纹图（interferogram），即有：sar影像1： u1=|u1|ej1sar影像2： u2=|u2|ej2干涉图： uint=|u1|u2|ej1-2接收到的雷达信号的相位主要有两部分组成：一是往返路径确定的相位；一是由地面分辨率单元不同的散射特性造成的随机相位，用公式可表示为：1=-4r1+arg12=-4r1+r+arg2上式中右边第一项表示由于往返路径确定的相位，arg1和arg2表示由地面分辨率单元不同的散射特性造成的随机相位，2 为系数，表示两副天线都发射和接收信号，即反映往返双程的相位差。若是一副天线用于发射信号两幅天线同时接收信息（双天线飞行模式），则系数应该为 1，在干涉图中只反映单程的相位差。一般地， 我们假设在两次获取的 sar 影像中由于散射特性引起的随机相位贡献基本相同，则干涉图的相位就只与雷达信号传播的路径差有关，有=1-2=4r根据图中的几何关系，有：b=bsin-=r1-r2=-r所以，可以得到与平行基线b的关系：=-4b由干涉图得到的只是相位主值，还需要进行相位解缠才能得到真实的相位。如图所示，在三角形s1ps2中，根据余弦定理可得：r1+r2=r12+b2-2r1bcos2-+那么r=-bsin-+b22r1-22r1r1=r2-b22bsin-r由于rr1，br1，因而上式可以近似表示为：r=-bsin-=-b其中，b为传感器s1到s2的距离，称为基线；为连线s1p与竖直方向的夹角，称为视角；为基线b与水平方向的夹角；b为基线b在连线s1p上的投影。干涉相位表示如下：=4r=-4b当地表相对于雷达视线向变化半个波长时，就带来一个2的相位变化，在干涉图上表现为一条完整干涉条纹的出现，所以雷达干涉测量在地表形变监测方面是相当灵敏的。此外还发现，形变增量和雷达视向形变相位成反比的关系，即顺着雷达视线延伸方向的移动对应的形变相位是负值，逆着雷达视线延伸方向的移动导致的形变相位是正值，地面沉降导致的形变相位是负值，而地面上升导致的形变相位为正值。同时根据图中由，p点高程z与传感器高度h有如下关系：z=h-r1cos将上述两式带入高程公式有：z=h-42-b22bsin-+2cos2.3 d-insar基础随着高分辨率sar数据的不断出现和差分干涉技术的逐渐成熟，合成孔径雷达差分干涉（d-insar）技术已经成为监测地表形变的重要手段之一，检测的精度可以达到厘米级甚至毫米量级，其应用领域也在不断扩大如地震、火山、冰川运动、地面沉降、滑波监测、活动板块和构造监测等方面。23.2 d-insar原理差分干涉测量就是应用重复轨道干涉测量的方法来进行地表形变的测量。详细地分析起来，insar 干涉测量的相位主要有几部分组成，可以用下面的公式表示：=flat+topo+defo+atm+noi式中为雷达波单程传播距离差。在干涉基线不为零的情况下，主要包括上式中五个部分。其中flat为平地相位；topo为地形起伏导致的地形相位；defo为两次成像期间地表位移引起的形变相位；atm为两次成像时大气状态不一致引起的相位延迟；noi为残余噪声相位。如果想要获取地表形变信息，则必须要消除区域地形相位、平球相位和大气延迟的相位信息。平球相位可通过卫星轨道和地球椭球严密的几何关系来形成严整的几何算法，大气延迟相位贡献份额是最少的，规律尚且不明显，有专门研究。现阶段，消除区域地形相位方法有四种。2.3.3常规d-insar方法一是利用两幅sar图像，再加上其它的dem数据，基于已有的成像参数模拟干涉条纹图，从而达到消除地形因素的效果，一般称之为“二轨法”（two-pass）。二是利用三副sar图像，采用干涉的方法，消除地形的影响，一般称之为“三轨法”（three-pass）。三是利用四幅sar图像，采用干涉的方法消除地形的影响，一般称之为“四轨法”（four-pass）。（一) 二轨法(two-pass)双轨法是1993年由massonnet等人首先提出的。此方法所用到的数据有：两景sar影像（一景形变前、一景形变后）、研究区域的dem数据。其具体的操作步骤是两景sar影像生成一幅干涉条纹图，研究区域dem模拟成地形相位条纹图，然后从雷达影像生成的干涉条纹图中减去dem模拟的地形相位条纹图，便可得到研究区域的形变信息。一般公式可表示如下： r=4d-dem2轨法的优缺点主要表现在dem数据容易获得，两景雷达影像的干涉图也比较容易做出，所以认为此方法是可行的。但是由于所获取的dem数据精度不高，且干涉基线过长，导致dem的误差直接传递到最后的形变结果中。(二) 三轨法（three-pass）1993年zebker等人提出了利用三轨法从sar干涉图中提取地表形变的方法，此方法需要三景sar影像，两景来自地表形变之前（称为主、辅影像），一景地表形变之后获取（称为从影像），通过对主、从与主、辅两对影像进行干涉处理，得到两幅干涉条纹图。前者只包含了地形信息，后者既有地形信息又有形变信息。然后对两幅干涉条纹图分别进行平地效应去除、相位解缠，最后对获得的两幅干涉图进行差分处理，获取地表的形变信息。三轨法差分干涉测量几何示意图所示，a1为主影像、a2为从影像、a3为辅影像。设定第一幅sar影像与第三幅sar影像成像时间较近，并且这期间没有形变发生，这两幅影像称为地形像对；第一幅和第二幅影像成像时间间隔较长，期间发生了形变，这两幅sar影像对称为形变相对。a1+a2deformation pair（含地形信息和形变信息）a1+a3topographic pair（仅含地形信息）获取的第一幅影像主影像的返回信号可以表示为：u1=|u1|exp-j4r1从影像获取的返回信号可以表示为： u2=|u2|exp-j4r2+r主影像与从影像之间相位差包含地表信息和形变信息，由干涉图相位公式可得：1=-4b12+r=-4b12sin2-2+r假定形变前获取的辅影像返回的信号表示为：u3=|u3|exp-j4r3主、辅影像所获得的干涉图中，既包含了研究区域的地形信息，根据公式有：2=-4b13=-4b13sin1-1由上述两个公式可以得到视线向形变量所引起的相位：defo=1-b12b132=-4r通过干涉条纹图的相位和初始影像的轨道参数计算defo，求得雷达视线向上的地表形变量r。其中基线距的比值为：b12b13=b12sin2-2b13sin1-1是视角的函数，它取决于成像参数和每个点的地形。为了解求出干涉相位的值，我们必须获得研究区域的高程资料。由于视角是由参考视角0和视角增量z组成。对于任意成像方程，0是确定的，随着地形变化改变。为了计算方便和确保结果的精度，我们采取去平地的方法，从干涉纹图中去除平地效应引起的相位变化，从而得到去平后的相位分量为：f1=-4b13sin1+z1-1-sin1-1=-4b13cos1-1z1f2=-4b12sin2+z2-2-sin2-2-4r=-4b12cos2-2z2-4r假设z1=z2，用去平地相位重新表示式可得：defo=f1-b12b13f2=-4r在这个式子中，不需要精确的和地形信息，就可以求解出视线方向形变量r。三轨法的运行过程需要三景sar影像，相对于双轨法它的优点在于不需要额外的dem数据，则所生成的两个干涉像对共用一景sar影像，因此干涉图间的配准处理过程也相对的比较容易。这种方法也存在着缺点，具体体现在相位解缠的精度将会直接影响最终的结果。(三) 四轨法（four-pass）四轨法是获取地表形变前的三幅影像和一副形变后的影像，如图2.7。两幅形变前图像组成了一对干涉对获取地面地形信息，另一副形变前图像与形变后图像组成另一对干涉对，该干涉对不仅包含了地面的地形信息，还包含了地面的形变信息，然后将二者进行差分，得到的差分相位便是形变相位，从而便可计算出地面的形变量。该方法也可以理解为两组主、辅sar影像，其中一对影像用于生成dem，作为另外一组sar影像消除地形相位的工具。当没有外部dem时，利用这种方法可以弥补，如ers-1、ers-2、radarsat、jers、envisat。如果主辅影像以及外部dem的主辅影像是相同的sar系统，成像系统参数就较一致，则配准就会比较容易一些，其实质就是两幅sar影像之间的配准。如果不是同样的雷达系统，那么主影像和外部主影像可能就会有一些差别，但可以运用模拟的sar影像来进行配准。例如ers-1/2的tandem数据便是生成dem较好的数据源。2.4 sbas方法2.4.1基于最小二乘法（ls）的d-insar方法usai(2001，2003)提出的多基线距d-insar方法是利用最小二乘方法求解sar影像集的时序形变。为保证组合影像的相干性，一般要求时间和空间基线都较小。在理想情况下，差分干涉相位可简单地表示为线性模型：x=4vxti+n上式中vx为干涉图象元的形变速率，ti为获取影像的时间间隔，n为随机噪声相位，利用最小二乘法可以获取象元处的时间序列形变情况。首先我们假定有n+1幅同一地区的sar图像，获取时间依次（为t0，.，tn），同时假设每一幅图像至少可以与另一幅图像构成干涉；这意味着每一短基线子集至少由 2幅图像组成。基于以上假设， 则生成的干涉图数量为 m 个， 于是可以推出 m 满足下列不等式（假设 n 为奇数）： 让我们简单总结差分雷达干涉图的主要特征并指出研究问题的关键因素。相应的，考虑位于方位-距离坐标系（x，r）处的第j幅干涉图，假定分别是在ta，tb时刻获得的sar图像产生的，并且已经去除了地形相位部分，则j在（x，r）处的干涉相位可以表示为： 式中为雷达波长，d( tb， x， r)和d( ta， x， r)分别tb和ta时刻相对于参考时刻t0的视线向 ( los) 累积形变量，因而有d( t0， x， r)0；自然的，我们可以用d( ti，x，r) ，i= 1，.，n，来表示我们所要得到的形变时间序列，并设对应的相位为(ti，x，r)，则有。我们将所分析那一点象元的形变量所对应的n个未知相位值用向量表示为： 将从差分干涉图上计算的m个值表示为向量，即： 式(4)可以用下面两个向量来定义： is和ie分别对应生成干涉图像对中从从图像和主图像获得的时间序列，假定主图像和从图像总按时间顺序排列，即iej isj，换言之，有如下等式： 因此表达式定义了含n个未知数的m个等式组成的方程组，用矩阵形式表示如下式： a是一个mn矩阵，时：若，则；否则为0。例如，如果，则a矩阵的前几项形式如下： . 式(8)表明a是一个近似关联矩阵，它直接取决于从可用数据中生成的一系列干涉图。由于该特点，如果所有数据都属于一个单一的小基线集，那么有mn，并且a是一个n阶矩阵。因此，当m=n时，方程组(7)是一个定解方程，当mn时，方程组(7)是一个超定解方程。通常，它的解是可求的，在最小二乘法约束下用矩阵形式表示如下： 其中. 首先，我们发现该方法没有考虑去相关现象，也没有考虑两幅图像获取时由于各层大气折射率变化引起的变化和由于没有精确去除地形相位部分而可能包含的原始相位，我们假定了简单的模型来做当前分析，这是理想情况下的最优解，但ls方法是sbas方法的基础。1.2.4 基于小基线集的d-insar方法bernardino等(2002)和lanari等(2004)提出了利用小基线集(sbas)方法探测地表形变的sbas方法将所有覆盖同一地区的sar影像组成若干个子集，子集内的影像基线距(包括时间基线距和空间基线距)较小，子集间的基线距较大。如果所有影像都属于一个小基线集，则可以利用最小二乘法得到形变相位，但实际上这种可能性很小，所有可用影像通常被分为几个子集。因此，为了增加形变信号的时间采样频率，我们面对的是数据存在于不同的基线集中。显然，根据方程(7)，矩阵a秩亏，ata是一个奇异矩阵。例如，假设有l个不同的小基线集，则a的秩为n-l+1，方程会有无穷个解(我们假设nm；其他例子也作此假设)。一个简单地求解方程(7)的方法是利用svd奇异值分解法。通过svd法求得矩阵a的广义逆，从而给出方程(7)的最小范数意义上的最小二乘解。特别的，我们通过svd将矩阵a分解如下： 此处，u是一个mm的正交矩阵，其前n行是aat的特征向量，称为a的左奇异向量；v是nm的酋矩阵，它的所有行是ata的特征向量，称为a的右奇异向量；s是一个mm矩阵，它的元素(奇异值i)是mm矩阵aat对应特征值的平方根。通常mn，有m-n个特征值为0；此外由于矩阵a的秩亏特性，有l-1个附加的0特征值，总结来看： . 在最小二乘约束下求，可以表示如下： 其中 式中.从而有： 此处ui和vi分别是u和v的行向量。为了得到具有物理意义的可信的解决方法，我们将方程(7)中的未知量转化成时间相邻观测量平均相位速率，相应的，新的未知量为： 式(6)变为： 用矩阵形式表示，可得下式。 b是一个mn矩阵，对一般项(j，k)有b(j，k)=tk+1-tk，当isj+1kiej，其他的d(j，k)=0。当然，对矩阵b进行奇异值分解，可以得到速度矢量v的最小范数解，且结果具有连续性。2.3.3附加约束条件的小基线集d-insar短基线d-insar法为了提高干涉相位的相干性，通常都需要依据空间基线的阈值将影像分成（）组，导致干涉图的主影像不一致，致使方程秩亏，即，其解不唯一，实际上，除了可以应用奇异值分解法求得影像序列间平均相位速率的最小范数最小二乘解，还可以利用形变相位在时间域上的某些约束条件，将未知数转化为这种模型中的参数，从而减少未知数的个数，直接获取模型参数的最小二乘解。假设式中的平均相位变化速率满足如下线性模型： 式中描述了的组成，为模型参数，将其带入得： 一般情况下，参数的个数较少，使得式的系数矩阵非奇异，可以利用最小二乘法直接求解，例如对任意像元，考虑差分相位的如下三次模型： 其中，分别表示形变的平均速度，平均加速度以及平均加速度变化率，则式中的和可以表示成： 实际计算中考虑到大气延迟以及dem误差的影响，式改写成如下形式： 式中，表示dem误差引起的相位，其大小与垂直基线成正比，与卫星到地面间的距离成反比，且与卫星视角有关；为大气延迟相位。第三章 insar数据处理合成孔径雷达干涉测量技术中用来提取相位信息的信息源是合成孔径雷达复数据，而精确测量出雷达图像上各点的三维位置和变化信息的数据依据则是根据雷达传感器高度、雷达波长、波束视向及天线基线距之间的几何关系。而d-insar 技术的发展基础及依据则是 insar 技术，其工作过程则是通过对地表形变前后的干涉图之间的相位差进行比较，目的是用来监测地面目标的位移。在 insar 的数据处理过程中，本节着重介绍了干涉复图像对的配准、去平地效应、干涉图滤波、相位解缠这几步的处理将直接影响输出结果影像的质量以及精度。3.1 影像配准处理在合成孔径雷达干涉测量中，图像对的精确配准是至关重要的。我们知道，insar 技术的工作思想是地表同一目标点在两次成像中的干涉相位差，因此，两次成像的地面点目标位置的掌握对于我们来说是至关重要的。如果两次成像中所各自对应的地面离散样点完全重合，那么就可以利用位置相同的像素构成像素对，进一步求取相对应像素的地面点相位差。实际操作过程中，由于两次飞行是各自进行的，这样就致使成像结果中像素对应的地面样点完全吻合的机率微乎其微。由于两次飞行的独立性，使得地面与各航过的斜平面存在不同的夹角，造成成像几何的差别，导致两次航过成像地面离散样点不可能完全重合。但是我们在影像处理过程中要求其应该完全重合，为了处理的顺利进行，我们就必要要进行某种方法的处理，这样就引入了复影像对的配准。简而言之，就是在空间域上对两幅或者多幅地面离散样点不重合的影像进行配准的过程。在 insar 数据处理过程中配准的质量结果将直接影响获取地面 dem的精度。对于精度至少达到 0.1 像素配准标准我们称之为而高质量的干涉图图像配准。对于未采取任何辐射分辨率改善、纹理模糊、具有严重的斑点噪声等insar 影像来说，高精度自动配准相对于光学影像配准，无论是人工方法还是自动配准方法都存在着很大的困难。而且图像配准本身还涉及到旋转操作、尺度变换和重采样等很多操作问题。整个干涉复图像对的配准分三个步骤来进行：首先是粗配准；其次是像元级配准；最后是亚像元级配准。在文中我们主要考虑后两个步骤，常见的方法是通过人工在主图像和辅图像上选取共同的地区进行像元级配准及亚像元级配准。这种方法需要花费很长的时间去找寻两幅图像的特征区域，从而无法实现图像对的自动配准。3.1.1 粗配准干涉复图像对的粗配准是根据 sar 数据获取时的卫星轨道参数所计算出主图像和辅图像中心的位置，然后为在统一的坐标系下计算出辅图像相对于主图像各点的偏移量。假设主图像上的一点为p，对应的行列值pml，p，而该点在辅图像上对应的行列值为ps(l，p)，则有psl，p=psl，p+offsetl，p，粗配准就是利用卫星的轨道参数来计算图像对的偏移量offsetl，p。其运行中的数学方程有：多普勒方程： xsdx =0斜距方程： dx dx vlighttrange=0椭球方程： x2a2+y2b2+z2c2-1=0上面式子中，xs表示卫星的位置矢量，x 是象元点的位置矢量，dx =x -xs。粗配准的步骤具体可以分为以下三步：步骤一：计算主图像中心点像元pcm(l，p)在轨道坐标系中的坐标值x，y，z。方位时间、地距时间分别利用行数和 pre 值、像元数和 rsr 值来计算，卫星位置的矢量是根据方位、地距时间以及 slc 头文件中的速度、位置矢量利用 cubic spline 方法计算。然后计算中心点像元的坐标值x，y，z，其方法是对方程求偏导利用多次迭代设置收敛阈值；步骤二：基于多普勒方程，计算该点在辅图像上的行列值pcs(l，p)；步骤三：根据p