基于方差分量估计与应力应变模型的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中求解各类观测值单位权中误差，直到每类观测值单位权中误差之间的差值不超过1mm²，停止迭代，求得待监测点的三维地表形变的最优解；

在迭代过程中目的是为了得到更为精确的单位权中误差，但是迭代的同时也会计算未知数的值，那么当迭代结束的时候，所对应的未知数的值即为最优方程解；

步骤4：对地表其余像素点重复上述步骤1-3，完成待监测区域的三维形变场的计算。

进一步地，所述步骤1中的观测方程如下：

lⁱ＝bⁱ·t

其中，lⁱ表示地表待监测点所在窗口范围内的第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)²，非主对角元素，上半部分中位于(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观测量之间的观测方程；

选取同一地区的升轨和降轨sar影像对进行预处理、配准、干涉、差分、滤波、解缠及地理编码等过程；利用dinsar和mai/oft方法得到该地区的一对los向和一对方位向形变值。

地表待监测点所在窗口大小为5×5-25×25像素；

应力应变模型是指将地表一定深度的陆地部分看作刚体，当发生地表地震、火山等自然灾害时可以看做有外力作用在该部分刚体上，该外力近似包含两部分作用，一是使刚体发生拉压形变，二是使刚体发生扭转形变。本实例中仅在小范围(约150m×150m，对于30米分辨率的数据即为5×5像素)内认为地球表面陆地为理想刚体。

步骤2：基于步骤1构建的观测方程，构建升降轨los向与方位向观测值关于待监测点三维形变之间的误差方程组，并结合方差分量估计求解4类观测值单位权中误差的初步解；

所述4类观测值依次包括升轨los向形变、升轨方位向形变、降轨los向形变与降轨方位向形变；

假设地表有邻近两点pⁱ[xⁱyⁱzⁱ]^t、p⁰[x⁰y⁰z⁰]^t,其形变量分别为令dxⁱ＝xⁱ-x⁰,dyⁱ＝yⁱ-y⁰,dzⁱ＝zⁱ-z⁰，那么根据弹性理论有下式：

dⁱ＝h·δ+d⁰(1)

其中δ＝[dxⁱdyⁱdzⁱ]^t，对于后续公式及变量，若无特殊说明上标数字或字母均代表不同的像素点；h代表应力应变模型参数矩阵，可表示为

e为地表拉压应力产生的拉压参数矩阵，r为地表剪切应力产生的扭转参数矩阵。

将上述式(1)可写成：

其中，

代表应力应变模型(deformationgradienttensor)参数矩阵，

代表未知参数向量，ξpq代表应力应变模型中的拉压参数(1≤p≤3,p≤q≤3)，ωr代表旋转参数(r＝1,2,3)。

再假设在像素i上对应的insar观测值为(下标1,2,3,4分别代表升轨los向、升轨方位向、降轨los向及降轨方位向四类不同观测值类型)，则根据sar卫星的成像几何可得：

lⁱ＝bgeo·dⁱ(3)

其中令

则

式中α^as,θ^as,αdes,θdes分别为升轨卫星方位角与入射角和降轨卫星方位角与入射角,bgeo为卫星成像几何(geometry)参数系数矩阵。

综合上述函数关系式(2)(3)可得：lⁱ＝bⁱ·t(4)

其中，

这样即建立了某像素点周围一点insar观测值与该像素点三维地表形变之间模型关系。

假设该像素点周围(以该像素为中心的5×5窗口内)有n个像素点具有insar观测值(包含该像素点本身)，则有：

其中，

步骤3：将观测值分为升轨los向、升轨方位向、降轨los向与降轨方位向4类，建立4类观测值的误差方程组，根据方差分量估计求解4类观测值单位权中误差的初步解。

简单起见可将各类观测值初始权均定为i，对于n个像素，则可得各类观测值权阵w1＝w2＝w3＝w4＝i，i为n×n的单位矩阵。

升轨los向观测值其对应的系数矩阵其中l1,b1分别为式(5)中l,b的第1,5,9…4n-3行数据的集合，同样也可得到观测值l2,l3,l4及相应系数矩阵b2,b3,b4。

令

则此时即可求得未知数初值t＝n^-1u；

进而可得

可得各类观测值改正数的二次型向量

则可得各类观测值单位权中误差向量

此时即可求得各类观测值单位权中误差初值。

步骤4：迭代计算各类观测值的精确权，求解高精度的三维地表形变。

在误差理论中各类观测值的单位权中误差近似相等时认为此时的权阵为最优，并且认为γ^-1δ是各类单位权中误差的无偏估计，由于初始假设各类观测值权相等均为1，所以在步骤3最后得到的各类观测值单位权中误差之间往往相差甚多，故本发明结合方差分量估计算法利用

进行权阵更新，其中代表更新后的权阵，并且假定第一类观测值(即升轨los向观测值)中误差所对应的权始终为单位权，即利用更新后的权阵重复步骤3直至各类观测值单位权中误差满足(即之间差别小于1mm²)方可停止迭代，此时中心点像素的高精度三维形变值即包含在步骤3解算的未知数t向量中。

针对不同像元重复上述步骤2-4，直至所有像素点的三维形变解算完成。

利用本发明所述方法对模拟数据进行三维地表形变求解，过程如下：

模拟数据描述：①在一定区域(图像大小400×450)模拟东西向、南北向及垂直向的三维形变场(如图3(a-c))；②假定升轨方位角349.24°，入射角38.75°，降轨方位角190.77°，入射角38.76°，根据成像几何投影关系，计算得到升降轨的los向与方位向投影形变；③在得到的升降轨的los向与方位向投影形变基础上添加高斯噪声，其标准差大小分别为5mm(升轨los向)、40mm(升轨方位向)、7mm(降轨los向)及35mm(降轨方位向)。此时即可得到用于模拟实验的原始数据，如图2所示，分别为升轨los向(图2(a))、升轨方位向(图2(b))、降轨los向(图2(c))及降轨方位向(图2(d))

传统求解三维地表形变场时是仅考虑地面单个点的观测数据，在假定形变场各态历经性的前提下，求解一定范围内(5×5像素区域)各类观测数据的方差以作为该类观测值的先验方差来进行定权。本次模拟实验分别利用传统解法(图3(g-i))与本发明提及的算法(图3(d-f))对模拟数据进行三维地表形变场的解算，通过求解出的三维形变场与原始模拟的三维形变场进行求差，得到三维地表形变场残差并统计其均方根(rms)，如表1：

表1三维地表形变场残差均方根

另外，在计算三维形变场的过程中，对两种方法得到的权重比例也进行了统计，统计数据如表2所示，具体对比数据图如图3所示：

表2不同方法得到的各类观测数据权重对比

在解算过程中，本发明默认升轨los向的权重为1。权重理论值是根据在模拟升降轨的los向与方位向数据所添加的高斯噪声标准差大小确定；在本发明算法及传统算法中，可得到400×450＝180000组权重比例，表2中的值是根据计算得出。

由表1、表2可知，本发明提及的算法相比于传统算法可得到更为精确权重比例及三维地表形变。

以上应用了具体个例对本发明进行阐述，只是为了帮助本领域中的普通技术人员很好的理解。在不偏离本发明的精神和范围的情况下，还可以对本发明的具体实施方式作各种推演、变形和替换。这些变更和替换都将落在本发明权利要求书所限定的范围内。