• 四川郎酒股份有限公司获第十二届人民企业社会责任奖年度环保奖 2019-05-13
  • 银保监会新规剑指大企业多头融资和过度融资 2019-05-12
  • 韩国再提4国联合申办世界杯 中国网友无视:我们自己来 2019-05-11
  • 中国人为什么一定要买房? 2019-05-11
  • 十九大精神进校园:风正扬帆当有为 勇做时代弄潮儿 2019-05-10
  • 粽叶飘香幸福邻里——廊坊市举办“我们的节日·端午”主题活动 2019-05-09
  • 太原设禁鸣路段 设备在测试中 2019-05-09
  • 拜耳医药保健有限公司获第十二届人民企业社会责任奖年度企业奖 2019-05-08
  • “港独”没出路!“梁天琦们”该醒醒了 2019-05-07
  • 陈卫平:中国文化内涵包含三方面 文化复兴表现在其中 2019-05-06
  • 人民日报客户端辟谣:“合成军装照”产品请放心使用 2019-05-05
  • 【十九大·理论新视野】为什么要“建设现代化经济体系”?   2019-05-04
  • 聚焦2017年乌鲁木齐市老城区改造提升工程 2019-05-04
  • 【专家谈】上合组织——构建区域命运共同体的有力实践者 2019-05-03
  • 【华商侃车NO.192】 亲!楼市火爆,别忘了买车位啊! 2019-05-03
    • / 27
    • 下载费用:30 金币  

    重庆时时彩官方购买平台: 一种适用于大梯度地表沉降监测的时序差分雷达干涉方法.pdf

    关 键 词:
    一种 适用于 梯度 地表 沉降 监测 时序 雷达 干涉 方法
      专利查询网所有资源均是用户自行上传分享,仅供网友学习交流,未经上传用户书面授权,请勿作他用。
    摘要
    申请专利号:

    CN201710019026.4

    申请日:

    2017.01.11

    公开号:

    CN106772342A

    公开日:

    2017.05.31

    当前法律状态:

    实审

    有效性:

    审中

    法律详情: 实质审查的生效IPC(主分类):G01S 13/08申请日:20170111|||公开
    IPC分类号: G01S13/08 主分类号: G01S13/08
    申请人: 西南石油大学
    发明人: 于冰; 马德英; 肖东升; 杨莹辉; 刘福臻; 贾宏亮; 苏勇; 王继燕
    地址: 610500 四川省成都市新都区新都大道8号
    优先权:
    专利代理机构: 四川君士达律师事务所 51216 代理人: 芶忠义
    PDF完整版下载: PDF下载
    法律状态
    申请(专利)号:

    CN201710019026.4

    授权公告号:

    |||

    法律状态公告日:

    2017.06.23|||2017.05.31

    法律状态类型:

    实质审查的生效|||公开

    摘要

    本申请公开了一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,该方法通过短时间基线差分干涉图筛选、离散点相位解缠、基于短时间基线差分干涉图的形变分量建模和解算、形变分量可靠性检验、差分干涉图相位梯度修正、对上述过程进行迭代以确保形变分量解算正确,以及基于修正后差分干涉图的形变时间序列建模和解算等过程,并形成整体的技术方案,解决原有时序差分雷达干涉技术中由于形变和相位梯度较大导致形变相位模糊度解算困难或失败的问题,最终达到正确提取大梯度地表形变速率和形变时间序列的目的,并起到降低大梯度形变建模和解算所需合成孔径雷达影像数量的效果,节约时序差分雷达干涉的应用经济成本。

    权利要求书

    1.一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,其特征在于,具体按照
    以下步骤实施:
    步骤1,对所有筛选后的SAR影像进行任意干涉组合并计算时间基线和空间基线;
    步骤2,设置时空基线阈值进行干涉对初始筛选,在保证干涉对数量的前提下限制时间
    失相干和空间失相干,按照干涉组合进行差分干涉处理得到差分干涉图集;
    步骤3,永久散射体PS和相干目标CT探测及点集合并,得到相干散射体CS点集,提取CS
    点集上的差分干涉相位,并以CS为节点构建不规则Delaunay三角网TIN;
    步骤4,基于Delaunay三角网和MCF方法的CS相位解缠;
    步骤5,线性形变速率和高程误差建模及解算;
    步骤6,基于数据模拟或地面测量数据的线性形变解算结果检验,若结果不可靠则执行
    步骤7,否则转向步骤8;
    步骤7,重新筛选参与计算的差分干涉图,并重新执行步骤5和步骤6,当达到计算终止
    条件时,执行步骤8;
    步骤8,从所有原始差分干涉图中减去线性形变相位分量,对残差相位进行重新解缠,
    然后再将线性形变相位分量加回重新解缠后的相位中,重新执行步骤5,得到更新后的线性
    形变速率和高程误差;
    步骤9,记录步骤8中所得到的线性形变速率,并从所有的原始差分干涉图中减去更新
    后的线性形变和高程误差相位分量,重新执行基于离散点的相位解缠,然后将步骤8中的线
    性形变分量重新加回新的解缠相位中,执行形变时间序列建模和解算过程,最终得到形变
    时间序列。
    2.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,其
    特征在于,步骤1具体实施方式为:排除受雨雪天气以及积雪影响的SAR影像后进行任意干
    涉组合配对,假设有N+1幅SAR影像,通过任意干涉组合形成N(N+1)/2个干涉对,每个干涉对
    由主、从两幅SAR影像构成,然后,根据每个干涉对主、从SAR影像的获取时间计算该干涉对
    的时间基线即SAR影像获取的时间差;根据主、从SAR影像的参数文件所记录的SAR传感器运
    行位置及相应的时间参数计算该干涉对的空间基线,主、从影像成像时SAR传感器的空间距
    离,取该空间距离在垂直于SAR传感器视线方向上的分量为空间基线,也即垂直基线。
    3.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,其
    特征在于,步骤2中在进行时序差分干涉处理时,为降低时间失相干和空间失相干的影响,
    采用时间基线阈值和空间基线阈值方法排除时间基线和空间基线大于某给定的阈值的干
    涉对,当SAR影像获取时间整体跨度小于2年时不考虑对时间基线进行限制。
    4.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,其
    特征在于,步骤3具体实施方式为:PS探测采用振幅阈值和振幅离差指数ADI阈值双阈值方
    法,对于SAR影像中某一特定的像元,其在N+1幅SAR影像中的振幅值直接从SAR振幅影像中
    提取,则其时序振幅均值和ADI值为:
    <mrow> <mover> <mi>a</mi> <mo>&OverBar;</mo> </mover> <mo>=</mo> <mfrac> <mrow> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mo>+</mo> <mn>1</mn> </mrow> </munderover> <msub> <mi>a</mi> <mi>i</mi> </msub> </mrow> <mrow> <mi>N</mi> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo> </mrow> </mrow>
    <mrow> <msub> <mi>D</mi> <mi>a</mi> </msub> <mo>=</mo> <mfrac> <msub> <mi>&sigma;</mi> <mi>a</mi> </msub> <mover> <mi>a</mi> <mo>&OverBar;</mo> </mover> </mfrac> <mo>,</mo> <msub> <mi>&sigma;</mi> <mi>a</mi> </msub> <mo>=</mo> <msqrt> <mrow> <mfrac> <mn>1</mn> <mrow> <mi>N</mi> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mo>&Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mo>+</mo> <mn>1</mn> </mrow> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>a</mi> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>a</mi> <mo>&OverBar;</mo> </mover> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo> </mrow> </mrow>
    其中,i为SAR影像索引号(序号);ai和分别为该像元在第i幅SAR影像中的振幅值及在
    所有影像中的振幅值的均值;Da和σa分别为该像元的ADI值及时序振幅标准差。当该像元的
    振幅信息满足如下条件时,认为其为PS:
    <mrow> <mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mrow> <msub> <mi>D</mi> <mi>a</mi> </msub> <mo>&le;</mo> <mn>0.25</mn> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mover> <mi>a</mi> <mo>&OverBar;</mo> </mover> <mo>&GreaterEqual;</mo> <mover> <mi>A</mi> <mo>&OverBar;</mo> </mover> <mo>&PlusMinus;</mo> <mn>2</mn> <msub> <mi>&sigma;</mi> <mi>A</mi> </msub> <mo>,</mo> <mover> <mi>A</mi> <mo>&OverBar;</mo> </mover> <mo>=</mo> <mfrac> <mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mo>+</mo> <mn>1</mn> </mrow> </munderover> <msub> <mi>A</mi> <mi>i</mi> </msub> </mrow> <mrow> <mi>N</mi> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <mo>,</mo> <msub> <mi>A</mi> <mi>i</mi> </msub> <mo>=</mo> <mfrac> <mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>c</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>C</mi> </munderover> <munderover> <mi>&Sigma;</mi> <mrow> <mi>l</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>L</mi> </munderover> <msub> <mi>a</mi> <mrow> <mi>c</mi> <mi>l</mi> </mrow> </msub> </mrow> <mrow> <mi>C</mi> <mi>L</mi> </mrow> </mfrac> <mo>,</mo> <msub> <mi>&sigma;</mi> <mi>A</mi> </msub> <mo>=</mo> <msqrt> <mrow> <mfrac> <mn>1</mn> <mrow> <mi>N</mi> <mo>+</mo> <mn>1</mn> </mrow> </mfrac> <munderover> <mi>&Sigma;</mi> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mrow> <mi>N</mi> <mo>+</mo> <mn>1</mn> </mrow> </munderover> <msup> <mrow> <mo>(</mo> <mrow> <msub> <mi>A</mi> <mi>i</mi> </msub> <mo>-</mo> <mover> <mi>A</mi> <mo>&OverBar;</mo> </mover> </mrow> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo> </mrow> </mrow>
    其中,和σA分别为N+1幅SAR影像所有像元时序振幅值的均值和标准差;c和l分别为像
    元的列号和行号,C和L分别为影像列数和行数;Ai和acl分别为第i幅SAR影像所有像元振幅
    值的均值和行、列号分别为c和l的某像元的振幅值,公式(3)中0.25即为ADI阈值,
    即为振幅阈值;
    CT探测采用相干系数阈值方法,假设由N+1幅SAR影像形成了L个干涉对,并通过差分干
    涉数据处理得到L幅差分干涉图,通过下式计算每个像元在所有干涉图中相应的相干系数:
    <mrow> <msub> <mover> <mi>&gamma;</mi> <mo>^</mo> </mover> <mi>l</mi> </msub> <mo>=</mo> <mfrac> <mrow> <mo>|</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>z</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>Z</mi> </munderover> <msubsup> <mi>i</mi> <mrow> <msub> <mi>IM</mi> <mi>l</mi> </msub> </mrow> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> </msubsup> <msup> <msubsup> <mi>i</mi> <mrow> <msub> <mi>IS</mi> <mi>l</mi> </msub> </mrow> <mo>*</mo> </msubsup> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> </msup> <mo>|</mo> </mrow> <msqrt> <mrow> <munderover> <mi>&Sigma;</mi> <mrow> <mi>z</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>Z</mi> </munderover> <mo>|</mo> <msubsup> <mi>i</mi> <mrow> <msub> <mi>IM</mi> <mi>l</mi> </msub> </mrow> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> </msubsup> <msup> <mo>|</mo> <mn>2</mn> </msup> <mo>+</mo> <munderover> <mi>&Sigma;</mi> <mrow> <mi>z</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>Z</mi> </munderover> <mo>|</mo> <msubsup> <mi>i</mi> <mrow> <msub> <mi>IS</mi> <mi>l</mi> </msub> </mrow> <mrow> <mo>(</mo> <mi>z</mi> <mo>)</mo> </mrow> </msubsup> <msup> <mo>|</mo> <mn>2</mn> </msup> </mrow> </msqrt> </mfrac> <mo>,</mo> <mi>l</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mi>L</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo> </mrow> </mrow>
    其中,为某像元在第l幅干涉图中的相干系数值;IMl和ISl分别为第l个干涉对的主影
    像和从影像;Z为相干系数估计窗口内像元数目,z为窗口内像元索引;l为干涉图索引;当某
    个像元满足以下条件时被识别为CT:
    <mrow> <mi>m</mi> <mi>i</mi> <mi>n</mi> <mrow> <mo>(</mo> <msub> <mover> <mi>&gamma;</mi> <mo>^</mo> </mover> <mi>l</mi> </msub> <mo>)</mo> </mrow> <mo>&gt;</mo> <msub> <mi>&gamma;</mi> <mrow> <mi>c</mi> <mi>r</mi> <mi>i</mi> <mi>t</mi> </mrow> </msub> <mo>,</mo> <mi>l</mi> <mo>=</mo> <mn>1</mn> <mo>,</mo> <mn>2</mn> <mo>,</mo> <mo>...</mo> <mo>,</mo> <mi>L</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo> </mrow> </mrow>
    其中,min(·)表示取变量的最小值;γcrit为判别某个像元是否为CT的相干系数阈值,
    取0.25至0.3;
    当探测出所有的PS和CT后,将二者进行合并,并去除重复点,得到CS点集,最后从所有
    差分干涉图中提取CS点上的相位值。
    5.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,其
    特征在于,步骤4中相位解缠采用基于离散点的相位解缠方法对CS上的差分干涉相位进行
    解缠处理;对于某个CS点P,其在差分干涉图中的相位和该点上的真实相位的关系为:

    其中,为P点上的真实相位,即相位解缠要得到的相位值;φP为P点上的缠绕相位值即
    P在差分干涉图中对应的相位值,位于[-π,π)区间内,只记录了不足整周(2π)的小数部分,
    存在整周模糊度;nP为整周模糊数;
    首先,根据Delaunay剖分法则,以所有的CS为顶点构造Delaunay三角网,然后以网络边
    端点对应的两个CS之间的缠绕相位差为观测量,估算两点间的绝对相位差,通过最小费用
    流MCF方法对与相位不连续性相关的网络费用流进行估算,并寻找最小费用流对应的积分
    路径,进行相位积分,完成相位解缠,此过程也即求解公式(6)中nP的过程。
    6.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,其
    特征在于,步骤5中在对CS上的差分干涉相位进行相位解缠后,根据步骤4中构建的
    Delaunay三角网重新计算相邻CS点间的解缠相位差,此处,相邻的CS点是指Delaunay三角
    网中三角形边的端点,以任意一边的两端点P和Q为例,二者在第l个差分干涉图中对应于的
    解缠相位为:


    其中,和分别为第l幅差分干涉图的时间基线和垂直基线;λ为雷达波波长;θP和θQ
    分别为P和Q点处的雷达波入射角;RP和RQ分别为雷达天线到P和Q之间的斜距;和
    分别为P和Q点上解缠后的差分干涉相位;vP和vQ分别为P和Q点的线性形变速率;δhP和
    δhQ分别为P和Q点的高程误差;和分别为P和Q点在第l幅差分干涉图中的非
    线性形变相位;和分别为P和Q点在第l幅差分干涉图中的轨道误差相位;
    和分别为P和Q点在第l幅差分干涉图中的大气延迟相位;和分别为
    P和Q点在第l幅差分干涉图中的噪声相位;
    线性形变和高程误差建模及解算的目的是对vP和vQ及δhP和δhQ进行估算,以P和Q上的解
    缠相位为观测量进行网络邻域差分建模,对于P和Q,二者之间的网络邻域差分相位增量为:

    其中,和分别为P和Q点处斜距的平均值和入射角的平均值;为P和Q之间的
    邻域差分相位增量;ΔvPQ为P和Q之间的线性形变速率增量;ΔδhPQ为P和Q之间的高程误差
    增量;为第l幅差分干涉图中P和Q之间的相位残差增量,即P和Q点上非线性形变、轨
    道误差、大气延迟和噪声相位之间的差值之和;
    以L幅差分干涉图为例,对于P和Q所对应的网络边而言,列出L个与式(9)相同的方程,
    组成相应的线性方程组,将其表达为矩阵的形式有:
    ΔΨ=AX+W (10)
    其中,

    <mrow> <mi>X</mi> <mo>=</mo> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <msub> <mi>&Delta;v</mi> <mrow> <mi>P</mi> <mi>Q</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <msub> <mi>&Delta;&delta;h</mi> <mrow> <mi>P</mi> <mi>Q</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>12</mn> <mo>)</mo> </mrow> </mrow>

    A=[κ,η] (14)
    其中,
    <mrow> <mi>&kappa;</mi> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mn>4</mn> <mi>&pi;</mi> </mrow> <mi>&lambda;</mi> </mfrac> <msup> <mrow> <mo>&lsqb;</mo> <msubsup> <mi>B</mi> <mn>1</mn> <mi>T</mi> </msubsup> <mo>,</mo> <msubsup> <mi>B</mi> <mn>2</mn> <mi>T</mi> </msubsup> <mo>,</mo> <mo>...</mo> <mo>,</mo> <msubsup> <mi>B</mi> <mi>l</mi> <mi>T</mi> </msubsup> <mo>,</mo> <mo>...</mo> <mo>,</mo> <msubsup> <mi>B</mi> <mi>L</mi> <mi>T</mi> </msubsup> <mo>&rsqb;</mo> </mrow> <mi>T</mi> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>15</mn> <mo>)</mo> </mrow> </mrow>
    <mrow> <mi>&eta;</mi> <mo>=</mo> <mo>-</mo> <mfrac> <mrow> <mn>4</mn> <mi>&pi;</mi> </mrow> <mi>&lambda;</mi> </mfrac> <msup> <mrow> <mo>&lsqb;</mo> <mrow> <mfrac> <msubsup> <mi>B</mi> <mn>1</mn> <mo>&perp;</mo> </msubsup> <mrow> <mover> <mi>R</mi> <mo>&OverBar;</mo> </mover> <mi>sin</mi> <mover> <mi>&theta;</mi> <mo>&OverBar;</mo> </mover> </mrow> </mfrac> <mo>,</mo> <mfrac> <msubsup> <mi>B</mi> <mn>2</mn> <mo>&perp;</mo> </msubsup> <mrow> <mover> <mi>R</mi> <mo>&OverBar;</mo> </mover> <mi>sin</mi> <mover> <mi>&theta;</mi> <mo>&OverBar;</mo> </mover> </mrow> </mfrac> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mfrac> <msubsup> <mi>B</mi> <mi>l</mi> <mo>&perp;</mo> </msubsup> <mrow> <mover> <mi>R</mi> <mo>&OverBar;</mo> </mover> <mi>sin</mi> <mover> <mi>&theta;</mi> <mo>&OverBar;</mo> </mover> </mrow> </mfrac> <mo>,</mo> <mn>...</mn> <mo>,</mo> <mfrac> <msubsup> <mi>B</mi> <mi>L</mi> <mo>&perp;</mo> </msubsup> <mrow> <mover> <mi>R</mi> <mo>&OverBar;</mo> </mover> <mi>sin</mi> <mover> <mi>&theta;</mi> <mo>&OverBar;</mo> </mover> </mrow> </mfrac> </mrow> <mo>&rsqb;</mo> </mrow> <mi>T</mi> </msup> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>16</mn> <mo>)</mo> </mrow> </mrow>
    此时,方程组中仅有ΔvPQ和ΔδhPQ两个未知数,观测值数量为L,通过最小二乘可得未知
    参数的估值为:
    <mrow> <mfenced open = "[" close = "]"> <mtable> <mtr> <mtd> <mrow> <mi>&Delta;</mi> <msub> <mover> <mi>v</mi> <mo>^</mo> </mover> <mrow> <mi>P</mi> <mi>Q</mi> </mrow> </msub> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>&Delta;</mi> <mi>&delta;</mi> <msub> <mover> <mi>h</mi> <mo>^</mo> </mover> <mrow> <mi>P</mi> <mi>Q</mi> </mrow> </msub> </mrow> </mtd> </mtr> </mtable> </mfenced> <mo>=</mo> <msup> <mrow> <mo>(</mo> <msup> <mi>A</mi> <mi>T</mi> </msup> <mi>A</mi> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msup> <mi>A</mi> <mi>T</mi> </msup> <mi>&Delta;</mi> <mi>&Psi;</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>17</mn> <mo>)</mo> </mrow> </mrow>
    对于所有通过网络边连接的CS点,均通过该过程解算相邻两点间的线性形变速率增量
    和高程误差增量,在得到所有CS点间的线性形变速率增量和高程误差增量后,以此作为观
    测量,以每个CS点的线性形变速率和高程误差为待估参数,在空间域进行最小二乘平差计
    算,解算所有CS点上的线性形变速率和高程误差,线性形变速率为:
    <mrow> <msub> <mover> <mi>v</mi> <mo>^</mo> </mover> <mi>P</mi> </msub> <mo>-</mo> <msub> <mover> <mi>v</mi> <mo>^</mo> </mover> <mi>Q</mi> </msub> <mo>=</mo> <mi>&Delta;</mi> <msub> <mover> <mi>v</mi> <mo>^</mo> </mover> <mrow> <mi>P</mi> <mi>Q</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>r</mi> <mrow> <mi>P</mi> <mi>Q</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>18</mn> <mo>)</mo> </mrow> </mrow>
    其中,和分别为点P和Q的线性形变速率估计值;为点P和Q之间的线性形变速率
    增量估值;rPQ为最小二乘残差,对于所有的CS点(设数量为G)和网络边(设数量为H)而言,根
    据公式(18)得到相应的观测方程,表达为矩阵的形式为:
    <mrow> <mi>B</mi> <mover> <mi>V</mi> <mo>^</mo> </mover> <mo>=</mo> <mi>&Delta;</mi> <mover> <mi>V</mi> <mo>^</mo> </mover> <mo>+</mo> <mi>r</mi> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>19</mn> <mo>)</mo> </mrow> </mrow>
    其中,B为系数矩阵,其元素由弧段所对应的方程式确定,由1、-1和0组成;矩阵中为待
    估参数,即每个CS点的线性形变速率;中为线性形变速率增量的估值;r为残差;由最小
    二乘平差得:
    <mrow> <mover> <mi>V</mi> <mo>^</mo> </mover> <mo>=</mo> <msup> <mrow> <mo>(</mo> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>B</mi> <mo>)</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msup> <mi>B</mi> <mi>T</mi> </msup> <mi>P</mi> <mi>&Delta;</mi> <mover> <mi>V</mi> <mo>^</mo> </mover> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>20</mn> <mo>)</mo> </mrow> </mrow>
    其中,P为线性速率增量的权阵,通过弧段最小二乘残差进行估算得到。
    7.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,其
    特征在于,步骤6中采用基于差分干涉图模拟的验证方法;
    使用估计所得到的年形变速率或累积形变量和相应的干涉基线参数模拟差分干涉图,
    即将形变量转换为形变相位,形变速率到相位的转换公式为:

    其中,为线性形变在第l幅差分干涉图中对应的相位;v为通过步骤5中的数据处理
    所得的线性形变速率;
    在大形变梯度区域,差分干涉图中的形变相位梯度随时间基线的增大而显著提升,在
    差分干涉图中表现为密集的条纹变化,若估计所得形变量符合或接近真实情况,使用形变
    结果模拟的差分干涉图与原始差分干涉图相似,则执行步骤8;若估计所得形变量不符合真
    实情况,则模拟所得条纹数目与原始差分干涉图会存在明显区别,则执行步骤7。
    8.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,其
    特征在于,步骤8具体为:
    在得到正确的线性形变后,从所有的原始差分干涉图中减去线性形变所对应的相位分
    量,设此时所得到的正确的线性形变为仍以第l幅差分干涉图为例,其所对应的相位分
    量为:

    此时,通过下式获取扣除线性形变分量后的相位:

    其中,为扣除线性形变分量后的缠绕相位;φdiff,l为原始的差分干涉相位;conj
    (·)为求复数共轭;对所有L幅差分干涉图执行此步运算,得L幅修正(减小)相位梯度后的
    差分干涉图;此时,差分干涉图中的相位梯度将显著降低,条纹数目减少,有利于相位解缠
    的正确执行;
    在得到L幅经过相位梯度修正的差分干涉图后,重新对干涉图集执行步骤4中所述的相
    位解缠过程,得到L幅相应的解缠相位图,新的解缠相位表达为:

    其中,为扣除线性形变分量后的解缠相位;为高程误差相位分量;
    为非线性形变相位分量;为轨道误差相位分量;为大气延迟相位分量;
    为噪声相位分量;在得到修正后差分干涉相位的解缠值后,将步骤7中所得线性形变分
    量与其相加,得到修正后的解缠相位:

    在得到修正的解缠相位后,以其为观测值重新执行步骤5中的解算流程,得到最终的线
    性形变速率和高程误差,设最终的线性形变速率和高程误差分别为和则二者所对应
    的相位为:


    9.如权利要求1所述的一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,其
    特征在于,步骤9具体为:
    在步骤8得到最终的线性形变速率和高程误差结果后,首先保留并输出线性形变速率
    结果,然后根据下式将线性形变和高程误差相位从原始差分干涉相位中扣除,得到新的差
    分干涉相位值:

    且有:
    <mrow> <msub> <mover> <mi>&phi;</mi> <mo>^</mo> </mover> <mrow> <mi>l</mi> <mi>i</mi> <mi>n</mi> <mi>e</mi> <mi>a</mi> <mi>r</mi> <mo>_</mo> <mi>t</mi> <mi>o</mi> <mi>p</mi> <mo>_</mo> <mi>c</mi> <mi>r</mi> <mi>c</mi> <mi>t</mi> <mo>,</mo> <mi>l</mi> </mrow> </msub> <mo>=</mo> <msub> <mi>&phi;</mi> <mrow> <mi>d</mi> <mi>e</mi> <mi>f</mi> <mo>_</mo> <mi>n</mi> <mi>ln</mi> <mi>r</mi> <mo>,</mo> <mi>l</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>&phi;</mi> <mrow> <mi>o</mi> <mi>r</mi> <mi>b</mi> <mo>_</mo> <mi>r</mi> <mi>e</mi> <mi>s</mi> <mo>,</mo> <mi>l</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>&phi;</mi> <mrow> <mi>a</mi> <mi>t</mi> <mi>m</mi> <mo>,</mo> <mi>l</mi> </mrow> </msub> <mo>+</mo> <msub> <mi>&phi;</mi> <mrow> <mi>n</mi> <mi>o</mi> <mi>i</mi> <mo>,</mo> <mi>l</mi> </mrow> </msub> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>29</mn> <mo>)</mo> </mrow> </mrow>
    在得到最终经线性形变和高程误差修正的相位后,重新执行步骤4中的相位解缠过程;
    在得到新的解缠相位后,重新将步骤8中得到的线性形变相位加回,得扣除高程误差后的解
    缠相位:

    其中,为经高程误差校正后的解缠相位;为由公式(26)计算得到的线性形
    变相位;
    对于某个CS而言,其在第l幅扣除了高程误差的解缠相位图中的相位表达为:

    其中,为形变相位,包括线性形变和非线性形变;为轨道误差相位;为大
    气延迟相位;为噪声相位。对于L个差分干涉图,经高程误差改正后的相位为:

    仍以上述N+1幅SAR影像为例,设其获取时刻为T=[t0,t1,…,tn,…,tN],由其所形成的L
    个干涉对的主(IM)、从(IS)SAR影像集分别为:
    IM=[IM1,IM2,…,IMl,…,IML];IS=[IS1,IS2,…,ISl,…,ISL] (33)
    此时,以t0时刻为参考时刻,其余任意时刻tn(n=1,2,…,N)相对于t0时刻的相位
    为待估参数,以扣除高程误差后的解缠相位
    为观测值,建立观测方程并恢复每个时刻的累积形变相位。为得到符合实际物理意义的形
    变估计参数,假设相邻两幅SAR影像获取时间间隔内相位的变化符合线性累积趋势,将对相
    位时间序列的求解转化为对相位变化速率vph的求解,此时,待估参数为:

    其中,为tn时刻相对于t0时刻的累积相位,且有式(34)实际的物理意义是相邻
    两幅SAR影像获取时间间隔内的相位变化速率,也称为分段相位速率;相应地,该模型被称
    为分段线性模型,此处的“段”指相邻两幅SAR影像之间的时间段,此时得观测方程为:

    将式(35)表达为矩阵的形式有:

    其中,B为L×N的矩阵,矩阵元素B[l,n]=tn-tn-1(ISl+1≤n≤IMl),其它元素值为零;
    首先对B进行奇异值分解SVD处理,然后估算出最小范数意义下各时间段的平均相位速
    率,通过在时间域上进行积分即可恢复与SAR影像获取时刻相对应的相位时间序列vph;B的
    奇异值分解为:
    B=USWT (37)
    其中,U为L×L阶正交矩阵,被称作B的左奇异向量,其前N列是BBT的特征向量;W是N×L
    阶正交矩阵,被称作B的右奇异向量;S是半正定L×L阶对角矩阵,其元素为相应于BBT特征
    向量的均方根,也即奇异值σi,且有,
    S=diag(σ1,…,σN-r+1,0,…,0) (38)
    其中,diag(·)代表对角矩阵,除对角线元素为σi外,其余元素值均为零;N-r+1为矩阵B
    的秩,r代表差分干涉图子集的数量,对公式(36)进行解算得:

    其中,S+=diag(1/σ1,…,1/σN-r+1,…,0,…,0);
    在解算出相邻时间段内的平均相位变化速率后,根据相位变化速率与时间之积并求和
    得到对应于SAR影像获取时刻的相位序列,即的值;然后,从相位序
    列中扣除由步骤8所得到的线性形变分量,并在空间域进行低通滤波去除噪声相位,在时间
    域进行高通滤波得到大气延迟和轨道误差相位序列;
    最后,从相位序列中扣除大气延迟和轨道误差相位序列,得到地表形变所对应的相位
    序列此时得形变时间序列为:

    其中,D为与每幅SAR影像获取时刻相对应的累积形变量所组成的向量;为与每幅
    SAR影像获取时刻相对应的累积形变相位,即:
    <mrow> <mi>D</mi> <mo>=</mo> <mo>&lsqb;</mo> <msub> <mi>d</mi> <msub> <mi>t</mi> <mn>1</mn> </msub> </msub> <mo>,</mo> <msub> <mi>d</mi> <msub> <mi>t</mi> <mn>2</mn> </msub> </msub> <mo>,</mo> <mo>...</mo> <mo>,</mo> <msub> <mi>d</mi> <msub> <mi>t</mi> <mi>N</mi> </msub> </msub> <mo>&rsqb;</mo> <mo>-</mo> <mo>-</mo> <mo>-</mo> <mrow> <mo>(</mo> <mn>41</mn> <mo>)</mo> </mrow> </mrow>

    说明书

    一种适用于大梯度地表沉降监测的时序差分雷达干涉方法

    技术领域

    本发明属于地表沉降监测技术领域,具体地说,涉及一种适用于大梯度地表沉降
    监测的时序差分雷达干涉方法。

    背景技术

    地表沉降(垂直向地表形变,为便于表述,后续采用“形变”代表“沉降”)是发生范
    围最广的地质灾害之一,具有持续时间长的特点,且多发于城市及其近郊等经济发达和人
    口聚集区,对经济发展、城市建设和人民生活均会构成持久危害。我国是世界上地表沉降灾
    害最为严重的国家之一,累积沉降大于200毫米的面积超过15万平方公里,主要集中在华北
    平原、长江三角洲和汾渭盆地等经济发达地区,且出现了严重的沉降漏斗,造成了严重的经
    济损失。对地表沉降(尤其是沉降漏斗)开展大范围精确监测,对沉降防控及避免相应的危
    害具有十分重要的现实意义。

    目前,时序差分合成孔径雷达(Synthetic Aperture Radar,SAR)干涉
    (Differential SAR Interferometry,DInSAR)技术已在地表形变监测中得到广泛应用,且
    是微波遥感、卫星大地测量以及地球物理学领域研究和应用的热点之一。时序DInSAR(Time
    Series DInSAR,TS-DInSAR)具有覆盖范围广、空间分辨率高、效率高、精度高且不易受云雾
    和阴雨天气影响的技术优势,非常适用于开展大范围地表形变监测。TS-DInSAR对缓慢累积
    性地表形变具有较好的监测效果和精度,但当形变较快和形变梯度较大时易造成解算精度
    较低或解算失败(如形变欠估计和形变漏斗区的“空洞”现象,即结果缺失)。

    发明内容

    有鉴于此,本申请针对TS-DInSAR难以满足快速或大梯度地表形变监测需求的问
    题,提供了一种适用于大梯度地表沉降监测的时序差分雷达干涉方法。

    为了解决上述技术问题,本申请公开了一种适用于大梯度地表沉降监测的时序差
    分雷达干涉方法,具体按照以下步骤实施:

    步骤1,对所有筛选后的SAR影像进行任意干涉组合并计算时间基线(形成一个干
    涉对的两幅SAR影像的获取时间差)和空间基线(两次成像时刻SAR传感器之间的空间距离,
    一般取该距离垂直于SAR视线方向的分量);

    步骤2,设置时空基线阈值进行干涉(干涉组合)对初始筛选,在保证干涉对数量的
    前提下限制时间失相干和空间失相干,按照干涉组合进行差分干涉处理得到差分干涉图
    集;

    步骤3,永久散射体(Persistent Scatterer,PS)和相干目标(Coherent Target,
    CT)探测及点集合并,得到相干散射体(Coherent Scatterer,CS)点集,提取CS点集上的差
    分干涉相位,并以CS为节点构建不规则三角网(Triangular Irregular Network,TIN),本
    发明技术方案中采用Delaunay三角网;

    步骤4,基于Delaunay三角网和最小费用流(Minimum Cost Flow,MCF)方法的CS相
    位解缠;

    步骤5,线性形变速率和高程误差建模及解算;

    步骤6,基于数据模拟或地面测量数据的线性形变解算结果检验,若结果不可靠则
    执行步骤7,否则转向步骤8;

    步骤7,重新筛选参与计算的差分干涉图,并重新执行步骤5和步骤6,当达到计算
    终止条件时,执行步骤8;

    步骤8,从所有原始差分干涉图中减去线性形变相位分量,对残差相位(从差分干
    涉图中减去线性形变相位后剩余的相位)进行重新解缠,然后再将线性形变相位分量加回
    重新解缠后的相位中,重新执行步骤5,得到更新后的线性形变速率和高程误差(这一步主
    要是让更多干涉对参与计算,提高线性形变速率和高程误差的解算精度);

    步骤9,记录步骤8中所得到的线性形变速率,并从所有的原始差分干涉图中减去
    更新后的线性形变和高程误差相位分量,重新执行基于离散点的相位解缠,然后将步骤8中
    的线性形变分量重新加回新的解缠相位中,执行形变时间序列建模和解算过程,最终得到
    形变时间序列。

    与现有技术相比,本申请可以获得包括以下技术效果:

    本发明技术方案依据TS-DInSAR时序差分干涉图中形变相位大小以及形变相位梯
    度大小与时间间隔(差分干涉图的时间基线,即形成差分干涉图的两幅SAR影像获取时间
    差)成正相关关系这一特点,提出短时间基线差分干涉图筛选、离散点相位解缠、基于短时
    间基线差分干涉图的形变分量建模和解算、形变分量可靠性检验、差分干涉图相位梯度修
    正、对上述过程进行迭代以确保形变分量解算正确,以及基于修正后差分干涉图的形变时
    间序列建模和解算这一整体的技术方案,解决原有TS-DInSAR技术中由于形变和相位梯度
    较大导致形变相位模糊度解算困难或失败的问题,最终达到正确提取大梯度地表形变速率
    和形变时间序列的目的。

    采用本发明技术方案,只要差分干涉图序列中存在至少6个短时间基线干涉对(由
    4幅SAR影像构成)即可实现线性形变分量的解算,然后开展长时间基线差分干涉图的形变
    相位梯度修正,促使更多的差分干涉图得到正确解缠,并参与形变分量的重新估算以及形
    变时间序列的建模和解算。因此,无需很高的SAR影像使用量便可实现大梯度形变的有效提
    取。

    当然,实施本申请的任一产品并不一定需要同时达到以上所述的所有技术效果。

    附图说明

    此处所说明的附图用来提供对本申请的进一步理解,构成本申请的一部分,本申
    请的示意性实施例及其说明用于解释本申请,并不构成对本申请的不当限定。在附图中:

    图1是一种适用于大梯度地表沉降监测的时序差分雷达干涉方法的实施流程图;

    图2是35幅SAR影像任意组合所形成的干涉对的时间基线和空间基线分布图;

    图3是空间基线阈值为30米时干涉对的时空基线分布图;

    图4是第一次解算所得CS点上线性形变速率结果图;

    图5是使用第一次解算所得的线性形变模拟的差分干涉图与原始差分干涉图的对
    比;其中,a为原始差分干涉图,b为使用第一次解算所得的线性形变模拟的差分干涉图;

    图6是经3次迭代后解算所得CS点上的线性形变速率结果图;

    图7是使用第三次解算所得的线性形变模拟的差分干涉图与原始差分干涉图的对
    比;其中,a为原始差分干涉图,b为使用第三次解算所得的线性形变模拟的差分干涉图;

    图8是2009年6月23日至2010年7月2日期间的累积形变量;

    图9是图8中所标示三个特征点(A、B和C)的形变时间序列。

    具体实施方式

    以下将配合附图及实施例来详细说明本申请的实施方式,藉此对本申请如何应用
    技术手段来解决技术问题并达成技术功效的实现过程能充分理解并据以实施。

    一种适用于大梯度地表沉降监测的时序差分雷达干涉方法,如附图1所示,具体按
    照以下步骤实施:

    步骤1,对所有筛选后的SAR影像进行任意干涉组合并计算时间基线和空间基线;

    在筛选出合适的SAR影像(排除受雨雪等天气以及积雪影响的SAR影像)后,进行任
    意干涉组合配对,假设有N+1幅SAR影像,通过任意干涉组合可形成N(N+1)/2个干涉对。每个
    干涉对由主、从两幅SAR影像构成。然后,根据每个干涉对主、从SAR影像的获取时间计算该
    干涉对的时间基线(即SAR影像获取的时间差),根据主、从SAR影像的参数文件所记录的SAR
    传感器运行位置及相应的时间参数计算该干涉对的空间基线(主、从影像成像时SAR传感器
    的空间距离,实际中一般取该空间距离在垂直于SAR传感器视线方向上的分量为空间基线,
    也即垂直基线)。

    本实施例中采用覆盖天津市精武镇的35幅SAR影像为数据进行展示。附图2所示为
    35幅SAR影像进行任意组合所形成的干涉对的时间基线和空间基线分布。图2中以“年年年
    年月月日日”的方式标注了每幅SAR影像的获取时间,如“20090327”代表2009年3月27日。通
    过虚线相连的两幅SAR影像为一个干涉对,纵轴代表干涉对的空间基线大小,横轴代表时
    间,时间的差值即为时间基线。

    步骤2,设置时空基线阈值进行干涉对的初始筛选,在保证干涉对数量的前提下限
    制时间失相干和空间失相干,按照干涉组合进行差分干涉处理得到差分干涉图集;

    在进行时序差分干涉处理时,为降低时间失相干和空间失相干的影响,常采用时
    间基线阈值和空间基线阈值方法排除时间基线和空间基线大于某给定阈值的干涉对(受时
    间失相干和空间失相干影响较显著的干涉对)。大量研究表明,时间基线不宜超过2年,空间
    基线不宜超过150米。但是,当SAR影像获取时间整体跨度较小时(如第一幅SAR影像和最后
    一幅SAR影像的获取时间间隔小于2年)可不考虑对时间基线进行限制。在本实施例中,考虑
    到SAR影像获取时间的整体跨度较短(1年零9个月),故未设置时间基线阈值,仅考虑对空间
    基线进行限制,设置空间基线阈值为30米,排除空间基线大于30米的干涉对。附图3所示为
    排除空间基线大于30米的干涉对后,得到保留的干涉对的时间基线和空间基线分布?;?br />此干涉组合,进行差分干涉处理,得到相应的差分干涉图集,处理过程如下:

    以其中一个干涉对为例,设组成此干涉对的两幅SAR影像分别为I1和I2,分别记为
    主影像和从影像,对于影像中所记录的一点P(即影像像元P),其在I1和I2中所对应的同名像
    元的值分别为:

    Ψ1=A1exp(jφ1) (1)

    Ψ2=A2exp(jφ2) (2)

    其中,exp(·)代表以e为底的指数函数;Ψ1为目标P在I1中的像元值,φ1为目标P
    在I1中的相位值;Ψ2为目标P在I2中的像元值,φ2为目标P在I2中的相位值;A1为目标P在I1
    中信号的振幅值;A2为目标P在I2中信号的振幅值。其中,Ψ1和Ψ2为复数。由复数的共轭相
    乘可得P在由I1和I2形成的干涉图中的像元值为:

    Ψint,P=Ψ1Ψ2*=A1A2exp(j(φ1-φ2)) (3)

    其中,“*”代表复数的共轭,即Ψ2*为Ψ2的共轭复数。

    由式(3)可得P点所对应的干涉相位(即两次观测的相位差表达形式)为:

    φint,P=φ1-φ2 (4)

    对SAR影像中的所有像元进行相同的处理,即可得到由主、从SAR影像I1和I2所形成
    的干涉图。以P点为例,此时,P点的干涉相位可表达为:

    φint,P=φref,P+φtop,P+φdef,P+φatm,P+φadd,P (5)

    其中,φint,P为P点在由SAR影像I1和I2所形成的干涉图中的干涉相位;φref,P为P点
    的参考椭球面相位,由SAR传感器两次对同一地区观测时所处位置不同以及地球曲面引起;
    φtop,P为P点的地形相位,由地表真实高程及其变化引起;φdef,P为P点所发生的地表形变对
    应的相位;φatm,P为P点上的大气延迟相位,由雷达波在大气层中传播发生折射引起;φadd,P
    是P点上的附加相位,主要包括各种随机噪声以及失相干噪声(两次成像时地面目标发生变
    化,导致雷达回波信号产生差异所致)。

    所谓差分干涉处理,即在上述干涉处理基础上,去除地球参考椭球面相位和地形
    相位,得到以形变相位为主的差分干涉相位,形成差分干涉图。首先,在不考虑大气延迟和
    噪声相位的情况下,P点上的干涉相位可表达为:


    其中,等式右边第一项即为参考椭球面相位的表达式,第二项为地形相位的表达
    式,第三项为形变相位的表达式;B为干涉对的空间基线;θ为雷达侧视角;α为空间基线与水
    平方向的夹角;λ为雷达波长;R1为雷达传感器到观测目标之间的斜距;为空间基线B在垂
    直于雷达侧视方向上的投影分量,即垂直基线;hP为P点的高程;Δr为P点所发生的地表形
    变。需要说明的是,Bsin(θ-α)即为干涉对的平行基线(空间基线B在沿雷达侧视方向上的投
    影分量)且有:


    其中,R1和R2分别为雷达传感器两次观测成像时到地面点P的斜距;ΔR=R1-R2为
    两次成像雷达传感器到P点斜距的差。

    差分干涉处理即要从干涉相位中减去公式(6)等号右边的前两项。首先需要计算
    并扣除参考椭球面相位。由公式(6)等号右边第一项以及公式(7)可知,参考椭球面相位与
    卫星两次观测至参考椭球面间斜距的差ΔR有关,即:


    由公式(7)可知,要获得ΔR的值,需首先已知R1和R2,此二者可通过两点间距离公
    式计算:


    其中,d(·)为两点间距离方程;M和S分别为对应于主、从影像中P点成像时刻的卫
    星三维空间坐标矢量,由SAR影像参数文件提供的参数计算得到;P为P点在参考椭球面上的
    空间三维坐标矢量,可通过卫星成像参数、斜距方程、多普勒方程和椭球方程计算得到。此
    时,可通过式(8)和(9)得到P点上的参考椭球面相位。

    对于地形相位,通过外部数字高程模型(Digital Elevation Model,DEM)可获取
    地面点P的高程值hP(美国地质调查局网站公布了全球免费的DEM,本实施例中即采用该
    DEM),R1为主影像斜距,可由P点在主、从影像上的成像参数计算得到。对于垂直基线
    有:


    其中,空间基线B=d(M,S),即主、从影像成像时雷达传感器的空间距离;
    在得到hP、R1、和后,可通过式(6)等号右边第二项计算P点的地形相位。

    在得到P点的参考椭球面相位和地形相位后,从式(5)中扣除这两部分相位分量即
    可得到差分干涉相位。对SAR影像中所有的像元点进行上述差分干涉处理,即可得到相应的
    差分干涉图。

    步骤3,PS和CT探测及点集合并,得到CS点集,提取CS点集上的差分干涉相位,并以
    CS为节点构建不规则三角网(Triangular Irregular Network,TIN),本发明技术方案中采
    用Delaunay三角网;

    PS探测采用振幅阈值和振幅离差指数(Amplitude Dispersion Index,ADI)阈值
    双阈值方法。对于SAR影像中某一特定的像元,其在N+1幅SAR影像中的振幅值可直接从SAR
    振幅影像中提取,则其时序振幅均值和ADI值为:



    其中,i为SAR影像索引号(序号);ai和分别为该像元在第i幅SAR影像中的振幅值
    及在所有影像中的振幅均值;Da和σa分别为该像元的ADI值及时序振幅标准差。当该像元的
    振幅信息满足如下两个条件时,认为其为PS:


    其中,和σA分别为N+1幅SAR影像所有像元时序振幅值的均值和标准差;c和l分别
    为像元的列号和行号,C和L分别为影像列数和行数;Ai和acl分别为第i幅SAR影像所有像元
    振幅值的均值和行、列号分别为c和l的某像元的振幅值。公式(13)中0.25即为ADI阈值。

    CT探测采用相干系数阈值方法。假设由N+1幅SAR影像形成了L个干涉对,并通过差
    分干涉数据处理得到L幅差分干涉图。通过下式计算每个像元在所有干涉图中相应的相干
    系数:


    其中,为某像元在第l幅干涉图中的相干系数值;IMl和ISl分别为第l个干涉对的
    主影像和从影像;Z为相干系数估计窗口内像元数目,z为窗口内像元索引;l为干涉图索引。
    当某个像元满足以下条件时被识别为CT:


    其中,min(·)表示取变量的最小值;γcrit为判别某个像元是否为CT的相干系数
    阈值,一般可取0.25至0.3,本实施例中采用0.25。

    当探测出所有的PS和CT后,将二者进行合并,并去除重复点,得到CS点集,最后从
    所有差分干涉图中提取CS点上的差分干涉相位值。

    步骤4,基于Delaunay三角网和MCF方法的CS相位解缠;

    相位解缠一般是基于规则格网进行,常规的相位解缠方法针对差分干涉图进行处
    理,易受失相干区域的影响,而基于离散点的相位解缠方法可有效避免这一问题。本发明技
    术方案中采用基于离散点的相位解缠方法对CS上的差分干涉相位进行解缠处理。对于某个
    CS点P,其在差分干涉图中的相位和该点上的真实相位的关系为:


    其中,为P点上的真实相位(即相位解缠要得到的相位值);φP为P点上的缠绕相
    位值(即P在差分干涉图中对应的相位值),位于[-π,π)区间内,只记录了不足整周(2π)的小
    数部分,存在整周模糊度;nP为整周模糊数。

    首先,根据Delaunay剖分法则,以所有的CS为顶点构造Delaunay三角网,然后以网
    络边端点对应的两个CS之间的缠绕相位差为观测量,估算两点间的绝对相位差。通过最小
    费用流(Minimum Cost Flow,MCF)方法对与相位不连续性相关的网络费用流进行估算,并
    寻找最小费用流对应的积分路径,进行相位积分,完成相位解缠。此过程也即求解公式(16)
    中nP的过程。

    步骤5,线性形变速率和高程误差建模及解算;

    在对CS上的差分干涉相位进行相位解缠后,根据步骤4中构建的Delaunay三角网
    重新计算相邻CS点间的解缠相位差。此处,相邻的CS点是指Delaunay三角网中三角形边的
    端点。以任意一边的两端点P和Q为例,二者在第l个差分干涉图中对应于的解缠相位为:



    其中,和分别为第l幅差分干涉图的时间基线和垂直基线;λ为雷达波波长;
    θP和θQ分别为P和Q点处的雷达波入射角;RP和RQ分别为SAR传感器到P和Q之间的斜距;
    和分别为P和Q点上解缠后的差分干涉相位;vP和vQ分别为P和Q点的线性形变速率;δhP
    和δhQ分别为P和Q点上的高程误差;和分别为P和Q点在第l幅差分干涉图
    中的非线性形变相位;和分别为P和Q点在第l幅差分干涉图中的轨道误差
    相位;和分别为P和Q点在第l幅差分干涉图中的大气延迟相位;和
    分别为P和Q点在第l幅差分干涉图中的噪声相位。

    线性形变和高程误差建模及解算的目的是对vP和vQ及δhP和δhQ进行估算。在本发
    明技术方案中,以P和Q上的解缠相位为观测量进行网络邻域差分建模。对于P和Q,二者之间
    的网络邻域差分相位增量(解缠相位差)为:


    其中,和分别为P和Q点处斜距的平均值和入射角的平均值;为P和Q之
    间的邻域差分相位增量;ΔvPQ为P和Q之间的线性形变速率增量(即线性形变速率差);Δδ
    hPQ为P和Q之间的高程误差增量(即高程误差之差);为第l幅差分干涉图中P和Q之间
    的相位残差增量(即相位残差之差),即P和Q点上非线性形变、轨道误差、大气延迟和噪声相
    位之间的差值之和。

    以L幅差分干涉图为例,对于P和Q所对应的网络边而言,可列出L个与式(19)相同
    的方程,组成相应的线性方程组。将其表达为矩阵的形式有:

    ΔΨ=AX+W (20)

    其中,




    A=[κ,η] (24)

    其中,



    此时,方程组中仅有ΔvPQ和ΔδhPQ两个未知数,观测值数量为L,通过最小二乘可
    得未知参数的估值为:


    对于所有通过网络边连接的CS点,均通过该过程解算相邻两点间的线性形变速率
    增量和高程误差增量。在得到所有CS点间的线性形变速率增量和高程误差增量后,以此作
    为观测量,以每个CS点的线性形变速率和高程误差为待估参数,在空间域进行最小二乘网
    络平差计算,解算所有CS点上的线性形变速率和高程误差。对于线性形变速率有:


    其中,和分别为点P和Q的线性形变速率估计值(待估参数);为点P和Q之
    间的线性形变速率增量估值;rPQ为最小二乘残差。对于所有的CS点(设数量为G)和网络边
    (设数量为H)而言,根据公式(28)得到相应的观测方程,表达为矩阵的形式为:


    其中,B为系数矩阵,其元素由弧段所对应的方程式确定,由1、-1和0组成;矩阵
    中为待估参数,即每个CS点的线性形变速率;中为所有CS点线性形变速率增量的估值;r
    为残差矩阵。由最小二乘平差可得:


    其中,P为线性速率增量的权阵,可通过弧段最小二乘残差进行估算得到。

    如附图4所示为使用附图3中所展示的干涉对所解算得出的CS点上的线性形变速
    率。图4中以分级设色的方式表达线性形变速率的量值,LOS DR表示雷达视线向形变速率,
    mm/yr为形变速率单位,即毫米/年。在步骤6中,将对此结果进行检验。

    步骤6,基于数据模拟或地面测量数据的线性形变解算结果检验,若结果不可靠则
    执行步骤7,否则转向步骤8;

    一般地,对TS-DInSAR形变结果进行验证主要有两种途径:基于外部测量数据的验
    证方法和基于差分干涉图模拟的验证方法。由于外部数据往往不容易获得,因此本发明技
    术方案中采用基于差分干涉图模拟的验证方法。

    该方法是使用估计所得到的年形变速率或累积形变量和相应的干涉基线参数模
    拟差分干涉图。即将形变量转换为形变相位,形变到相位的转换公式为:


    其中,为线性形变在第l幅差分干涉图中对应的相位;v为通过步骤5中的数据
    处理所得的线性形变速率。

    在大形变梯度区域(如沉降漏斗区),差分干涉图中的形变相位梯度随时间基线的
    增大而显著提升,在差分干涉图中表现为密集的条纹变化,若估计所得形变量符合或接近
    真实情况,使用形变结果模拟的差分干涉图与原始差分干涉图应该相似,尤其是形变条纹
    数目应相等或差别较小。若估计所得形变量不符合真实情况(如欠估计),则模拟所得条纹
    数目与原始差分干涉图会存在明显区别。据此,可简单有效地对形变解算结果的有效性进
    行验证,且这种验证方法能够有效识别形变欠估计现象。

    如附图5所示为由步骤5解算所得线性形变模拟所得的差分干涉图与原始差分干
    涉图的对比。其中,左侧(图5a)为原始差分干涉图,右侧(图5b)为模拟的差分干涉图。二者
    条纹数目差别明显,模拟所得差分干涉图中条纹数目低于原始差分干涉图,表明发生了形
    变欠估计现象。因此,第一次解算的线性形变并不可靠。此时,需要执行步骤7。

    步骤7,重新筛选参与计算的差分干涉图,并重新执行步骤5和步骤6,当达到计算
    终止条件时,执行步骤8;

    该步骤是实现正确解算线性形变的关键。经过步骤6的检验,确定所解算结果不可
    靠时,则通过设置更为严格的时间基线重新筛选出短时间基线差分干涉图,重新执行步骤5
    和步骤6,并进行检验。具体地,重新筛选短时间基线差分干涉图是指设置更小的时间基线,
    筛选出比上一次解算所用干涉图时间基线更短的干涉图。这是由于在一般的地表沉降区,
    时间基线越短则形变越小,且形变梯度越小,差分干涉图中的形变条纹数目就越少,更容易
    对差分干涉图进行正确的相位解缠,为形变的建模和解算提供正确的观测量。例如,本实施
    例中第一次解算时未对时间基线进行限制,解算所得到的线性形变不可靠,则设置较短的
    时间基线值(本实施例中采用180天),筛选出时间基线小于180天的干涉图重新解算线性形
    变,然后再次进行结果检验,若不可靠,则进行第三次解算,本实施例中第三次解算时选取
    90天作为时间基线阈值,即筛选出时间基线小于90天的干涉图重新解算线性形变。时间基
    线逐渐缩短,即所谓的设置更为严格的时间基线限制。

    本实施例中,经过3次迭代计算可获取正确的线性形变结果,如附图6所示。附图7
    所示为利用第三次解算所得线性形变模拟的差分干涉图与原始差分干涉图的对比,其中左
    侧(图7a)为原始差分干涉图,右侧(图7b)为模拟的差分干涉图。此时,二者条纹数目几乎一
    致。得到可靠的线性形变后,可执行步骤8。

    需要指出的是,在本实施例中,通过3次迭代计算获取了正确的线性形变速率,但
    对于不同的研究区域,迭代次数是可变的,如增加或减少,以能够获取正确的线性形变为
    准。

    步骤8,从所有原始差分干涉图中减去线性形变相位分量,对残差相位(从差分干
    涉图中减去线性形变相位后剩余的相位)进行重新解缠,然后再将线性形变相位分量加回
    重新解缠后的相位中,重新执行步骤5,得到最终的线性形变速率和高程误差(这一步主要
    是让更多干涉对参与计算,提高线性形变速率和高程误差的解算精度);

    在通过第7步骤得到正确的线性形变后,从所有的原始差分干涉图中减去线性形
    变所对应的相位分量。为便于阐述,设此时所得到的正确的线性形变为仍以第l幅差分
    干涉图为例,其所对应的相位分量为:


    此时,通过下式获取扣除线性形变分量后的相位:


    其中,为扣除线性形变分量后的缠绕相位;φdiff,l为原始的差分干涉相
    位;conj(·)为求复数共轭。对所有L幅差分干涉图执行此步运算,得L幅修正时空相位梯度
    后的差分干涉图。此时,差分干涉图中的相位梯度将显著降低,条纹数目减少,有利于相位
    解缠的正确执行。

    在得到L幅经过时空相位梯度修正的差分干涉图后,重新对干涉图集执行步骤4中
    所述的相位解缠过程,得到L幅相应的解缠相位图,新的解缠相位表达为:


    其中,为扣除线性形变分量后的解缠相位;为高程误差相位分量;
    为非线性形变相位分量;为轨道误差相位分量;为大气延迟相位分量;
    为噪声相位分量。在得到修正后差分干涉相位的解缠值后,将步骤7中所得线性形变分
    量(即所对应的形变相位)与其相加,得到修正后的解缠相位:


    在得到修正的解缠相位后,以其为观测值重新执行步骤5中的解算流程,得到最终
    的线性形变速率和高程误差。为便于后续阐述,此处设最终的线性形变速率和高程误差分
    别为和则二者所对应的相位为:



    现有研究表明,参与运算的干涉对数目越多,线性形变速率和高程误差的解算精
    度越高,这也是对相位梯度进行修正使更多干涉对参与运算的重要目的之一。此外,在后续
    的数据处理中,最终的线性形变和高程误差将用于最终的相位梯度修正。

    步骤9,记录步骤8中所得到的线性形变速率,并从所有的原始差分干涉图中减去
    更新后的线性形变和高程误差相位分量,重新执行基于离散点的相位解缠,然后将步骤8中
    的线性形变分量重新加回新的解缠相位中,执行形变时间序列建模和解算过程,最终得到
    形变时间序列。

    在得到最终的线性形变速率和高程误差结果后,首先保留并输出线性形变速率结
    果,然后可根据下式将线性形变和高程误差相位(可通过公式(36)和(37)计算)从原始差分
    干涉相位中扣除,得到新的差分干涉相位值:


    其中,为扣除线性形变和高程误差相位后的差分干涉相位,且有:


    在得到最终经线性形变和高程误差修正的相位后,重新执行第4步骤中的相位解
    缠过程。在得到新的解缠相位后,重新将步骤8中得到的线性形变相位加回,得到扣除高程
    误差后的解缠相位:


    其中,为经高程误差校正后的解缠相位;为由公式(36)计算得到的线
    性形变相位。

    对于某个CS而言,其在第l幅扣除了高程误差的解缠相位图中的相位表达为:


    其中,为形变相位(包括线性形变和非线性形变);为轨道误差相位;
    为大气延迟相位;为噪声相位。对于L个差分干涉图,经高程误差改正后的相位为:


    仍以上述N+1幅SAR影像为例,设其获取时刻为T=[t0,t1,…,tn,…,tN],由其所形
    成的L个干涉对的主(IM)、从(IS)SAR影像集分别为:

    IM=[IM1,IM2,…,IMl,…,IML];IS=[IS1,IS2,…,ISl,…,ISL] (43)

    此时,以t0时刻为参考时刻,其余任意时刻tn(n=1,2,…,N)相对于t0时刻的相位
    为待估参数,以扣除高程误差后的解缠相位
    为观测值,建立观测方程并恢复每个时刻对应的相位(即相位时间序列)。为得到符合实际
    物理意义的形变估计参数,假设相邻两幅SAR影像获取时间间隔内相位的变化符合线性累
    积趋势,将对相位时间序列的求解转化为对相位变化速率(相邻时间段内的平均相位变化
    速率)vph(注意与步骤5、6、7和8中形变速率v的区别)的求解,此时,待估参数为:


    其中,为tn时刻相对于t0时刻的累积相位,且有式(44)实际的物理意义
    是相邻两幅SAR影像获取时间间隔内的相位变化速率,也称为分段相位速率。相应地,该模
    型可被称为分段线性模型。此处的“段”指相邻两幅SAR影像之间的时间段。此时可得观测方
    程为:


    将式(45)表达为矩阵的形式有:


    其中,B为L×N的矩阵,矩阵元素B[l,n]=tn-tn-1(ISl+1≤n≤IMl),其它元素值为
    零。由于干涉图集可能是不连续的,即在某个SAR影像获取时刻断开,进而造成系数
    矩阵B发生奇异(此时方程组为欠定状态,不存在唯一解)。因此,实际处理中首先对B进行奇
    异值分解(Singular Value Decomposition,SVD)处理,然后估算出最小范数意义下各时间
    段的平均相位速率,通过在时间域上进行积分即可恢复与SAR影像获取时刻相对应的相位
    时间序列vph。B的奇异值分解为:

    B=USWT (47)

    其中,U为L×L阶正交矩阵,被称作B的左奇异向量,其前N列是BBT的特征向量;W是
    N×L阶正交矩阵,被称作B的右奇异向量。S是半正定L×L阶对角矩阵,其元素为相应于BBT
    特征向量的均方根,也即奇异值σi,且有,

    S=diag(σ1,…,σN-r+1,0,…,0) (48)

    其中,diag(·)代表对角矩阵,除对角线元素为σi外,其余元素值均为零;N-r+1为
    矩阵B的秩,r代表差分干涉图子集的数量(即由于设置了时间基线和空间基线阈值限制,导
    致由N+1幅SAR影像所形成的干涉图集合被分成了多个子集,造成干涉图集不连续)。最终,
    可对公式(46)进行解算得:


    其中,S+=diag(1/σ1,…,1/σN-r+1,…,0,…,0)。

    在解算出相邻时间段内的平均相位变化速率后,根据相位变化速率与时间之积并
    求和得到对应于SAR影像获取时刻的相位序列,即的值。然后,从相
    位序列中扣除由步骤8所得到的线性形变分量,并在空间域进行低通滤波去除噪声相位,在
    时间域进行高通滤波得到大气延迟和轨道误差相位序列。最后,从相位序列中扣除大气延
    迟和轨道误差相位序列即可得到地表形变所对应的相位序列此时可得形变时间序列
    为:


    其中,D为与每幅SAR影像获取时刻相对应的累积形变量所组成的向量;为与每
    幅SAR影像获取时刻相对应的累积形变相位。即:



    如附图8所示为2009年6月23日至2010年7月2日期间的累积形变量。图8中标记了
    A、B和C三个特征点,附图9所示为三个特征点的形变时间序列。至此,便完成了大梯度形变
    区域(沉降漏斗区域)形变速率和形变时间序列的解算。

    本发明技术方案带来的有益效果:

    (1)解决原有TS-DInSAR技术难以进行大梯度地表形变建模和解算的问题

    原有TS-DInSAR技术难以适用于大梯度地表形变的建模和解算,本发明技术方案
    依据TS-DInSAR时序差分干涉图中形变相位大小以及形变相位梯度大小与时间间隔(差分
    干涉图的时间基线,即形成差分干涉图的两幅SAR影像获取时间差)成正相关关系这一特
    点,提出短时间基线差分干涉图筛选、离散点相位解缠、基于短时间基线差分干涉图的形变
    分量建模和解算、形变分量可靠性检验、差分干涉图相位梯度修正、对上述过程进行迭代以
    确保形变分量解算正确,以及基于修正后差分干涉图的形变时间序列建模和解算这一整体
    的技术方案,解决原有TS-DInSAR技术中由于形变和相位梯度较大导致形变相位模糊度解
    算困难或失败的问题,最终达到正确提取大梯度地表形变速率和形变时间序列的目的。

    针对本发明技术方案中具体的思路进行有益效果分析如下:

    (a)短时间基线差分干涉图筛选。短时间基线差分干涉图中形变相位梯度较小,较
    容易实现正确的相位模糊度解算(即相位解缠),得到正确的解缠相位后,可用于后续形变
    分量建模和解算。此外,短时间基线差分干涉图筛选是一个可以重复执行的过程。

    (b)离散点相位解缠。相位解缠也即解算相位模糊度,离散点相位解缠仅针对高质
    量的观测目标进行处理,避免低质量目标的影响。解缠后的相位将被用于形变分量建模。相
    比于原有技术基于非解缠相位进行形变分量解算而言,基于解缠相位进行形变分量建???br />避免对SAR影像时间采样率的依赖。

    (c)基于短时间基线差分干涉图的形变分量建模和解算,并检验形变分量的可靠
    性。在(a)和(b)基础上进行形变分量(线性形变速率)建模和解算(方法见技术方案的详细
    阐述部分),并采用差分干涉图模拟方法检验形变分量解算结果的可靠性。先保证形变分量
    解算正确,后续将被应用于对差分干涉图的形变相位梯度进行修正,即从差分干涉图中减
    去形变分量,起到降低形变相位梯度的作用,进而实现所有差分干涉图的正确解缠。

    (d)迭代计算。在对形变分量进行可靠性检验时,若发现结果不可靠,则重新进行
    (a)、(b)和(c),确保形变分量解算结果可靠。迭代计算是连接各个步骤的关键。

    (e)采用正确的形变分量结果重新对差分干涉图进行相位梯度修正,然后重新进
    行相位解缠,并基于此对形变时间序列进行建模和解算。

    综上所述,以上各个过程紧密结合,共同形成了本发明的整体技术方案,最终实现
    大梯度形变区域形变分量和形变时间序列的精确解算。

    (2)可降低SAR影像使用量,降低TS-DInSAR技术的应用经济成本

    在大梯度形变解算方面,现有研究指出SAR时间采样率的提高可有效缓解梯度较
    大造成的问题,但这往往受现有SAR系统观测周期(对同一地区进行重复观测的周期)以及
    研究区域可用影像数量的限制,并且数据成本会显著提高(如针对大梯度形变,以往的TS-
    DInSAR技术往往需要很高的SAR影像时间采样率,即较高的数据量需求)。尤其是高分辨率
    SAR影像,价格较高,数据成本将十分显著。

    采用本发明技术方案,理论上只要差分干涉图序列中存在至少6个短时间基线干
    涉对(由4幅SAR影像构成)即可实现线性形变分量的解算,然后开展长时间基线差分干涉图
    的形变相位梯度修正,促使更多的差分干涉图(尤其是长时间基线差分干涉图)得到正确解
    缠,并参与形变分量的重新估算以及形变时间序列的建模和解算。因此,无需很高的SAR影
    像时间采样率便可实现大梯度形变的有效提取。

    上述说明示出并描述了发明的若干优选实施例,但如前所述,应当理解发明并非
    局限于本文所披露的形式,不应看作是对其他实施例的排除,而可用于各种其他组合、修改
    和环境,并能够在本文所述发明构想范围内,通过上述教导或相关领域的技术或知识进行
    改动。而本领域人员所进行的改动和变化不脱离发明的精神和范围,则都应在发明所附权
    利要求的?;し段?。

    关于本文
    本文标题:一种适用于大梯度地表沉降监测的时序差分雷达干涉方法.pdf
    链接地址://www.4mum.com.cn/p-6001550.html
    关于我们 - 网站声明 - 网站地图 - 资源地图 - 友情链接 - 网站客服 - 联系我们

    [email protected] 2017-2018 www.4mum.com.cn网站版权所有
    经营许可证编号:粤ICP备17046363号-1 
     


    收起
    展开
  • 四川郎酒股份有限公司获第十二届人民企业社会责任奖年度环保奖 2019-05-13
  • 银保监会新规剑指大企业多头融资和过度融资 2019-05-12
  • 韩国再提4国联合申办世界杯 中国网友无视:我们自己来 2019-05-11
  • 中国人为什么一定要买房? 2019-05-11
  • 十九大精神进校园:风正扬帆当有为 勇做时代弄潮儿 2019-05-10
  • 粽叶飘香幸福邻里——廊坊市举办“我们的节日·端午”主题活动 2019-05-09
  • 太原设禁鸣路段 设备在测试中 2019-05-09
  • 拜耳医药保健有限公司获第十二届人民企业社会责任奖年度企业奖 2019-05-08
  • “港独”没出路!“梁天琦们”该醒醒了 2019-05-07
  • 陈卫平:中国文化内涵包含三方面 文化复兴表现在其中 2019-05-06
  • 人民日报客户端辟谣:“合成军装照”产品请放心使用 2019-05-05
  • 【十九大·理论新视野】为什么要“建设现代化经济体系”?   2019-05-04
  • 聚焦2017年乌鲁木齐市老城区改造提升工程 2019-05-04
  • 【专家谈】上合组织——构建区域命运共同体的有力实践者 2019-05-03
  • 【华商侃车NO.192】 亲!楼市火爆,别忘了买车位啊! 2019-05-03
  • 北京pk10赛车怎样稳赢 龙虎和马 彩票可以买大小吗 秒速飞艇计划软件 大乐透杀号16法 最新后三不定位独胆稳赚 无错24码特围 如何买六肖稳赚 双色球走势图3000 pk10人工免费计划app 时时彩技巧100稳赚 功夫时时彩计划软件免费版 十一选五走势图分析技巧 一分快三在线稳赚计划 3d稳赚不赔绝招 大乐透质合走势图表图