基于方差分量估计与应力应变模型的InSAR三维地表形变监测方法技术

本发明专利技术公开了一种基于方差分量估计与应力应变模型的InSAR三维地表形变监测方法，首先，利用待监测地区的LOS向和方位向形变值，基于应力应变模型建立地表待监测点的三维形变与其所在窗口范围内的多个点的InSAR观测量之间的观测方程；接着，根据上一步骤确定的观测方程按照方差分量估计进行迭代定权；最后，将所得权重代入法方程即可求得该像素点精确的三维地表形变，对所有像素进行上述步骤，可得到该区域精确的三维地表形变场。该方法在少量InSAR观测量的情况下实现了三维地表形变场的精确估计。

【技术实现步骤摘要】
基于方差分量估计与应力应变模型的InSAR三维地表形变监测方法
本专利技术涉及一种基于方差分量估计与应力应变模型的InSAR三维地表形变监测方法。
技术介绍
合成孔径雷达干涉测量(InterferometricSyntheticApertureRadar,InSAR)在监测地表形变方面具有全天候，全天时，大范围，高精度与高空间分辨率等特点，然而根据目前的技术仅能获取雷达视线方向(LightofSight，LOS)与方位向(即卫星飞行方向)投影形变，要想获得三维地表形变则至少需要将两个卫星轨道得到的LOS向形变与方位向形变进行融合处理。利用差分InSAR(DInSAR)技术获得的LOS向形变投影是通过两幅全孔径的SAR单视复数影像处理计算得到，其精度可达到毫米级或厘米级；利用多孔径InSAR(MAI)或者偏移量跟踪法(OFT)获得的方位向投影形变，精度仅在厘米级或者分米级。因此，不同类型的InSAR形变观测量在精度方面相差几倍甚至十几倍。此外，不同卫星平台、不同波段数据得到的InSAR形变观测量，其观测误差也不尽相同。因此，在融合不同类型、不同平台的InSAR数据时，精确确定它们之间的权重关系对于获取高精度的三维地表形变结果十分重要。目前，不同类型数据定权的主要途径是利用一个合适的模板窗口(如5×5)在不同类型数据矩阵中遍历，将各类数据中模板窗口内数值的方差大小作为该类数据的先验方差，再进行定权。但是此类定权方法需要假设窗口内各点的形变具有各态历经性，这个条件一般难以满足，因此其定权的精度有限。定权的另外一个途径是根据方差分量估计来进行验后定权，此类方法理论上可得到各类观测量的精确权重，但是利用方差分量进行验后定权时需要较多的多余观测，单凭传统的解算方法则需要大量的时序数据来获取多余观测，对于瞬时形变(如火山、地震等)而言并不适用。
技术实现思路
本专利技术提出一种基于方差分量估计与应力应变模型的InSAR三维地表形变监测方法，该方法在解算三维地表形变时考虑地表应力应变关系，更加准确地反映三维地表形变之间的关系，同时，应用方差分量估计使得在融合多平台、多源InSAR数据时能够根据不同数据的方差水平来决定其在解算三维地表形变中的权重问题，得到了更为精确的三维地表形变场。一种基于方差分量估计与应力应变模型的InSAR三维地表形变监测方法，包括以下步骤：步骤1：利用待监测地区的LOS向和方位向形变值，基于应力应变模型建立地表待监测点的三维形变与其所在窗口范围内的多个点的InSAR观测量之间的观测方程；地表待监测点所在窗口大小为5×5-25×25像素；应力应变模型是指将地表一定深度的陆地部分看作刚体，当发生地表地震、火山等自然灾害时可以看做有外力作用在该部分刚体上，该外力近似包含两部分作用，一是使刚体发生拉压形变，二是使刚体发生扭转形变。步骤2：基于步骤1构建的观测方程，构建升降轨LOS向与方位向观测值关于待监测点三维形变之间的误差方程组，并结合方差分量估计求解4类观测值单位权中误差的初步解；所述4类观测值依次包括升轨LOS向形变、升轨方位向形变、降轨LOS向形变与降轨方位向形变；步骤3：利用方差分量估计算法对各类观测值权阵进行更新，将更新后的各类观测值权阵重新代入步骤2中求解各类观测值单位权中误差，直到每类观测值单位权中误差之间的差值不超过1mm2，停止迭代，求得待监测点的三维地表形变的最优解；在迭代过程中目的是为了得到更为精确的单位权中误差，但是迭代的同时也会计算未知数的值，那么当迭代结束的时候，所对应的未知数的值即为最优方程解；步骤4：对地表其余像素点重复上述步骤1-3，完成待监测区域的三维形变场的计算。进一步地，所述步骤1中的观测方程如下：Li＝Bi·t其中，Li表示地表待监测点所在窗口范围内的第i个像素点对应的InSAR观测值，及分别代表地表待监测点所在窗口范围内的第i个像素点的升轨LOS向、升轨方位向、降轨LOS向及降轨方位向的观测值；Bgeo为是卫星成像的几何投影系数矩阵，Bgeo_1、Bgeo_2、Bgeo_3及Bgeo_4依次为4类观测值对应的卫星成像的几何投影系数，为地表待监测点与其所在窗口范围内的第i个像素点之间的应力应变模型系数矩阵；t表示未知数向量，及分别表示地表监测点的三维形变量，ξpq代表应力应变模型中的拉压参数，1≤p≤3，p≤q≤3，ωr表示旋转参数，r＝1,2,3。进一步地，所述4类观测值单位权中误差的初步解的求解过程如下：步骤A：首先求解未知数向量t的初值：t＝N-1U；其中，j的取值范围为1-4，i的取值范围为1-n，n表示从地表待监测点所在窗口范围内选取的像素点个数；Wj表示第j类观测值权阵，初始值为n×n的单位矩阵I；步骤B：利用公式vj＝Bj·t-Lj得到各类观测值对应的改正数v1,v2,v3,v4，并获取基于各类观测值改正数的二次型向量δ，步骤C：计算各类观测值单位权中误差的初步解：其中，Γ为系数对称矩阵。系数矩阵是通过观测值改正数的统计信息从概率论的角度推导得出，若有m类观测值，则系数矩阵即为一个m×m的对称矩阵，对于主对角线元素，位于(i，i)位置处，元素取值为n-2tr(N-1Ni)+tr(N-1Ni)2，非主对角元素，上半部分中位于(i，j)处，元素取值为tr(N-1NiN-1Nj)，非主对角元素下半部分与上半部分对称相等。进一步地，所述利用方差分量估计算法对各类观测值权阵进行更新是按照以下公式进行更新计算：其中，表示更新后的第j类观测值权阵。有益效果本专利技术提供了一种基于方差分量估计与应力应变模型的InSAR三维地表形变监测方法，该方法在解算三维地表形变时考虑了地表应力应变关系，更加准确地反映三维地表形变之间的关系，与此同时，建立了关于中心像素点更多的观测方程(即多余观测增多)，为方差分量估计的应用提供了契机。相比于传统解法，本专利技术考虑了应力应变模型，可更加准确地反映地表三维形变场，与此同时应用方差分量估计使得在融合多平台、多源InSAR数据时能够根据不同数据的方差水平来决定其在解算三维地表形变中的权重问题，实现了不同类型观测量权重的精确确定，大大提高了InSAR三维地表形变监测的精度。此外，本专利技术在求解某一像素点的三维形变时，利用到了周围像素点的InSAR观测值，因此即使在中心像素点数据缺失的情况下依旧可以解算求得该像素点的三维形变。附图说明图1是本专利技术所述方法的流程图；图2是升降轨LOS向与方位向的模拟形变数据示意图，其中，(a)为升轨LOS向，(b)为升轨AZI向，(c)为降轨LOS向，(d)为降轨AZI向；图3是本专利技术方法、传统算法解算得到的三维形变场与原始模拟数据的对比图，其中，(a)为东西向原始模拟数据，(b)为南北向原始模拟数据，(c)为垂直向原始模拟数据，(d)为应用本专利技术得到的东西向三维形变场数据，(e)为应用本专利技术得到的南北向三维形变场数据，(f)为应用本专利技术方法得到的垂直向三维形变场数据，(g)为应用传统算法得到的东西向三维形变场数据，(h)为应用传统算法得到的南北向三维形变场数据，(i)为应用传统算法得到的垂直向三维形变场数据。具体实施方式下面将结合附图和实施例对本专利技术做进一步的说明。如图1所示，一种基于方差分量估计与应力应变模型本文档来自技高网...

【技术保护点】
一种基于方差分量估计与应力应变模型的InSAR三维地表形变监测方法，其特征在于，包括以下步骤：步骤1：利用待监测地区的LOS向和方位向形变值，基于应力应变模型建立地表待监测点的三维形变与其所在窗口范围内的多个点的InSAR观测量之间的观测方程；地表待监测点所在窗口大小为5×5‑25×25像素；步骤2：基于步骤1构建的观测方程，构建升降轨LOS向与方位向观测值关于待监测点三维形变之间的误差方程组，并结合方差分量估计求解4类观测值单位权中误差的初步解；所述4类观测值依次包括升轨LOS向形变、升轨方位向形变、降轨LOS向形变与降轨方位向形变；步骤3：利用方差分量估计算法对各类观测值权阵进行更新，将更新后的各类观测值权阵重新代入步骤2中求解各类观测值单位权中误差，直到每类观测值单位权中误差之间的差值不超过1mm

【技术特征摘要】
1.一种基于方差分量估计与应力应变模型的InSAR三维地表形变监测方法，其特征在于，包括以下步骤：步骤1：利用待监测地区的LOS向和方位向形变值，基于应力应变模型建立地表待监测点的三维形变与其所在窗口范围内的多个点的InSAR观测量之间的观测方程；地表待监测点所在窗口大小为5×5-25×25像素；步骤2：基于步骤1构建的观测方程，构建升降轨LOS向与方位向观测值关于待监测点三维形变之间的误差方程组，并结合方差分量估计求解4类观测值单位权中误差的初步解；所述4类观测值依次包括升轨LOS向形变、升轨方位向形变、降轨LOS向形变与降轨方位向形变；步骤3：利用方差分量估计算法对各类观测值权阵进行更新，将更新后的各类观测值权阵重新代入步骤2中求解各类观测值单位权中误差，直到每类观测值单位权中误差之间的差值不超过1mm2，停止迭代，求得待监测点的三维地表形变的最优解；步骤4：对地表其余像素点重复上述步骤1-3，完成待监测区域的三维形变场的计算。2.根据权利要求1所述的方法，其特征在于，所述步骤1中的观测方程如下：Li＝Bi·t其中，Li表示地表待监测点所在窗口范围内的第i个像素点对应的InSAR观测值，及分别代表地表待监测点所在窗口范围内的第i个像素点的升轨LOS向、升轨方位向、降轨LOS向及降轨方位向的观测值；Bgeo为是卫星成像的几何投影系数矩阵，Bgeo...

【专利技术属性】
技术研发人员：胡俊，刘计洪，李志伟，朱建军，
申请(专利权)人：中南大学，
类型：发明
国别省市：湖南,43