CN113091600B - 一种利用时序InSAR技术监测软土地基形变的监测方法 - Google Patents
一种利用时序InSAR技术监测软土地基形变的监测方法 Download PDFInfo
- Publication number
- CN113091600B CN113091600B CN202110368285.4A CN202110368285A CN113091600B CN 113091600 B CN113091600 B CN 113091600B CN 202110368285 A CN202110368285 A CN 202110368285A CN 113091600 B CN113091600 B CN 113091600B
- Authority
- CN
- China
- Prior art keywords
- deformation
- soft soil
- insar
- phase
- time sequence
- 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.)
- Active
Links
- 239000002689 soil Substances 0.000 title claims abstract description 69
- 238000000034 method Methods 0.000 title claims abstract description 66
- 238000012544 monitoring process Methods 0.000 title claims abstract description 30
- 238000005516 engineering process Methods 0.000 title claims abstract description 23
- 230000002068 genetic effect Effects 0.000 claims abstract description 15
- 239000011159 matrix material Substances 0.000 claims description 24
- 238000001914 filtration Methods 0.000 claims description 8
- 230000035772 mutation Effects 0.000 claims description 8
- 238000012937 correction Methods 0.000 claims description 7
- 238000005457 optimization Methods 0.000 claims description 6
- 238000004062 sedimentation Methods 0.000 claims description 6
- 238000012545 processing Methods 0.000 claims description 4
- 210000000349 chromosome Anatomy 0.000 claims description 3
- 230000008929 regeneration Effects 0.000 claims description 3
- 238000011069 regeneration method Methods 0.000 claims description 3
- 238000012952 Resampling Methods 0.000 claims description 2
- 238000000605 extraction Methods 0.000 claims description 2
- 238000010587 phase diagram Methods 0.000 claims description 2
- 230000017105 transposition Effects 0.000 claims 1
- 230000007774 longterm Effects 0.000 abstract description 9
- 238000010276 construction Methods 0.000 abstract description 6
- 238000011161 development Methods 0.000 abstract description 4
- 238000007430 reference method Methods 0.000 abstract description 2
- 230000008569 process Effects 0.000 description 13
- 238000010586 diagram Methods 0.000 description 5
- 230000008859 change Effects 0.000 description 4
- 238000007596 consolidation process Methods 0.000 description 4
- 230000000694 effects Effects 0.000 description 4
- 238000004088 simulation Methods 0.000 description 4
- 230000007547 defect Effects 0.000 description 3
- 230000018109 developmental process Effects 0.000 description 3
- 230000002159 abnormal effect Effects 0.000 description 2
- 238000006243 chemical reaction Methods 0.000 description 2
- 201000010099 disease Diseases 0.000 description 2
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 2
- 239000006185 dispersion Substances 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- 238000005259 measurement Methods 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012216 screening Methods 0.000 description 2
- 230000006641 stabilisation Effects 0.000 description 2
- 238000011105 stabilization Methods 0.000 description 2
- 238000012360 testing method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000007812 deficiency Effects 0.000 description 1
- 238000013461 design Methods 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 238000013178 mathematical model Methods 0.000 description 1
- 238000000691 measurement method Methods 0.000 description 1
- 239000008239 natural water Substances 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 238000007781 pre-processing Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000001932 seasonal effect Effects 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
- 230000000007 visual effect Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01B—MEASURING LENGTH, THICKNESS OR SIMILAR LINEAR DIMENSIONS; MEASURING ANGLES; MEASURING AREAS; MEASURING IRREGULARITIES OF SURFACES OR CONTOURS
- G01B7/00—Measuring arrangements characterised by the use of electric or magnetic techniques
- G01B7/16—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge
- G01B7/24—Measuring arrangements characterised by the use of electric or magnetic techniques for measuring the deformation in a solid, e.g. by resistance strain gauge using change in magnetic properties
-
- 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/9004—SAR image acquisition techniques
-
- 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
- G01S19/00—Satellite radio beacon positioning systems; Determining position, velocity or attitude using signals transmitted by such systems
- G01S19/38—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system
- G01S19/39—Determining a navigation solution using signals transmitted by a satellite radio beacon positioning system the satellite radio beacon positioning system transmitting time-stamped messages, e.g. GPS [Global Positioning System], GLONASS [Global Orbiting Navigation Satellite System] or GALILEO
- G01S19/42—Determining position
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Radar, Positioning & Navigation (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Computer Networks & Wireless Communication (AREA)
- Electromagnetism (AREA)
- Geophysics And Detection Of Objects (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
本发明提供一种利用时序InSAR技术监测软土地基形变的监测方法，利用时序InSAR技术，结合高速公路软土地基沉降曲线发展规律，将公路沉降预测领域广泛应用的泊松曲线引入时序InSAR形变模型，建立了一种基于泊松曲线的时序InSAR形变模型，并提出一种将遗传算法与正则化牛顿迭代法串行的方法对构建的时序InSAR形变模型参数进行求解，从而估计软土覆盖区域的时序形变结果，可作为软土地基区域长期变形监测的借鉴方法，拓宽了InSAR技术的应用领域，也可为软土地基工后长期变形分析解译和形变预测提供一种参数估计的新方法。
Description
技术领域
本发明涉及软土地基地表形变监测及变形预计领域，特别地，涉及泊松曲线的利用时序InSAR技术监测软土地基形变的监测方法。
背景技术
由于软土具有天然含水量大、压缩性高、强度低、结构差等特性，软土地基更易发生沉降，导致变形失稳，其稳定性控制已成为公路路基工程中遇到的主要技术难题。开展软土地基工后长期沉降监测，对防治公路变形失稳，预防安全事故的发生，保证公路工程建设质量具有重要的工程意义。传统的地面沉降监测方法主要有水准测量、GPS测量等，尽管这些监测手段精度较高，但是其低空间分辨率、造价昂贵、且需要监测人员定期到监测现场的缺陷限制了其在公路形变监测中的应用。差分雷达干涉测量技术(DifferentialInterferometric Synthetic Aperture Radar，D-InSAR)虽可弥补传统测量手段的不足，但因其众所周知的时空失相关和大气延迟相位的影响却严重影响了其形变监测的精度。MT-InSAR技术是以永久散射体(Persistent Scatterer InSAR,PS-InSAR)、小基线集(Small Baseline Subsets InSAR,SBAS-InSAR)等技术为代表的高级InSAR技术，已被证实在测量人工线状地物、大型基础设施方面更具应用潜力，在国内外已取得了诸多先进的应用成果。
在InSAR数据处理过程中，形变模型的建立是至关重要的环节。目前，InSAR形变模型大多为单一数学经验模型(如线性模型、季节性模型、多项式模型等)或组合而成，并未考虑监测对象的具体变形机制，这就导致这些模型的拟合效果并不能反应软土地基公路沉降的真实时变规律，不但严重影响变形监测的精度，也不利于工后期的形变分析解译。软土沉降过程一般经过发生、发展、稳定后达一定极限，其表现形式呈现反“S”型曲线特征。泊松曲线为典型的“S”型增长曲线，具有单调递增性、不过原点性、有界性以及存在反弯点等特点，能较准确地反映软土路基固结沉降的变化规律，已被广泛应用于软土地基公路的工后沉降预测中。因此，将公路沉降预测领域广泛应用的泊松曲线引入时序InSAR形变模型，建立了一种基于泊松曲线的时序InSAR形变模型，可更合理的反应软土区形变随时间演变的规律，提高InSAR技术监测公路工后长期变形的精度，更合理地解译形变估计结果。
发明内容
本发明提供一种利用时序InSAR技术监测软土地基形变的监测方法，利用时序InSAR技术，结合高速公路软土地基沉降曲线发展规律，建立一种基于泊松曲线的时序InSAR形变模型。由于传统时序InSAR形变监测中采用的形变模型为纯经验数学模型，没有考虑软土地基沉降随时间发展演化的真实形变规律，因此无法实现软土形变的准确估计。软土沉降过程一般经过发生、发展、稳定后达一定极限，其表现形式呈现反“S”型曲线特征。泊松曲线为典型的“S”型增长曲线，具有单调递增性、不过原点性、有界性以及存在反弯点等特点，能较准确地反映软土路基固结沉降的变化规律，因此将泊松曲线引入时序InSAR，建立一种基于泊松曲线的时序InSAR形变模型，可改进传统时序InSAR形变模型的缺陷，可更合理的反应软土区形变随时间演变的规律，提高InSAR技术监测公路工后长期变形的精度，更合理地解译形变估计结果。
现有的普通牛顿迭代法的迭代公式为：
其中Gk为x(k)处的迭代矩阵。
普通牛顿迭代法要求矩阵Gk非奇异。在实际问题解算中，当Gk接近奇异或严重病态时，其求逆不稳定，此时式1)的迭代格式则不再适用，即普通牛顿迭代法无法进行求解。这种情况下，矩阵Gk的病态性可以通过正则化方法来改善，将迭代矩阵加上αkI(其中αk为正则因子)，即Gk+αkI，以替换表达式1)式中的迭代矩阵。通过适当选取αk，就可以改善原迭代矩阵Gk的病态性，有效保证Gk+αkI求逆运算的稳定性，从而促使迭代过程顺利进行。正则化的牛顿迭代公式为：
x(k+1)＝x(k)+dx(k)＝x(k)-(Gk+αkI)-1(g(k))T；
矩阵条件数是衡量矩阵病态的一个重要指标。上式中αkI的加入，使得Gk的条件数有效降低，改善了普通牛顿迭代法迭代过程中迭代矩阵的病态性导致的迭代不收敛的情况。由于GA搜索结果含有较大的随机误差，可能存在不成熟收敛，即没有收敛到全局最优解就结束迭代了；而单独采用RN算法进行参数的迭代解算，往往依赖于对解初值的设定，如果设置的初值与真实值偏差大，则迭代过程慢，且很可能不收敛。因此提出一种将遗传算法(Genetic Algorithm，GA)与正则化牛顿迭代法(Regularized Newton IterativeAlgorithm，RN)串行的方法对构建的时序InSAR形变模型参数进行求解，从而估计软土覆盖区域的时序形变结果，可作为软土地基区域长期变形监测的借鉴方法，拓宽了InSAR技术的应用领域，也可为软土地基工后长期变形分析解译和形变预测提供一种参数估计的新方法。
为实现上述目的，本发明提供一种利用时序InSAR技术监测软土地基形变的监测方法，包括以下步骤：
步骤一、基于高相干点的时序InSAR差分干涉图生成，具体是：利用时间序列InSAR技术获取软土地基解缠后的差分干涉相位图，并利用三重阈值法(平均相干系数、强度值和振幅离差指数)选取高相干点；
步骤二、构建InSAR干涉相位与泊松曲线的时间序列函数模型；
步骤三、软土InSAR形变模型参数估计，具体是：先利用遗传算法估计软土InSAR形变模型参数的初始值，再利用正则化牛顿迭代算法对遗传算法得到的软土InSAR形变模型参数进行优化；
步骤四、基于泊松曲线的软土InSAR形变模型参数的时序形变估计，具体是：估计InSAR影像干涉时段内的软土地基时序形变。
以上技术方案中优选的，步骤一包括如下步骤：
步骤1.1、对M+1幅SAR数据进行干涉组合、超级主影像选取、时空基线估计与阈值设定、影像配准及重采样，生成干涉图和相干图；M≥1；
步骤1.2、对步骤1.1中干涉图进行去轨道、平地相位和地形相位处理，生成差分干涉图，去除了轨道误差、平地和地形相位的差分干涉图生成；
步骤1.3、对步骤1.2中的差分干涉图进行相位解缠(部分干涉对进行3D解缠)，得到解缠后的差分干涉图；
步骤1.4、利用步骤1.1中生成的相干图和步骤1.3获得的解缠后的差分干涉图进行高相干点提取，生成基于高相干点目标的时序差分干涉相位矩阵。
以上技术方案中优选的，步骤二包括如下步骤：
泊松曲线为典型的“S”型增长曲线，能较准确地反映软土地基固结沉降的变化规律。其形式为表达式1)：
式中：W(t)为软土路基表面在t时刻的垂直向沉降值；W0为地表的最大沉降值；a和b为模型的形状参数，其中参数a越大，地表点的起始下沉时间越长，达到最大下沉值的时间也越久；参数b越大，地表点下沉速度越快，曲线形状之间差异也越大。在此，忽略水平移动的影响，将InSAR监测的LOS向形变转换为垂直向的时间序列沉降，即
(其中θ为雷达入射角)，则表达式1)写为表达式2)：
将传统线性模型的低通形变分量改写为线性形变分量和泊松曲线分量的组合，可得表达式3)：
其中：φ为时序InSAR得到的差分干涉相位；λ为信号的中心波长，tB为B时刻相对于起始时刻的时间，tA为A时刻相对于起始时刻的时间；v为线性速率，Δh为高程改正，B⊥为垂直基线长度；R和θ分别为传感器到目标的距离和雷达入射角；残余相位
由噪声相位、大气延迟相位和高频形变相位组成。
以上技术方案中优选的，步骤三包括以下步骤：
步骤3.1、利用遗传算法估计软土InSAR形变模型参数的初始值，包括以下步骤：
步骤3.1.1、根据残差最小原则，建立适应度函数
本发明中的残差可表示为表达式4)：
步骤3.1.2、产生初始种群，具体是：根据种群规模随机产生初始种群，每个个体表示染色体的基因型；
步骤3.1.3、计算适应度，具体是：计算每个个体的适应度函数，并判断是否满足优化准则，若满足则输出最佳个体及其代表的最优解，并结束算法；若不满足，则进入下一步；
步骤3.1.4、选择，具体是：依据适应度选择再生个体，适应度低的个体被选中的概率高，反之被选中的概率低，甚至被淘汰；
步骤3.1.5、交叉，具体是：根据一定的交叉概率和交叉方法生成子代个体；
步骤3.1.6、变异，具体是：根据一定的变异概率和变异方法，生成子代个体；
步骤3.1.7、循环计算适应度，具体是：由交叉和变异产生的新一代种群，返回步骤3.1.3，直到满足优化准则，输出最佳个体及其代表的最优解，作为高斯牛顿迭代法的初始值；
步骤3.2、利用正则化牛顿迭代算法对遗传算法得到的软土InSAR形变模型参数进行优化，具体是：将步骤3.1.7中遗传算法获取的估计值作为初值，利用正则化牛顿迭代法进行模型未知参数的求解，此处：
普通牛顿迭代法的迭代公式为表达式5)：
其中：Gk为x(k)处的迭代矩阵。普通牛顿迭代法要求Gk非奇异。在实际问题解算中，Gk接近奇异或严重病态时，其求逆不稳定，此时表达式5)的迭代格式则不再适用。矩阵的病态性可以通过正则化方法来改善，即单位矩阵I乘以正则因子αk后为αkI，迭代矩阵加上αkI，即Gk+αkI。适当选取αk，就可以保证Gk+αkI求逆运算的稳定性，从而促使迭代过程顺利进行。正则化的牛顿迭代公式为表达式6)：
x(k+1)＝x(k)+dx(k)＝x(k)-(Gk+αkI)-1(g(k))T 6)；
设x(0)为第一个迭代值，则x(1)＝x(0)+dx(0)；dx(k)为迭代步长，给定的允许误差率K(此处设为10-3)，当dx(k)＜K时输出结果，否则重复表达式6)，进行重复迭代，直到满足条件，结束算法，输出参数值。
以上技术方案中优选的，步骤四包括以下步骤：
步骤4.1、利用步骤3.2获取的未知参数W0、a、b、v及Δh带入表达式3)，得到高相干点低通形变分量部分；
根据大气延迟相位和噪声相位在时空域上的特性，对表达式3)中残差相位在时间维上进行高通滤波，空间维上进行低通滤波，从而提取出残余相位中的高通形变部分；
步骤4.2、将步骤4.1输出的低通形变分量与高通形变分量累加，获取高相干点视线向总形变；
步骤4.3、将步骤4.2生成的视线向形变进行地理编码，生成软土地基垂直向时序形变场。
应用本发明的技术方案，具有以下有益效果：
本发明利用时间序列InSAR技术，将软土沉降预测领域广泛应用的泊松曲线引入时序InSAR形变模型，构建InSAR干涉相位与软土形变参数的函数关系模型，并求解了模型中的未知参数，进而估计软土地基区域的时间序列形变，可应用于软土地基区域长期的地表变形监测中，克服了单一数学经验模型未考虑监测对象的具体变形机制所导致模型的拟合效果并不能反应软土地基公路沉降的真实时变规律的缺陷，可更合理的反应软土区形变随时间演变的规律，提高InSAR技术监测公路工后长期变形的精度，更合理地解译形变估计结果。
除了上面所描述的目的、特征和优点之外，本发明还有其它的目的、特征和优点。下面将参照图，对本发明作进一步详细的说明。
附图说明
构成本申请的一部分的附图用来提供对本发明的进一步理解，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
图1为实施例中模型构建及参数估计解算流程图；
图2为模拟出的地表最大沉降值W0的真实值；
图3为模拟出的形状参数a的真实值；
图4为模拟出的形状参数b的真实值；
图5为基于泊松曲线模拟出的软土地基时序形变结果图；
图6为高相干点上参数真实值与预测值对比图，其中：(a)为最大沉降值W0的计算值和真值比对图；(b)为形状参数a的计算值和真值比对图；(c)为形状参数b的计算值和真值比对图；(d)为线性速率v的计算值和真值比对图；(e)为高程改正Δh的计算值和真值比对图；
图7为高相干点上时序形变量真实值与预测值对比图；
图8为高相干点上时序形变量真实值与预测值残差的均方根误差图。
具体实施方式
以下结合附图对本发明的实施例进行详细说明，但是本发明可以根据权利要求限定和覆盖的多种不同方式实施。
实施例：
一种利用时序InSAR技术监测软土地基形变的监测方法，即软土地基InSAR形变监测方法，具体实施步骤如图1，详情如下：
第一步、基于高相干点的时序InSAR差分干涉图生成，即利用时序InSAR技术获取待监测矿区解缠后的差分干涉图并选取高相干点，详情如下：
步骤1.1、前期数据处理。将收集的多幅(可根据实际需求选择)SAR卫星数据批量进行格式转换，目的是生成可由处理软件读取的单视复数(SLC)格式，并进行初始影像的预滤波，减弱影像的噪声相位。格式转换中，每景影像都会生成配套的参数文件。处理软件会自动选取一景影像作为所有影像的超级主影像，设置时空阈值，将生成所有符合设定的干涉对。在此基础上，所有的影像会配准并重采样到超级主影像，配准过程中采用多项式拟合。利用之前生成的可用干涉组合依次生成干涉图，在此过程中会附带生成相干系数图和强度图。然后借助精密轨道数据去除平地相位、轨道误差；利用外部30m分辨率的DEM数据去除地形相位；再进行残余相位滤波，去除影像中的噪声影响。最后对经过前期处理完的干涉相位进行相位解缠(解缠方法的原则在于：满足3D解缠条件的干涉图进行3D解缠，其余则采用最小费用流法)，进而生成重去平后的时间序列差分干涉相位，本实施例中模拟生成的干涉图为重去平的干涉图。
步骤1.2、提取测试区域的高相干点并组建所有优点的时间序列差分相位矩阵。SBAS技术是在严格控制时空基线阈值，挑选出所有可用的干涉对，并通过筛选高质量点进行建模分析，从而提取出形变分量，反演整个测区时序形变结果的过程。其中，采用三重阈值法(强度、振幅离差指数和相干系数)反复筛选高相干点，并进行人工筛选，去除一些不符合常理的高相干点。本次实例中高相干点的坐标位置利用模拟实验中的随机模拟器生成，在此基础上，提取出对应在步骤1.1中重去平后的差分干涉相位，组建成时序高相干点差分干涉矩阵。
第二步、构建InSAR干涉相位与泊松曲线的时间序列函数模型，详情如下：
泊松曲线为典型的“S”型增长曲线，能较准确地反映软土地基固结沉降的变化规律。其形式为表达式1)：
式中，W(t)为软土路基表面在t时刻的垂直向沉降值；W0为地表的最大沉降值；a和b为模型的形状参数，其中参数a越大，地表点的起始下沉时间越长，达到最大下沉值的时间也越久；参数b越大，地表点下沉速度越快，曲线形状之间差异也越大。
在此，忽略水平移动的影响，将InSAR监测的LOS向形变转换为垂直向的时间序列沉降，即
(其中θ为雷达入射角)，则表达式1)可写为表达式2)：
将传统线性模型的低通形变分量改写为线性形变分量和泊松曲线分量的组合，可得表达式3)：
其中：φ为时序InSAR得到的差分干涉相位；λ为信号的中心波长，tB为B时刻相对于起始时刻的时间，tA为A时刻相对于起始时刻的时间；v为线性速率，Δh为高程改正，B⊥为垂直基线长度；R和θ分别为传感器到目标的距离和雷达入射角；残余相位
由噪声相位、大气延迟相位和高频形变相位组成。
对上式，其中未知参数W0，a，b，线性速率v，高程改正Δh。当方程个数大于或者等于5时，即可解算出影像上所有高相干点的未知参数。
第三步、遗传算法结合正则化牛顿迭代法进行软土InSAR形变模型参数估计，具体是：先利用遗传算法估计软土InSAR形变模型参数的初始值，再利用正则化牛顿迭代算法对遗传算法得到的软土InSAR形变模型参数进行优化，详情如下：
步骤3.1、利用遗传算法估计软土InSAR形变模型参数的初始值，包括以下步骤：
步骤3.1.1、根据残差最小原则，建立适应度函数
本发明中的残差可表示为表达式4)：
步骤3.1.2、(产生初始种群)根据种群规模。随机产生初始种群，每个个体表示染色体的基因型。
步骤3.1.3、(计算适应度)计算每个个体的适应度函数，并判断是否满足优化准则，若满足则输出最佳个体及其代表的最优解，并结束算法；若不满足，则转入下一步。
步骤3.1.4、(选择)依据适应度选择再生个体，适应度低的个体被选中的概率高，反之被选中的概率低，甚至被淘汰。
步骤3.1.5、(交叉)根据一定的交叉概率和交叉方法生成子代个体；
步骤3.1.6、(变异)根据一定的变异概率和变异方法，生成子代个体；
步骤3.1.7、(循环计算适应度)由交叉和变异产生的新一代种群，返回步骤3.1.3，直到满足优化准则，输出最佳个体及其代表的最优解，作为高斯牛顿迭代法的初始值；
步骤3.2、(正则化牛顿迭代法)将步骤3.1.7中遗传算法获取的估计值作为初值，利用正则化牛顿迭代法进行模型未知参数的求解。普通牛顿迭代法的迭代公式为表达式5)：
其中：Gk为x(k)处的迭代矩阵。普通牛顿迭代法要求Gk非奇异。在实际问题解算中，Gk接近奇异或严重病态时，其求逆不稳定，此时表达式5)的迭代格式则不再适用。由相关文献可知，矩阵的病态性可以通过正则化方法来改善。即单位矩阵I乘以正则因子αk后为αkI。迭代矩阵加上αkI，即Gk+αkI。适当选取αk，就可以保证Gk+αkI求逆运算的稳定性，从而促使迭代过程顺利进行。正则化的牛顿迭代公式为表达式6)：
x(k+1)＝x(k)+dx(k)＝x(k)-(Gk+αkI)-1(g(k))T 6)；
设x(0)为第一个迭代值，则x(1)＝x(0)+dx(0)；dx(k)为迭代步长，给定的允许误差率K(本专利设为10-3)，当dx(k)＜K时输出结果，否则重复表达式6)，进行重复迭代，直到满足条件，结束算法，输出参数值。
第四步、基于泊松曲线的软土InSAR形变模型参数的时序形变估计，具体包括如下步骤：
步骤4.1、利用步骤3.2获取的未知参数W0，a，b，v，Δh带入表达式3)，得到高相干点低通形变分量部分；
步骤4.2、根据大气延迟相位和噪声相位在时空域上的特性，对表达式3)中残差相位在时间维上进行高通滤波，空间维上进行低通滤波，从而提取出残余相位中的高通形变部分；
步骤4.3、将步骤4.1输出的低通形变分量与步骤4.2输出的高通形变分量累加，获取高相干点视线向总形变；
步骤4.4、将步骤4.3生成的视线向形变进行地理编码，生成软土地基垂直向时序形变场，详见图5。
应用本实施例的方案，详情如下：
通过对测试区域的设计资料及结构形态调查，在本实施例中，设置最大沉降量W0取值区间为[-0.1，0]m，形状参数a的取值区间为[0，100]，形状参数b的取值区间为[2，8]，线性速率v取值位于区间[-0.2，0.1]m/yr。四个参数的真实参数场均利用二维高斯函数模型进行模拟，高程改正值Δh则利用二维高斯随机函数模拟器模拟，取值范围为[-15，15]m。在模拟实验中，采用的卫星参数均依据真实测区(TerraSAR-X波段，Stripmap降轨模式)的影像参数设定。根据现有真实实验结果，选取其中12个干涉效果最佳的干涉对参数作为本实施例中的时空基线参数。在预先模拟的真实参数场中随机提取出200个高相干点的参数值作为真值，用于对前述模型及GARN算法估计出的模型参数结果进行验证。图2为模拟出的地表最大沉降值W0的真实值；图3为模拟出的形状参数a的真实值；图4为模拟出的形状参数b的真实值；图5为基于泊松曲线模拟出的软土地基时序形变结果图；(参考时间为2015年1月1日)；图6表示实验获取的200个高相干点参数真值与预测值的比较结果(噪声水平为0.5rad)，其中(a)示意了最大沉降值W0；(b)示意了形状参数a；(c)示意了形状参数b；(d)示意了线性速率v；(e)示意了高程改正Δh，由图2-图6可知，在一定噪声的干扰下，本发明估计出的泊松曲线参数与真值吻合情况仍然良好。图7为高相干点上时序形变量真实值与预测值对比图，图8为高相干点上时序形变量真实值与预测值残差的均方根误差图，由图7和图8可知，在一定噪声的干扰下，本发明估计的形变量与真实形变量基本保持一致。
以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (4)
1.一种利用时序InSAR技术监测软土地基形变的监测方法，其特征在于，包括以下步骤：
步骤一、基于高相干点的时序InSAR差分干涉图生成，具体是：利用时间序列InSAR技术获取软土地基解缠后的差分干涉相位图，并利用三重阈值法选取高相干点；
步骤二、构建InSAR干涉相位与泊松曲线的时间序列函数模型；
步骤三、软土InSAR形变模型参数估计；
步骤四、基于泊松曲线的软土InSAR形变模型参数的时序形变估计，具体是：估计InSAR影像干涉时段内的软土地基时序形变；
其中，步骤三包括以下步骤：
步骤3.1、利用遗传算法估计软土InSAR形变模型参数的初始值，包括以下步骤：
步骤3.1.1、根据残差最小原则，建立适应度函数
残差
表示为表达式4)：
步骤3.1.2、产生初始种群，具体是：根据种群规模随机产生初始种群，每个个体表示染色体的基因型；
步骤3.1.3、计算适应度，具体是：计算每个个体的适应度函数，并判断是否满足优化准则，若满足则输出最佳个体及其代表的最优解，并结束算法；若不满足，则进入下一步；
步骤3.1.4、选择，具体是：依据适应度选择再生个体，适应度低的个体被选中的概率高，反之被选中的概率低，甚至被淘汰；
步骤3.1.5、交叉，具体是：根据一定的交叉概率和交叉方法生成子代个体；
步骤3.1.6、变异，具体是：根据一定的变异概率和变异方法，生成子代个体；
步骤3.1.7、循环计算适应度，具体是：由交叉和变异产生的新一代种群，返回步骤3.1.3，直到满足优化准则，输出最佳个体及其代表的最优解，作为高斯牛顿迭代法的初始值；
步骤3.2、利用正则化牛顿迭代算法对遗传算法得到的软土InSAR形变模型参数进行优化，具体是：将步骤3.1.7中遗传算法获取的估计值作为初值，利用正则化牛顿迭代法表达式6)进行模型未知参数的求解：
x(k+1)＝x(k)+dx(k)＝x(k)-(Gk+αkI)-1(g(k))T 6)；
其中，x(k+1)为第k+1次迭代值，x(k)为第k次迭代值，dx(k)为迭代步长；g(k)为雅可比矩阵，即模型对未知参数的一阶偏导；Gk为海森矩阵，即模型对未知参数的二阶偏导；αk为正则因子，I为单位矩阵，T表示为将g(k)转置；设x(0)为第一个迭代值，则x(1)＝x(0)+dx(0)，给定的允许误差率K，当dx(k)<K时输出结果，否则重复表达式5)进行重复迭代，直到满足条件，结束算法，输出每一个最优参数值。
2.根据权利要求1所述的监测方法，其特征在于，步骤一包括如下步骤：
步骤1.1、对M+1幅SAR数据进行干涉组合、超级主影像选取、时空基线估计与阈值设定、影像配准及重采样，生成干涉图和相干图；M≥1；
步骤1.2、对步骤1.1中干涉图进行去轨道、平地相位和地形相位处理，生成差分干涉图，去除了轨道误差、平地和地形相位的差分干涉图生成；
步骤1.3、对步骤1.2中的差分干涉图进行相位解缠，得到解缠后的差分干涉图；
步骤1.4、利用步骤1.1中生成的相干图和步骤1.3获得的解缠后的差分干涉图进行高相干点提取，生成基于高相干点目标的时序差分干涉相位矩阵。
3.根据权利要求1所述的监测方法，其特征在于，步骤二具体是：
步骤2.1、构建InSAR干涉相位与软土形变参数的函数关系模型：
泊松曲线的形式为表达式1)：
式中：W(t)为软土路基表面在t时刻的垂直向沉降值；W0为地表的最大沉降值；a和b为模型的形状参数，其中参数a越大，地表点的起始下沉时间越长，达到最大下沉值的时间也越久；参数b越大，地表点下沉速度越快，曲线形状之间差异也越大；
忽略水平移动的影响，将InSAR监测的LOS向形变转换为垂直向的时间序列沉降，即
则表达式1)写为表达式2)：
其中：θ为雷达入射角；dLOS(t)为t时刻LOS向形变值；
将传统线性模型的低通形变分量改写为线性形变分量和泊松曲线分量的组合，得到表达式3)：
其中：φ为时序InSAR得到的差分干涉相位；λ为信号的中心波长，tB为B时刻相对于起始时刻的时间，tA为A时刻相对于起始时刻的时间；v为线性速率，Δh为高程改正，B⊥为垂直基线长度；R和θ分别为传感器到目标的距离和雷达入射角；残余相位
由噪声相位、大气延迟相位和高频形变相位组成。
4.根据权利要求1所述的监测方法，其特征在于，步骤四包括以下步骤：
步骤4.1、利用步骤3.2中获取的最优参数估值W0、a、b、v及Δh带入表达式3)，得到高相干点低通形变分量部分；
根据大气延迟相位和噪声相位在时空域上的特性，对表达式3)中残差相位在时间维上进行高通滤波，空间维上进行低通滤波，从而提取出残余相位中的高通形变部分；
步骤4.2、将步骤4.1输出的低通形变分量与高通形变分量累加，获取高相干点视线向总形变；
步骤4.3、将步骤4.2生成的视线向形变进行地理编码，生成软土地基垂直向时序形变场。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110368285.4A CN113091600B (zh) | 2021-04-06 | 2021-04-06 | 一种利用时序InSAR技术监测软土地基形变的监测方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202110368285.4A CN113091600B (zh) | 2021-04-06 | 2021-04-06 | 一种利用时序InSAR技术监测软土地基形变的监测方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN113091600A CN113091600A (zh) | 2021-07-09 |
| CN113091600B true CN113091600B (zh) | 2022-12-16 |
Family
ID=76674598
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202110368285.4A Active CN113091600B (zh) | 2021-04-06 | 2021-04-06 | 一种利用时序InSAR技术监测软土地基形变的监测方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN113091600B (zh) |
Families Citing this family (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN114325693A (zh) * | 2021-12-19 | 2022-04-12 | 南京市测绘勘察研究院股份有限公司 | 一种基于InSAR时序形变结果的采空区中心形变预测方法 |
| CN115201822B (zh) * | 2022-07-07 | 2023-03-14 | 长沙理工大学 | 一种钻井水溶岩盐矿区采卤量估计方法 |
Family Cites Families (5)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| IT1393118B1 (it) * | 2009-03-04 | 2012-04-11 | Eni Spa | Apparato e metodo di misura di spostamenti spaziali di strutture impiantistiche |
| CN101770027B (zh) * | 2010-02-05 | 2012-05-16 | 河海大学 | 基于InSAR与GPS数据融合的地表三维形变监测方法 |
| CN111323776B (zh) * | 2020-02-27 | 2021-04-13 | 长沙理工大学 | 一种矿区形变的监测方法 |
| CN111622274A (zh) * | 2020-05-08 | 2020-09-04 | 中交路桥建设有限公司 | 山区巨粒土高填方地基基础沉降预测方法及*** |
| CN111856465B (zh) * | 2020-07-30 | 2021-12-28 | 电子科技大学 | 一种基于稀疏约束的前视海面目标角超分辨方法 |
-
2021
- 2021-04-06 CN CN202110368285.4A patent/CN113091600B/zh active Active
Also Published As
| Publication number | Publication date |
|---|---|
| CN113091600A (zh) | 2021-07-09 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN104123464B (zh) | 一种高分辨率InSAR时序分析反演地物高程与地面沉降量的方法 | |
| CN113091600B (zh) | 一种利用时序InSAR技术监测软土地基形变的监测方法 | |
| CN108802727B (zh) | 一种顾及流变参数的时序InSAR公路形变监测模型及解算方法 | |
| CN111323776B (zh) | 一种矿区形变的监测方法 | |
| CN114966685B (zh) | 基于InSAR和深度学习的大坝形变监测及预测方法 | |
| CN108645994B (zh) | 一种基于多点地质统计学的地质随机反演方法及装置 | |
| CN104597490A (zh) | 基于精确Zoeppritz方程的多波AVO储层弹性参数反演方法 | |
| CN112113545B (zh) | 一种基于多维海面信息的内波振幅反演方法 | |
| CN111398959A (zh) | 基于地表应力应变模型的InSAR时序地表形变监测方法 | |
| CN108919266A (zh) | 一种基于PSInSAR技术的桥梁安全预警方法 | |
| CN113189559B (zh) | 一种星载成像高度计遥感数据海底地形反演方法 | |
| CN113281749A (zh) | 一种顾及同质性的时序InSAR高相干点选取方法 | |
| CN112051572A (zh) | 一种融合多源sar数据三维地表形变监测方法 | |
| Behzadi et al. | Introducing a novel digital elevation model using artificial neural network algorithm | |
| CN116068511B (zh) | 一种基于深度学习的InSAR大尺度***误差改正方法 | |
| CN116224327B (zh) | 一种基于学习网络的矿区大梯度形变区相位解缠方法 | |
| WO2022153984A1 (ja) | 学習データ生成方法、モデル生成方法および学習データ生成装置 | |
| CN115795952A (zh) | 一种基于有限元和贝叶斯克里金的边坡变形场重构方法 | |
| Abdalla et al. | Augmentation of vertical accuracy of digital elevation models using Gaussian linear convolution filter | |
| CN110794469B (zh) | 基于最小地质特征单元约束的重力反演方法 | |
| CN112268517B (zh) | 一种利用PSInSAR监测输电铁塔设备变形的方法 | |
| CN110857999B (zh) | 一种基于全波形反演的高精度波阻抗反演方法及*** | |
| CN115201822B (zh) | 一种钻井水溶岩盐矿区采卤量估计方法 | |
| CN116363057B (zh) | 一种融合PCA和时序InSAR的地表形变提取方法 | |
| CN117111057A (zh) | 一种煤矿采空区形变敏感性评价方法 |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| CB03 | Change of inventor or designer information |
Inventor after: Xing Xuemin Inventor after: Zhu Lingjie Inventor after: Peng Wei Inventor after: Zhu Jun Inventor after: Liu Bin Inventor after: Zhang Tengfei Inventor after: Huang Li Inventor after: Zhang Jihang Inventor before: Xing Xuemin Inventor before: Zhu Lingjie Inventor before: Zhang Tengfei Inventor before: Huang Li Inventor before: Zhang Jihang |
|
| CB03 | Change of inventor or designer information | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant | ||
| TR01 | Transfer of patent right |
Effective date of registration: 20231019 Address after: Room 1001, Comprehensive Building, Futian Xingye Building, No. 377 Jinma Road, Furong North Road Street, Kaifu District, Changsha City, Hunan Province, 410000 Patentee after: Hunan Zhongkan Beidou Research Institute Co.,Ltd. Address before: Changsha University of science and technology Patentee before: CHANGSHA University OF SCIENCE AND TECHNOLOGY |
|
| TR01 | Transfer of patent right |