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

小基线集技术监测北京沉降 Land subsidence Monitoring in Beijing based on SBAS method Supervisor: XXX XXXXX University June XXXX 第一章绪论 1.1问题提出及研究意义 地面沉降广泛存在于全世界的城市中，这是全球在城市化进程中出现的不可忽视的环境地质问题，这种连续且渐进的地质灾害一旦形成很难恢复，不仅给城市带来巨大危害，同时还影响到人类生命财产安危和区域的可持续发展。目前，我国已经有96个大中型城市发生了地面沉降，主要分布在长江三角洲地区，松辽、黄淮海、东南沿海和内陆河谷等平原地区，沉降总面积超48655 km2;其中全国累计沉降量超过200 mm的地区达到7.9万km2，比较严重的地区有上海、天津、台北、西安、宁波、苏州等地，地面累计沉降量达460~2780 mm，地面年沉降速率达10~56 mm/a，已经造成上千亿元的经济损失[1]。当下长江流域、黄河流域、珠三角、松辽平原和环渤海地区及东南沿海城市的地面沉降由于过量开采地下流体资源已经形成巨型漏斗地面沉降在区域上呈扩张趋势，并且对地区基础设施建设造成极坏的影像，其形成和发展速度令世界震惊[2]。 地面沉降现象日益恶化和严重，并在部分地区伴生地裂缝、岩溶塌陷 等环境地质问题，其破坏性较大，给人民生命财产与社会经济发展带来严重威胁，已经成为制约区域环境、社会经济可持续发展的主要地质灾害之一。严重的地面沉降可产生许多危害，特别是一些大中城市可能存在着较大的潜在的危险，主要表现在以下几个方面: ?地面标高损失明显，地面低洼积水难排，河道泄洪能力下降，洪涝灾害加重; ?地面沉降、地裂缝使建筑物基础下沉、 工程 路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理 设施损坏，长时期严重的地面沉降会造成地下排水管道的破坏，严重影响城市生活污水和雨后积水的排放; ?地面沉降导致地下水环境恶化，影响居民生活用水。 为了达到有效预防地面沉降以保证城市可持续发展和为城市发展规划提供 可靠地质资料数据的目的，必须建立长期有效的地面沉降调查与监测机制并制订合理措施对沉降灾害进行预防与治理。常规的地面沉降监测一般包括高等级水准测量、地下精细观测方法(基岩标、分层标等)和GPS网络技术监测等方法，这些方法虽然能够达到较高的地表形变观测精度，但这些方法往往耗费大量的人力物力，观测周期长，实测费用高，同时这些技术仅能获取离散点的形变量，而且在量测大范围的沉降灾害监测时这些技术劣势明显[3]。 下表是几种地面沉降监测方法的比较。 地面沉降给社会经济和人民的生命安全等方面造成巨大的危害，直接影响区域经济的可持续发展，而实时监测地面沉降可以提前做出风险预警、 灾害预报与评估，避免重大的灾难事故的发生;可以为合理开采利用地下水资源、煤、油、气等矿产资源做出指导;可以作为城市的合理规划、发展、立法的依据。因此，实时监测地面的沉降对于维护人民的生命安全和经济的可持续发展具有很高的意义。及时开展 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=7(3)的同震形变场，与测量数据及弹性形变模型获得的结果吻合很好，这一结果发表在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数据进行 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 比较，得到较一致的形变场和形变速率。 为了抑制失相关和大气影响对形变信息提取的负面影响，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影像处理软件 1.2.3 几种典型的星载SAR卫星系统 国际上主要卫星SAR系统及相关参数 第二章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)单轨双天线横向(Cross–track interferometry，XTI)模式 图3(1是SAR交轨干涉测量的示意图。交轨干涉测量是利用一副天线发射雷达波，由两幅天线同时接受回波。由于两幅天线接收的回波具有一定的相干性，两幅天线与地面目标物之间的路径差为影像的相位差。相位差是由两幅天线与地面点之间的路径差造成的，路径差是与地面点的高程相联系的。因此，获取干涉测量系统的几何参数，就可以将相位信息转换为高程信息。即除了图像点的方位向和距离向的坐标外，相位差提供了高程信息，这样就能够重建地面目标点的三维特征。 由于进行SAR交轨干涉测量需要同时安装两幅SAR天线，因此通常应用于机载系统或航天飞机任务。 (2)单轨双天线纵向(Along-track interferometry，ATI)模式 SAR顺轨干涉测量技术是由Goldstein和Zebker在1987年提出的。实现方 法也是在一般的机载SAR系统上加装另一副天线，但与XTI不同的是，这两副天线是顺着平台飞行轨道来安装，前后基线之间的长度通常为2-20m。假设天线之间的间距是20m，飞机的飞行速度是200m/s，则两副天线所观测到同一地物的时问间隔是0.1s，在这个时间间隔内时钟漂移和传播延迟的变化都可以忽略不计的，因此测量到的相位差反映的是地面点的高程变化信息。这项技术被称为“SAR顺轨干涉测量”(ATI)。安装的同 一架飞机上的两幅天线是沿着飞行方向的。ATI常用于水流制图、动目标 检测 工程第三方检测合同工程防雷检测合同植筋拉拔检测方案传感器技术课后答案检测机构通用要求培训 以及定向波谱的测量。 上面两种方式都属于单轨双天线干涉，这种方式的优点有: 1由于可以同时获取两幅SAR影像，在两幅影像中地物具有相同的的后向? 2这种方式是一个天线发送，由两个天线散射系数，获取数据的相干性比较好;? 同时接收，即两幅影像的回波是来自同一区域的同一次扫描，地面上的各点在两幅影像中位置能很好的对应起来，这样大大减少了配准的难度，从而可以避免配 3由于天线的位置是固定的，两天线的距离即基线也是固定的，准产生的误差;? 这样也减少了基线求解引入的误差。 但是这种方式也有缺点，例如需要安装两幅天线和两套接收通道，使得雷达硬件系统变得比较复杂;对导航等硬件的技术水平和对飞机飞行状况都要求较高;而且空间基线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的基本原理。重复轨道干涉测量的卫星轨道与目标的相对几何关系如图3(2所示。A1、A2是在不同时间同一天线对同一地区进行观测，??1和??2分别为雷达两次成像时天线中心到目标点P的斜距，B是两天线之间的基线，θ为雷达天线A1对目标点P成像时的侧视角，α为基线与水平方向的夹角，卫星传感器高度为h，地面点P高程为z，此外，令??2=??1+Δ?? 当两副天线所形成的复数影像精确配准后，对应象素值共轭相乘就形成干涉条纹图(interferogram)，即有: SAR影像1:??1=|??1|??????1 SAR影像2: ??2=|??2|??????2 干涉图:????????=|??1||??2|???? ??1???2 接收到的雷达信号的相位主要有两部分组成:一是往返路径确定的相位;一是由地面分辨率单元不同的散射特性造成的随机相位，用公式可表示为: ??1=? ??2=?4????+?????? ??1 14?? ??1+Δ?? +?????? ??2 上式中右边第一项表示由于往返路径确定的相位，?????? ??1 和?????? ??2 表示由地面分辨率单元不同的散射特性造成的随机相位，2 为系数，表示两副天线都发射和接收信号，即反映往返双程的相位差。若是一副天线用于发射信号两幅天线同时接收信息(双天线飞行模式)，则系数应该为 1，在干涉图中只反映单程的相位差。 一般地，我们假设在两次获取的 SAR 影像中由于散射特性引起的随机相位贡献基本相同，则干涉图的相位就只与雷达信号传播的路径差有关，有 ?=??1???2= 根据图中的几何关系，有: =Bsin ???α =??1???2=?ΔR ??? 所以，可以得到?与平行基线???的关系: ?=?4???? ?4??Δ?? 由干涉图得到的只是相位主值，还需要进行相位解缠才能得到真实的相位?。 如图所示，在三角形S1PS2中，根据余弦定理可得: ?? ??1+Δ?? 2=??12+??2?2??1Bcos ?θ+α 那么 ??2Δ??2 Δ??=?Bsin ???α +?11Δ??2???2 ??1=由于ΔR???1，B???1，因而上式可以近似表示为: Δ??=?Bsin ???α =???? 其中，B为传感器S1到S2的距离，称为基线;??为连线S1P与竖直 方向的夹 角，称为视角;α为基线B与水平方向的夹角;???为基线B在连线S1P上的投影。 干涉相位表示如下: ?=4??4??Δ??=??? ? 当地表相对于雷达视线向变化半个波长时，就带来一个2π的相位变化，在干涉图上表现为一条完整干涉条纹的出现，所以雷达干涉测量在地表形变监测方面是相当灵敏的。此外还发现，形变增量和雷达视向形变相位成反比的关系，即顺着雷达视线延伸方向的移动对应的形变相位是负值，逆着雷达视线延伸方向的移动导致的形变相位是正值，地面沉降导致的形变相位是负值，而地面上升导致的形变相位为正值。 同时根据图中由，P点高程z与传感器高度h有如下关系: z=h???1cos?? 将上述两式带入高程公式有: ?4π?B2 2Bsin θ?α +2πλ?2z=h? cosθ 2.3D-InSAR基础 随着高分辨率SAR数据的不断出现和差分干涉技术的逐渐成熟，合成孔径雷达差分干涉(D-InSAR)技术已经成为监测地表形变的重要手段之一，检测的精度可以达到厘米级甚至毫米量级，其应用领域也在不断扩大如地震、火山、冰川运动、地面沉降、滑波监测、活动板块和构造监测等方面。 2(3.2D-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模拟的地形相位条纹图，便可得到研究区域的形变信息。一般公式可表示如下: ?? ?????dem ?? 2轨法的优缺点主要表现在DEM数据容易获得，两景雷达影像的干涉图也比较ΔR=容易做出，所以认为此方法是可行的。但是由于所获取的DEM数据精度不高，且干涉基线过长，导致DEM的误差直接传递到最后的形变结果中。 (二)三轨法(three-Pass) 1993年Zebker等人提出了利用三轨法从SAR干涉图中提取地表形变的方法，此方法需要三景SAR影像，两景来自地表形变之前(称为主、辅影像)，一景地 表形变之后获取(称为从影像)，通过对主、从与主、辅两对影像进行干涉处理，得到两幅干涉条纹图。前者只包含了地形信息，后者既有地形信息又有形变信息。然后对两幅干涉条纹图分别进行平地效应去除、相位解缠，最后对获得的两幅干涉图进行差分处理，获取地表的形变信息。 三轨法差分干涉测量几何示意图所示，A1为主影像、A2为从影像、A3为辅影像。设定第一幅SAR影像与第三幅SAR影像成像时间较近，并且这期间没有形变发生，这两幅影像称为地形像对;第一幅和第二幅影像成像时间间隔较长，期间发生了形变，这两幅SAR影像对称为形变相对。 A1+A2———————————deformation pair(含地形信息和形变信息) A1+A3———————————topographic pair(仅含地形信息) 获取的第一幅影像主影像的返回信号可以表示为: 4????1=|??1|exp ????? 1 从影像获取的返回信号可以表示为: 4?? ??2+Δr ??2=|??2|exp ???主影像与从影像之间相位差包含地表信息和形变信息，由干涉图相位公式可得: 4??124??12??1=? ???+Δr =? ??sin ??2???2 +Δr 假定形变前获取的辅影像返回的信号表示为: 4????3=|??3|exp ????? 3 主、辅影像所获得的干涉图中，既包含了研究区域的地形信息，根据公式有: ??2=?4??134??13???=???sin ??1???1 12???4??=??1???2=?Δr ???由上述两个公式可以得到视线向形变量Δρ所引起的相位: ????efo 通过干涉条纹图的相位和初始影像的轨道参数计算????efo，求得雷达视线向上的地表形变量Δr。 其中基线距的比值为: 12?????12sin ??2???2 = ???11是视角θ的函数，它取决于成像参数和每个点的地形。为了解求出干涉相位的值，我们必须获得研究区域的高程资料。由于视角θ是由参考视角??0和视角增量Δ????组成。对于任意成像方程，??0是确定的，Δθ随着地形变化改变。为了计算方便和确保结果的精度，我们采取去平地的方法，从干涉纹图中去除平地效应引起的相位变化，从而得到去平后的相位分量为: 4??13????1=??? sin ????1+Δ????1???1 ?sin ????1???1 4??13=???cos ????1???1 Δ????1 4??124??=??? sin ????2+Δ????2???2 ?sin ????2???2 ?Δr 4??124??=???cos ????2???2 Δ????2?Δr ????2 假设Δ????1=Δ????2，用去平地相位重新表示式可得: ????efo12???4??=????1?????2=?Δr ???在这个式子中，不需要精确的θ和地形信息，就可以求解出视线方向形变量Δ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 =4????? ?? ?????+???? 上式中?? ?? 为干涉图象元的形变速率，????为获取影像的时间间隔，????为随机噪声相位，利用最小二乘法可以获取象元处的时间序列形变情况。 首先我们假定有N+1幅同一地区的SAR图像，获取时间依次(为t0，....，tn)，同时假设每一幅图像至少可以与另一幅图像构成干涉;这意味着每一短基线子集至少由 2幅图像组成。基于以上假设，则生成的干涉图数量为M个，于是可以推出M满足下列不等式(假设N为奇数): N?1N?1?M?N()22 让我们简单总结差分雷达干涉图的主要特征并指出研究问题的关键因素。相应的，考虑位于方位-距离坐标系(x，r)处的第j幅干涉图，假定分别是在tA，tB时刻获得的SAR图像产生的，并且已经去除了地形相位部分，则j在(x，r)处的干涉相位可以表示为: ??j(x,r)??(tB,x,r)??(tA,x,r) 4??d(tB,x,r)?d(tA,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)，则有?(ti,x,r)?4?d(ti,x,r)。 我们将所分析那一点象元的形变量所对应的N个未知相位值用向量表示为: ?T???(t1),...,?(tN)? 将从差分干涉图上计算的M个值表示为向量，即: ??T????1,...,??M? 式(4)可以用下面两个向量来定义: IS??IS1,...,ISM?IE??IE1,...,IEM? IS和IE分别对应生成干涉图像对中从从图像和主图像获得的时间序列，假定主图像和从图像总按时间顺序排列，即IEj>ISj，?j?1,...,M.换言之，有如下等式: ??j??(t)??(tIS),IEjj?j?1,...,M. 因此表达式定义了含N个未知数的M个等式组成的方程组，用矩阵形式表示如下式: A???? IS?0A?j,ISj???1,A(j,IEj)??1?j?1,...,M时:A是一个M×N矩阵，若j，则; 否则为0。例如，如果??1??4??2,??2??3??0，则A矩阵的前几项形式如下: ?0?10?1?00?10A???............??...............?...??...??...?. 式(8)表明A是一个近似关联矩阵，它直接取决于从可用数据中生成的 一系列干涉图。由于该特点，如果所有数据都属于一个单一的小基线集，那么有M?N，并且A是一个N阶矩阵。因此，当M=N时，方程组(7)是一个定解方程，当M>N时，方程组(7)是一个超定解方程。通常，它的解是可求的，在最小二乘法约束下用矩阵形式表示如下: #TT??A#????A?AAA?其中. ?1 首先，我们发现该方法没有考虑去相关现象，也没有考虑两幅图像获取时由于各层大气折射率变化引起的变化和由于没有精确去除地形相位部分而可能包含的原始相位，我们假定了简单的模型来做当前分析，这是理想情况下的最优解，但LS方法是SBAS方法的基础。 1.2.4 基于小基线集的D-InSAR方法 Bernardino等(2002)和Lanari等(2004)提出了利用小基线集(SBAS)方法探测地表形变的SBAS方法将所有覆盖同一地区的SAR影像组成若干个子集，子集内的影像基线距(包括时间基线距和空间基线距)较小，子集间的基线距较大。 如果所有影像都属于一个小基线集，则可以利用最小二乘法得到形变相位，但实际上这种可能性很小，所有可用影像通常被分为几个子集。因此，为了增加形变信号的时间采样频率，我们面对的是数据存在于不同的基线集中。显然，根据方程(7)，矩阵A秩亏，ATA是一个奇异矩阵。例如，假设有L个不同的小基线集，则A的秩为N-L+1，方程会有无穷个解(我们假设N?M;其他例子也作此假设)。 一个简单地求解方程(7)的方法是利用SVD奇异值分解法。通过SVD法求得矩阵A的广义逆，从而给出方程(7)的最小范数意义上的最小二乘解。 特别的，我们通过SVD将矩阵A分解如下: A?USVT 此处，U是一个M×M的正交矩阵，其前N行是AAT的特征向量，称为A的左奇异向量;V是N×M的酋矩阵，它的所有行是ATA的特征向量，称为A的右奇异向量;S是一个M×M矩阵，它的元素(奇异值бi)是M×M矩阵AAT对应特征值的平方根。通常M>N，有M-N个特征值为0;此外由于矩阵A的秩亏特性，有L-1个附加的0特征值，总结来看: S?diag??1,...,?N?L?1,0...,0?. 在最小二乘约束下求Φ，可以表示如下: ??A?????T?其中A?VSU 式中 ???S??diag(1,...,N?L?1,0,...,0).从而有: N?L?1?i?1??Tui?ii 此处ui和vi分别是U和V的行向量。 为了得到具有物理意义的可信的解决方法，我们将方程(7)中的未知量转化成时间相邻观测量平均相位速率，相应的，新的未知量为: v??v1?T?1 t1?t0,...,v??N??N?1??tN?tN?1? 式(6)变为: k?ISj?1?(tIEjk?tk?1)vk???j,?j?1,...,M 用矩阵形式表示，可得下式。 Bv??? B是一个M×N矩阵，对一般项(j，k)有B(j，k)=tk+1-tk，当ISj+1?k?IEj，?j?1,...,M，其他的D(j，k)=0。当然，对矩阵B进行奇异值分解，可 以得到速度矢量v的最小范数解，且结果具有连续性。 2.3.3附加约束条件的小基线集D-InSAR 短基线D-InSAR法为了提高干涉相位的相干性，通常都需要依据空间基线的阈值将影像分成L(L?1)组，导致干涉图的主影像不一致，致使方程错误～未找到引用源。秩亏，即R(B)?N?L?1，其解不唯一，实际上，除了可以应用奇异值分解法求得影像序列间平均相位速率的最小范数最小二乘解，还可以利用形变相位在时间域上的某些约束条件，将未知数转化为这种模型中的参数，从而减少未知数的个数，直接获取模型参数的最小二乘解。 假设式错误～未找到引用源。中的平均相位变化速率v满足如下线性模型: v?Mp 式错误～未找到引用源。中M描述了v的组成，p为模型参数，将其带入错误～未找到引用源。得: BMp??? 一般情况下，参数p的个数较少，使得式错误～未找到引用源。的系数矩阵非奇异，可以利用最小二乘法直接求解，例如对任意像元，考虑差分相位的如下三次模型: ?(ti)?(ti?t0)??(ti?t0)2???(ti?t0)31 216 其中，，?分别表示形变的平均速度，平均加速度以及平均加速度变化率，则式错误～未找到引用源。中的p和M可以表示成: pT?[,,?]T ?t1?t01?2??t2?t1?2t0?1M??2?????tN?tN?1?2t01?2?(t1?t0)3?(t0 ?t0)3??6(t1?t0)?33?(t2?t0)?(t1?t0)?6(t2?t1)????33(tN?t0)?(tN?1?t0)??6(tN?t N?1)? 实际计算中考虑到大气延迟以及DEM误差的影响，式错误～未找到引用源。改写成如下形式: BMp?c??h??A??? ，c??h表示DEM误差式中cT?[(4?/?)(B?1/rsin?),?,(4?/?)(B?M/rsin?)] B?h引起的相位，其大小与垂直基线?成正比，与卫星到地面间的距离r成反比， 且与卫星视角?有关; ?A为大气延迟相位。 第三章InSAR数据处理 合成孔径雷达干涉测量技术中用来提取相位信息的信息源是合成孔径雷达复数据，而精确测量出雷达图像上各点的三维位置和变化信息的数据依据则是根据雷达传感器高度、雷达波长、波束视向及天线基线距之间的几何关系。而D-InSAR 技术的发展基础及依据则是 InSAR 技术，其工作过程则是通过对地表形变前后的干涉图之间的相位差进行比较，目的是用来监测地面目标的位移。在 InSAR 的数据处理过程中，本节着重介绍了干涉复图像对的配准、去平地效应、干涉图滤波、相位解缠这几步的处理将直接影响输出结果影像的质量以及精度。 3.1 影像配准处理 在合成孔径雷达干涉测量中，图像对的精确配准是至关重要的。我们知道，InSAR 技术的工作思想是地表同一目标点在两次成像中的干涉相位差，因此，两次成像的地面点目标位置的掌握对于我们来说是至关重要的。如果两次成像中所各自对应的地面离散样点完全重合，那么就可以利用位 置相同的像素构成像素对，进一步求取相对应像素的地面点相位差。实际操作过程中，由于两次飞行是各自进行的，这样就致使成像结果中像素对应的地面样点完全吻合的机率微乎其微。由于两次飞行的独立性，使得地面与各航过的斜平面存在不同的夹角，造成成像几何的差别，导致两次航过成像地面离散样点不可能完全重合。但是我们在影像处理过程中要求其应该完全重合，为了处理的顺利进行，我们就必要要进行某种方法的处理，这样就引入了复影像对的配准。简而言之，就是在空间域上对两幅或者多幅地面离散样点不重合的影像进行配准的过程。在InSAR 数据处理过程中配准的质量结果将直接影响获取地面 DEM的精度。对于精度至少达到 0.1 像素配准标准我们称之为而高质量的干涉图图像配准。对于未采取任何辐射分辨率改善、纹理模糊、具有严重的斑点噪声等InSAR 影像来说，高精度自动配准相对于光学影像配准，无论是人工方法还是自动配准方法都存在着很大的困难。而且图像配准本身还涉及到旋转操作、尺度变换和重采样等很多操作问题。 整个干涉复图像对的配准分三个步骤来进行:首先是粗配准;其次是像元级配准;最后是亚像元级配准。在文中我们主要考虑后两个步骤，常见的方法是通过人工在主图像和辅图像上选取共同的地区进行像元级配准及亚像元级配准。这种方法需要花费很长的时间去找寻两幅图像的特征区域，从而无法实现图像对的自动配准。 3.1.1 粗配准 干涉复图像对的粗配准是根据 SAR 数据获取时的卫星轨道参数所计算出主图像和辅图像中心的位置，然后为在统一的坐标系下计算出辅图像 相对于主图像各点的偏移量。假设主图像上的一点为P，对应的行列值???? ??，?? ，而该点在辅图像上对应的行列值为??则有???? ??，?? =??粗??(??，??)，?? ??，?? +???????????? ??，?? ，配准就是利用卫星的轨道参数来计算图像对的偏移量???????????? ??，?? 。其运行中的数学方程有: 多普勒方程: ????????? =0 斜距方程:???? ????? – v?????????????????????? =0 椭球方程:??2+??2+??2?1=0 上面式子中， ????表示卫星的位置矢量，?? 是象元点的位置矢量，???? =?? ?????。 粗配准的步骤具体可以分为以下三步: 步骤一:计算主图像中心点像元??cm(??，??)在轨道坐标系中的坐标值 ??，??，?? 。方位时间、地距时间分别利用行数和 PRE 值、像元数和 RSR 值来计算，卫星位置的矢量是根据方位、地距时间以及 SLC 头文件中的速度、位置矢量利用 cubic spline 方法计算。然后计算中心点像元的坐标值 ??，??，?? ，其方法是对方程求偏导利用多次迭代设置收敛阈值; 步骤二:基于多普勒方程，计算该点在辅图像上的行列值??cs(??，??); 步骤三:根据??cm(??，??)和??cs(??，??)来计算主图像和辅图像之间的偏移量???????????? ??，?? 。 x2??2??2 3.1.2 像元级配准 像元级配准其精度大约为 1 个像元，在文中我们所采用的方法是基于窗口的自动配准方法，其配准的运行时在主、辅图像的频率域中进行的。窗口匹配方法的匹配窗口一般选择在主图像上大小可根据实际情况具体 设置，通常用到的窗口大小有 256 像素×256 像素或者 128 像素×128 像素，而搜索窗口则通常选择在辅图像上，运行过程中搜索窗口按照整像元偏移量在不同的行列之中来计算两个窗口的匹配质量评价指标。示意图 3.1 像元级配准需利用相应的评价参数来评价其匹配质量，目前较常用的评价方法有以下四种: 1. 频率极大法:以寻求两幅复图像最大频率为基准、以相干班的稀少分布为前提，且频率小、复值大，则可以表示大频率; 2. 平均波动函数法:以寻求两幅复图像相位差后的某点与在两个方向相邻点的分别相位差的和为最小为基准、以相位噪声的最小或者相位畸异点稀少分布作为前提; 3. 相干系数法:以寻求两幅图像相干系数最大作为基准、以零作为均值、高斯分布为前提，此方法是最为常用的评价方法。 4. 最小二乘法:以相位差平方和最小为基准，同时要考虑几何畸变，以相位畸变点的稀少分布作为前提，但是几何畸变不一定遵循二维的几何变换。 由于文中实验室对整幅图像来进行配准，所以我们选择了比较大匹配窗口和搜索窗口。由于相干估算器所造成的影响，在窗口搜索方法中所计算出来的主、辅图像的相对偏移量估算值并不是完全可靠的。因此在实验过程中要估算相对偏移量的大小，不能只是依据最大相干，还需要剔除明显的误差，来确保多数点相对偏移量。 3.1.3 亚像元级的精配准 文中亚像元级精确配准也是在频率域中实现的，配准精度要求达到1/8个像元，其实现方法与像元级配准类似，具体的操作步骤可以分为以下几步: 步骤 1:首先对主、辅图像作 oversample 处理，然后再利用插值方法对图像对进行插值。目前我们常用的插值方法与方位和斜距方向是无关的，插值的间隔为 0.01 个像元。 步骤 2:进行最大相干估算。其运用的方法是基于窗口的搜索方法。它的计算偏移量的过程像元级配准相似，但也有不同。其不同在于，数据拟合所选用的窗口数量和大小不一样。相对于像元级配准的窗口，亚像元级配准的窗口相对较小，而为了防止最大相干偏差的出现，所采取的措施是适当的增大了搜索窗口的大小，并增加了计算次数。而估算值的可靠性则是根据相对偏移量的稳定性来判断。如果偏移量保持稳定则估算值可靠，反之则不可用，应该重新搜索。 步骤 3:偏移量拟合与辅图像的重采样。文中匹配数据拟合的方法是最小二乘法，其运用了二阶多项式来计算复图像对的坐标转换关系。根据转换关系对辅图像进行重采样。 如果有SAR成像区域的DEM，而且能够轻易的和精确的在SAR影像与DEM之间找到一些共同的地面控制点，利用这些控制点进行SAR图像配准是摄影测量和遥感领域非常成熟的方法。但是如何在SAR影像中辨别和选取地面控制点是该配准方法的难点，较好的做法是在影像所在区域布设一些角反射器，用GPS精确测定它们的坐标和高程，在SAR图像中，这些角反射器的强度信号远远高于背景的强度，因此可以很方便的将这些角反射 器选取出来。如果没有角反射器，那么只有依靠经验来选取地面控制点了，这样的做法误差较大，而且比较费时，往往得不到好的结果。利用地面控制点可以很快实现雷达图像的粗配准，这在有些情况下还是有其特殊的作用的。 3.2重采样 获取了从图像相对于主图像的配准参数之后，需要对从图像进行重采样，使从图像上的每一个像元值精确的对应于主图像上的每一个像元。重采样的方法通 常有最近临法(nearest neighbor)、线性内插法(1inear interpolation)、(4点，6点)立方体卷积插值法(cubic convolution)、(6，8或16点)有限长sinc函数法(truncated sinc)。对于ERS-1/2卫星，6点立方卷积插值可以提供较好的采样精度，但是计算较费时间。 3.3干涉图生成、去平地效应和相干图生成 3.3.1干涉图生成 如第二章提到的，干涉图是由两幅SAR影像进行 实际上，由于SAR图像在距离向和方位向上的分辨率不同，在干涉后一般需将方位向像元做多视处理，使得干涉图符合实际的纵横比例。ERS系列的SAR图像通常选取5:l的多视比例，这是因为其距离向分辨率约为20m左右，而方位向分辨率约为4m左右，那么多视处理后的干涉图分辨率约为20m×20m左右。 3.3.2去平地效应 所谓的“平地效应”是指在干涉纹图中，理论上高程不变的平地不存 在相位差所产生的干涉条纹，而实际的干涉纹图中，受 InSAR 系统本身设计结构的关系，在高程不变的干涉纹图上相位差发生了变化，进而使得干涉纹图相位也产生了偏移。通常我们便把上述所说的这种周期性的变化影响称为“平地效应”。消除地形相位信号中平地效应的过程称之为“去平地效应”。“平地效应”为干涉纹图的解释带来了很大的不变，会让人们对相位差做出错误的判读，所以必须要消除这种影响。目前常见的平地效应消除方法有:l、通过估算距离向的占优势的条纹频率来计算平地相位，然后去除平地效应。该方法对平坦地区最有效，但对于地形变化剧烈的地区，部分地形变化造成的条纹频率和平地条纹频率相同，也被去除，造成错误。2、基于轨道参数和成像区域中心点的位置来去除平地效应。 3、利用粗精度DEM数据去除平地效应。 3.3.3相干图生成 在得到干涉纹图后，需要对相位数据进行分析，以确定在哪些区域相位是一致的，在哪些区域相位不一致，根据缠绕相位的质量来确定相位解缠策略。一般情况下，用相位质量图来评价相位数据的质量。相位质量图主要有4种:相干图、伪相干图、相位倒数变化图和最大相位梯度图。在这些相位质量图中，应用最多的是相干图。相干图是由Parti等人于1993年提出的一种度量SAR影像变化质量的方法。两幅SAR影像??1和??2，的相干性在理论上定义如下: γ=? ?? ????2211在这里，E代表干涉相位期望值，实际上是采样的平均值。但是在实际过程中，不可能获得干涉的数学期望。因此，一般采用一定区域内像素的相位算数平均值来估计数学期望值。于是γ的估计值?? ， 表示如下: ?? =??? ?? ?? ??，?? ?? ??，?? ????2??2 ????=1 ??=1|??1 ??，?? | ??=1 ??=1|??2 ??，?? |其中，M和N为计算相干系数时采用的窗口大小，计算出来的相干系数值在0(完全失相干)和1之间(无噪声)。相干系数的值越大，表明这些点上的干涉图质量越好，反之越差。相干系数可以在计算干涉图之后计算。其对应的图称为相干系数图。 相干系数与信噪比(SNR)之间有对应关系: SNRγ=相干系数越高，信噪比越高。 3.3.4干涉图滤波 雷达成像过程中，雷达系统本身可产生噪声。而且由于不同地物反射电磁渡，在影像的同一像素中产生叠加，也会产生斑点噪声。所以在干涉图行解缠前，还 需要进行滤波，来击除噪声，便于进行解缠计算。进行滤波的主要原因不仅要提高信噪比，为解缠计算提供较好的数据，而且要在去除干涉图噪声的基础上保持边缘信息，防止产生误差。滤波的方法也较多。总的可分为成像前的去躁处理和成像后的滤波去噪。成像后滤波的方法有:空间域滤波、频率域滤波，几何滤波和能够去除高程不变区域噪声的同时保留图像的边缘、纹理等细节信息的基于小波变换的滤波算法，等等。 成像前的滤波操作，可简称为前置滤波，是在生成干涉图钱对原始的附属干涉图像进行滤波。往往从图像对频率特性出发，针对频谱进行滤波。根据机理不同和方法分类，分成距离向和方位向滤波。 后置滤波是在生成干涉图之后对干涉图进行滤波，分为空间滤波和频域滤波。空间域滤波是直接以图像为处理对象，不做任何变换，利用各种图像平滑模板对图像进行卷积处理，以达到压抑或消除噪声的目的。频域滤波是将原始图像进行快速傅里叶变换(fast Fourier transform，FFT)，根据斑点噪声和信号的频谱特性进行滤波。由于快速傅里叶变换耗时长，还会产生低频效应，逐渐被小波变换取代。常用的滤波方法是基于能量谱的滤波方法，该方法是一种局部自适应滤波方法，其结果是增强了滤波窗口中的最强信号部分。实质上，该滤波算子是增强了滤波窗口内最强谱的波长，同时还保留了那些高相位梯度区域;而不像某些线性滤波算子，对这些区域进行平均处理。同样地，该滤波算子由于无信号(不相干)区域内的短波部分具有高能量而对其也进行了保留，使用较多的滤波方法还有低通滤波等。 3.3.5相位解缠 对相位图进行滤波来去除噪声的影响后，就可以对相位进行解缠。相位解缠是InSAR数据处理过程中关键的步骤之一，相位解缠质量的好坏直接影响地面沉降监测的精度。相位差值(或相位条纹)是经过2π周期折叠后的主值，从已有的相位差值恢复出折叠前的原始相位值的处理过程，称之为“相位解缠”。相位解缠是制约InSAR精度的一个瓶颈，主要因为SAR侧视成像方式以及地形起伏引起的图像几何畸变(雷达阴影、透视收缩、叠掩);干涉相位信号的信噪比太低; 地表不连续导致干涉相位存在显著跳跃。 从数学上讲，给定一个二维相位矩阵??i，j，为了解缠该矩阵，需要对 每个点 i，j 加上或者减去2π的整数倍从而得到一个连续的函数??i，j，即: 0，???1 ，??? 0，???1 ??i，j=????，??+2??????，?? i? 上式中????，??是整数，且?π<????，??<??，相位解缠的实质思想就是利用干涉条纹线，从某一参考点对每一个积分路径上的相位缠绕点进行加、减2π的整数倍，该整数值取决于缠绕点与参考点之间的条纹数目。由于多数干涉图中存在着相位的不连续和噪声，所以需要采取一定的数学方法来进行处理，使解缠效果尽可能的准确。 现有的相位解缠方法主要可分为两类，即路径跟踪类法(Path Following)和最小二乘类法(Least Square)。路径跟踪类法基于像元到像元的局部运算来解缠，而最小二乘类法是通过使解绕后与解绕前相位的梯度差整体上最小来进行求解。路径跟踪类主要有四种方法用来处理相位解缠，分别称为枝状切口(Branch Cutting)，质量图法(Quality Algorithm)、掩膜切口法(Mask Cut Algorithm)和Flynn的最小不连续法。最小二乘类有FFT/DCT法、PCG法、多级网格法和最小Lp-norm法等，其他方法还有条纹检测(Fringe detection)、网格自动化(Cellular-Automata)和知识介入(Knowledge Injection) 等。 3.3.6地理编码 地理编码是把相位解缠后的坐标转为椭球坐标系下的坐标。这一步的关键是确定解缠后影像上每个像元对虚的三维坐标。地理编码一般有两种方法:地面控制点法和的方法。地面控制点法顾名思义(是采用解缠后的影像上已知的控制点坐标推算其他地面点的坐标:基于成像参数的方法是 根据SAR的系统参数、影像成像时几何关系和影像分辨率来推算各点的坐标。通过InSSR的一系列流程，即干涉图像的配准、干涉图滤波、去除平地效应、相位解缠及相位到高程的转换，就可得到高精度的数字高程模型。 第四章 D-InSAR方法处理北京地面沉降 4.3 北京地区沉降概况 北京地面沉降最早于1935年在西单到东单一带被发现，截至1952年的17年间最大累计沉降量仅为58mm。上世纪五六十年代，随着东郊地区电子工业区、纺织工业区等迅速发展，地下水大量开采，逐步形成东郊沉降区。上世纪七十年代是北京市地面沉降快速发展时期，沉降区集中，沉降量大。东郊沉降区快速发展的同时，昌平区、顺义区、大兴区等均出现新的沉降中心，个别监测点的沉降量甚至达到81mm。80年代初至90年代初，随着顺义区第八水厂自来水引入市区，并采取了大力开展节约用水和加强地下水开发管理等措施，城区地面沉降明显减缓，但在东郊边缘地带及远郊区地下水集中开采区，地面沉降区面积仍在迅速扩大，并形成了新的通州沉降区。进入21世纪以来，北京平原区地面沉降仍处于快速发展阶段，随着城市规模的不断扩大，加上近十年的干旱少雨，北京市平原区广泛存在着地面沉降灾害。 据统计，截止到2011年底，北京地区沉降量大于100毫米的地区面积已经超过3969平方公里，依据历年累计沉降量和近几年年度沉降情况，将北京平原区地面沉降区域划分为北部沉降区和南部沉降区。北部沉降区包括昌平沙河—八仙庄、朝阳来广营、东郊八里庄—大郊亭(三间房、通州城区、黑庄户—台湖)沉降区，以及顺义平各庄沉降区。南部沉降区主 要分布在大兴榆垡,礼贤地区。昌平沙河—八仙庄中心区最大累积沉降量达到1302毫米，朝阳区来广营沉降区和大兴榆垡—礼贤沉降区的最大累计沉降量均超过1000mm，北京平原区部分地区近几年地面沉降发展迅速，年沉降速率超过50mm/a，局部地区连续三年沉降速率超过100mm/a，部分地面沉降快速发展区内的最大沉降速率均超过128.2 mm/a，早期形成的老沉降区依然在缓慢下沉，且逐渐向东、东北、北、南移动，形成了新的沉降区。 根据北京市地面沉降危害分区参考标准，城市建设发展相对较缓的郊区县，年沉降速率30~50mm/a或累计沉降量500~1000mm区域为地面沉降发育较严重区;大于50mm/a或1000mm区域为严重区。市区内及重点城镇区，年沉降 10~30mm/a或累计沉降量300~500mm区域为地面沉降发育较严重区;大于30mm/a或500mm区域为严重区。北京市地面沉降较严重区和严重区面积不断扩大，由最初发展时期的几十km2扩大到2009年底的1225.41 km2，覆盖面涉及海淀、朝阳、昌平、顺义、通州、大兴等区县。据监测，在朝阳区金盏乡、平谷区、房山区、石景山区及海淀区西北角、昌平区西北部、顺义区北部地区、通州区南部地区，近年来地面沉降发展态势迅猛，地面沉降速率有逐年加大的趋势。沉降区内的工厂、居民楼、地下管道等遭受了不同程度的破坏。 4.2实验区数据介绍 4.2.1 SRTM DEM 数据 在差分干涉测量中，一般需要用到高精度的 DEM，如果无法获取高 精度 DEM，可以使用 USGS NIMA(National Imagery and Mapping Agency)免费发布的 SRTM DEM 数据。 2000 年 2 月 11 日，由美国宇航局(NASA)、美国国家影像与测绘署(NIMA)、德国空间局等单位联合进行了第一次载人航天飞行的雷达地形测图计划(Shuttle Radar Topography Mission-SRTM)，通过装载在奋进号航天飞机上的雷达和从航天飞机上伸出 60m长臂端的另一副雷达天线所构成的航天干涉成像雷达系统的连续工作，完成了对地球北纬 60?至南纬 56?之间陆地表面的三维地形测绘，生成了全球 80%的陆地表面的数字高程模型(DEM)。目前 NIMA 发布的 SRTM DEM 在欧亚板块的经纬度方向分辨率为 3″(约 90m)，其高程方向的精度约为 10 米，在地形起伏较大的地区，精度会差一些，用户可以免费 下载 课程表模板下载资产负债表下载英语单词下载学习机资料下载励志文章下载 使用。 NIMA发布的SRTM DEM每个数据文件文件覆盖范围为 1?×1?，将地面划分为1200×1200 个小区域，数据格式为RAW格式的，文件的命名依据DEM左下角的点所在的经纬度坐标编号，如N30E120.hgt，表示该DEM的左下角的点的经纬度坐标为北纬 30?、东经 120?。根据DEM的文件名以及研究区域的经纬度范围，可以很快地确定所需要的DEM数据。 SRTM DEM 在 InSAR 干涉测量中，将起到十分重要的作用，特别是在地 形起伏不大的地区，由于精度较高，一方面，可用于干涉基线的精确估计，另一方面，在差分干涉测量中可以消除的地形相位的影响。 4.2.2 ASAR数据介绍 ENVISAT为欧洲太空总署(EAS)为延续ERS-1/2的地球观测任务，于2002 年3月所发射的卫星。ENVISAT为太阳同步卫星，飞行高度约800km，重复观测周期为35天。ENVISAT上共载有10个载具，其中包括一个合成孔径雷达(SAR)系统，名为ASAR，为多极化雷达，可提供HH/VV，HH/HV，HH/VH资料。下表为ENVISAT卫星参数表，详细资料请参考 本文采用的是欧空局发射的极轨对地观测卫星ENVISAT所载的ASAR(Advanced Synthetic Aperture Radar)传感器在交叉极化模式(Alternating Polarization Mode)下生成的SAR影像。ENVISAT-1卫星是欧洲迄今建造的最大的环境卫星，也是费用最高的地球观测卫星(总研制成本约25亿美元)。星上载有10种探测设备，其中4种是ERS-1/2所载设备的改进型，所载最大设备是先进的合成孔径雷达(ASAR)，可生成海洋、海岸、极地冰冠和陆地的高质量图象，为科学家提供更高分辨率的图象来研究海洋的变化。其他设备将提供更高精度的数据，用于研究地球大气层及大气密度(作为ERS-1/2合成孔径雷达卫星的延续，ENVISAT-1数据主要用于监视环境，即对地球表面和大气层进行连续的观测，供制图、资源勘查、气象及灾害判断之用。 ASAR测量模式简介 ENVISAT上搭载的ASAR(Advanced Synthetic Aperture Radar)运行波段为C波段，ASAR是一种高分辨率的侧视成像雷达，一共有五种不同的测量模式，分别为图像模式(Image Mode)、交叉极化模式(Alternating Polarization Mode)、宽带模式(Wide Swath Mode)、全球监测模式(Global Monitoring Mode)以及波模式(Wave Mode)(各种模式的工作特性见下表。 成像宽度 最大 100km 最大100km 约400km 约400km 5km 下行速率 100Mbit/s 极化方式 VV或 HH 分辨率 30m 100Mbit/s VV/HH或VV/VH 或HH/HV 30m 100Mbit/s VV 或HH 0.9Mbit/s VV或HH 0.9Mbit/s VV或HH 150m 1000m 10m 其中交叉极化测量模式每次获得2幅配准的图像，带宽为7个可选择带宽之一，两幅配准图像对的极化方式为HH/VV、HH/HV或者VV/VH之一，空间分辨率大约为30m。 ASAR数据特点简介 ASAR工作在C波段，可为每个轨道连续地获取30min图像。它继承了ERS-1/2AMI中的成像模式和波模式，增强了在覆盖率、入射角范围、极化和工作模式上的功能(其主要优点表现如下: ?用ScanSAR可以达到500km的幅照宽度; ?可以获得垂直和水平极化的信息; ?交替极化模式可使目标同时以垂直极化和水平极化方式成像 ?有不同的空间分辨率和数据率; ?可提供7个条带，入射角在15?,45?的雷达数据。 SLC数据是基于CEOS(Cornmiteefor Earth ObservationSatellite)的标准格式提供给用户的。单视复数图像(Single Look Complex，SLC)，是用最新的各种卫星参数(辅助数据)加以修正之后的原始SAR图像，SLC数据除了提供复 影像数据外，随影像数据一起提供的辅助数据(auxiliary)中还包含了卫星传感器的参数、遥感平台在遥感过程中的位置矢量、相应的时间参数、影像四角和中心的地理坐标等。这种图像是采用斜距投影的(slant range)。 4.3干涉数据的选取 在确定研究目的和研究区域后，进行差分干涉测量的第一步是选择合适的干 涉SAR数据，干涉SAR数据的选择至关重要，它的选择正确与否不仅关系到干涉测量的精度问题，还影响到干涉图能否生成，如果选取的数据不合适，则很有可能形成不了干涉。 关于 InSAR 干涉数据的选取，需考虑以下几个方面: 1、传感器的类型传感器的类型以及平台的特性 传感器的类型包括一些主要的参数如波长、发射器的带宽、信噪比、轨道倾角以及重复周期等。星载SAR重复轨道干涉测量所发射的波长一般都从L波段(λ=23cm)到X波段(λ=3cm)。若波长大于L波段时，电离层会严重影响观测效果，波长小于X波段时，雷达信号会对天气状况非常敏感。波长同时也影响着干涉条纹图的密度。另外，微波对各种地物的穿透深度还与波长有很大关系，波长越长，穿透能力越强。平台的倾角和重复周期决定了地面的覆盖情况、两极盲区的出现、以及获取SAR图像的重访时间。最后，精确的跟踪装置和轨道维护程序的可用性影响了干涉基线的准确度。现有的星载干涉SAR传感器主要有ERS-1/2、JERS、Radarsat、ENVISAT ASAR、ALOS PALSAR等，它们都有各自的特性。 2、在确定传感器和平台的类型后需要考虑数据的实用性 主要是考虑存档数据(或历史数据)是否可用，由于数据量很大，即使卫星飞过某些地区，地面站也可能没有接收雷达信号，致使该地区没有可用的存档数据。在进行干涉测量时，不可能由不同的航迹号(Track)获取的数据生成干涉图像对。另外需要考虑卫星飞行的方向，由于 SAR 是侧视成像，在同一区域获取的升轨(Ascending-卫星由南向北飞行)与降轨(Descending-卫星由北向南)的图像信息差别很大，用于干涉的图像只能全是升轨，或者全是降轨的。在我国，降轨数据指的是雷达卫星由北向南飞行过程中采集数据，而升轨数据指的是雷达卫星由南向北飞行过程中采集数据，需要注意的是，降轨数据在处理完后需要对其进行左右翻转，而升轨数据在处理完后需要对其进行上下翻转。另外，在我国，降轨数据多是在夜间获取的，而升轨数据是在白天获取的。 3、空间基线的影响 空间基线(spatial baseline):即重复轨道的两个轨道之间的距离，空间基线应该小于临界基线并且大于最小基线，应该取一个最优长度的基线(即最优基线)。 空间基线太长，会增大图像间的去相关性，而太小又难以形成干涉。对InSAR 有着重要影响的应该是垂直基线，一方面它决定了空间去相关的大小，另一方面又影响着生成的干涉图的条纹密度。下面以ERS 卫星为例进行说明，对于地形测图为目的的干涉测量，要产生高质量的 DEM，其垂直基线的选取一般在 100 米到 500米，根据经验 700 米是容许上限。而对于变形测量为目的的差分干涉测量，应该选取较小的垂直基线，垂直基线距最好在 50 米到 100 米之间，这样可以减小地形因素、噪声等的 影响。 4、时间基线的影响 时间基线(temporal baseline)引起的去相关(decorrelation)现象对干涉效果的影响较重。因为，重复轨道获取的干涉图像对通常要间隔几天或几十天，在这期间，地表状况的微小变化，如:湿度、植被、风速等，均可能会附加上干扰相位。时间去相关的解决很大程度上取决传感器的设计和运行方式，目前 ERS-1 和 ERS-2 的双星串联飞行模式可以将时间基线缩短至 1 天，在一定程度上解决了时间去相关的影响。对于地形测图为目的的干涉测量，其时间基线的选取要尽可能的短来减小时间去相关的影响。而对于变形测量为目的的差分干涉测量，理想的时间基线取决于具体的应用，也就是需要考虑变形的性质，是连续变形还是瞬时变形。对于瞬时变形，应该选取较短的时间基线;而对于连续变形，应该保证在变形和干涉图的其它影响(如噪声、大气效应、残余地形影响)之间有足够高的信噪比的情况下，根据预计变形量来选取时间基线。 5、成像区域地形等因素的影响 成像时地形的特性包括:粗糙度、高度范围、植被覆盖类型及覆盖情况、土壤湿度、以及人为活动的影响。另外是与季节相关的影响，主要由于体散射引起的，体散射与介质内的多次折射有关。由于一般地表都有植被覆盖，不同地区不同季节地面植被的生成情况差异很大，如农作物、矮小的植被、高大的树木等，它们的影响各不相同，因此，应该选择植被枯萎期时获取雷达影像，一般在我国北方，最好选取早春、晚秋、冬季时获取 SAR 图像，这时受地面植被的干扰最小。 6、成像时大气状况的影响 主要考虑对流层(troposphere)、电离层(ionosphere)对传播信号的影响，大气 效应等引起附加相位是选取数据时需要考虑的另一个问题。由于SAR图像是不同时间获取的，大气对两幅图像的影响就不同，特别的，如果遇上风暴等引起的涡流，可能会引起信号的中断，影响信号在雷达与地面之间的传播。另外，当气压、温度、湿度发生变化时，会影响大气层的厚度和密度，所有这些都会影响到干涉图中的相位差的变化，所以要选择气候比较稳定的时间获取图像，试验表明，选取夜间获取的图像要比白天获取的干涉效果好。 4.4 利用二轨法对实验区数据进行处理 4.4.1 二轨法技术流程与步骤 步骤: 4.4.2 二轨法处理实例 4.5 利用SBAS方法对实验区进行处理 4.5.1 SBAS流程与步骤 步骤:1: 4.5.2 SBAS实例 第五章结论与展望 参考文献 结束语