CN108446516B - 一种基于力源反演的构造断裂与地面沉降形变分解法 - Google Patents

Abstract

本发明公开了一种基于力源反演的构造断裂与地面沉降形变分解的方法，首先对研究区的SAR影像进行时序干涉处理，获取对应的年平均形变速率或形变时间序列；然后选取受地面沉降、构造断裂综合影响的形变区进行反演模型相关参数设置；进而对反演模型进行格网化，并通过格林函数与附加限制条件的线性最小二乘方法计算出各格网的压缩量，并结合格林函数正演出因地下水开采引起的地面沉降；最后将恢复所得的地面沉降形变从InSAR技术获取的地表形变中去除，得到构造断裂引起的地表形变。该方法对构造断裂形变与地面沉降形变综合出现、共同影响的研究区域的地质灾害的成灾模式以及防灾减灾具有指导意义。

Description

一种基于力源反演的构造断裂与地面沉降形变分解法

技术领域

本发明涉及地表形变监测技术领域，特别涉及一种基于力源反演的构造断裂与地面沉降形变分解法。

背景技术

在地球的很多区域同时受到构造断裂、岩浆活动、地下水开采、地下采矿、滑坡、地裂缝等不同组合的多源机制影响，表现在地表为多源形变叠加后的综合形变场，给人民的人身安全与财产安全带来了巨大的威胁。不同形变场信息的发生机制不同，地表形变表现不同，引起的后果不同，采取的防治措施也不同，因此对多源机制引起的综合形变场进行合理分解是一项兼具科学意义与工程价值的研究。如何有效地分离出不同尺度形变对区域地表形变的贡献越来越受到重视。

目前，采用空间大地测量技术进行地表形变监测时获取的地表形变大都是研究区的综合形变场，即不同尺度力源引起的形变综合作用的结果。绝大部分学者进行地表形变监测与分析时，忽略了影响量级较小但影响尺度较大的形变，仅考虑单一因素引起的形变(如地面沉降等)。常用的构造断裂形变监测方法主要为高精度的水准测量技术，该技术费时费力、覆盖范围有限，无法实现大范围高精度的形变监测。

发明内容

本发明的目的是提供一种基于力源反演的构造断裂与地面沉降形变分解法，利用地面沉降、构造断裂形变各自的机理特性，采用合适的地球物理模型进行地面沉降与构造断裂形变的分离，分析不同机制对地表形变的贡献，从而实现大范围高精度的形变监测。

本发明提供了一种基于力源反演的构造断裂与地面沉降形变分解的方法，该方法具体步骤如下：

首先对研究区的SAR影像进行时序干涉处理，获取对应的年平均形变速率或形变时间序列；然后选取受地面沉降、构造断裂综合影响的形变区进行反演模型的相关参数设置，进而对反演模型进行格网化，并通过格林函数与附加限制条件的线性最小二乘方法计算出各格网的压缩量；再采用格林函数正演出因地下水开采引起的地面沉降；最后将正演所得的地面沉降形变从InSAR技术获取的综合地表形变中扣除，得到构造断裂引起的地表形变；

该方法的具体步骤如下：

步骤1：模型参数的最优化确定，采用分布式矩形位错模型进行地表沉降形变反演；

矩形位错反演模型参数与地表形变观测数据之间满足如下公式：

d＝G(m)+ε (1)

其中，d表示地表形变观测数据资料，m表示反演模型参数，包括：矩形位错模型参考点的坐标(X,Y)、长度(Length)、宽度(Width)、走向(Strike)、倾角(Dip angle)、深度(Depth)；G(*)为联系观测数据与模型参数的非线性函数，ε表示观测数据的误差；在具体反演中，参考点的坐标、走向、长度和宽度可通过空间大地测量技术(如合成孔径雷达干涉测量(InSAR)技术、全球导航卫星系统(GNSS)技术)监测的地表形变来确定；深度参数可结合待研究区地下水开采深度给定，也可通过参数的最优化搜索确定；由于地下水开采引起的地面沉降主要表现为地层在垂直方向的压缩，故可将模型的倾向参数设置为0；

步骤2：分布式矩形压缩量的计算

首先对设置的矩形位错模型按照给定的格网大小进行格网化，各子块滑动量与地表形变观测值之间满足线性关系，对应的函数模型为：

d＝GS+ε (2)

式中，d为地表形变观测值向量，S为滑动参数向量，G为格林函数矩阵，其元素为各个矩形位错模型子块单位张量引起的地表点的位移，具体形式如下：

其中，上标代表矩形子块的编号，下标代表地表形变观测值个数；

然后，利用线性最小二乘方法反演各格网内的压缩量(滑动量)S：

S＝(G^TG)^-1G^Td (4)

步骤3：地下水引起的地表形变恢复

对步骤1确定的反演模型参数(矩形位错模型的长度(Length)、宽度(Width)、走向(Strike)、倾角(Dip angle)、深度(Depth))，利用公式(2)计算各高相干点的地表形变的模拟值，如公式(5)所示，地表任意一点在垂直方向的位移大小可表示为：

其中，(x,y)表示地表高相干点在局部坐标系中的坐标，L为矩形模型长度的一半，W为矩形模型的宽度，d为矩形模型的深度，R₁＝(x+L)²+y²+d²，R₂＝(x+L)²+(y-W)²+d²，R₃＝(x-L)²+y²+d²，R₄＝(x-L)²+(y-W)²+d²；

步骤4：构造断裂形变的获取

利用步骤3所计算的由地下水开采引起的地面沉降形变量，从空间大地测量技术获取的综合形变场中扣除，所得的残差形变即为构造断裂的形变，如公式(6)所示；

d_fault＝d_obs-u_z。 (6)

本发明的有益效果是：

本发明提供的一种基于力源反演的构造断裂与地面沉降形变分解法，有效地分离出不同尺度形变对地表形变的贡献比例，更好地为不同类型地质灾害的成灾机理分析提供基础信息。根据地面沉降、构造断裂形变的不同形变特征与力源机制，设置合理的弹性半空间矩形位错模型对因地下水开采引起的地面沉降形变进行恢复，然后从残余形变中获取构造断裂的形变信息，对于分析与揭示同时受到构造断裂、地下水开采等多源机制影响，构造断裂形变与地面沉降形变综合出现、共同影响的研究区域的地质灾害的成灾模式以及防灾减灾具有重要指导意义。

附图说明

图1为本发明提供的基于力源模型反演的不同尺度形变分解原理框图；

图2为山西清徐-交城区域地面沉降形变恢复；

图3为通过反演获取的矩形位错模型位移场；

图4为断裂构造形变。

具体实施方式

下面结合附图，对本发明的一个具体实施方式进行详细描述，但应当理解本发明的保护范围并不受具体实施方式的限制。

本发明具体实施方案以受地面沉降、构造断裂综合影响的山西清徐-交城区域为例来详细描述。

首先对研究区的SAR影像进行时序干涉处理，获取对应的年平均形变速率或形变时间序列(如附图2(a))；然后选取受地面沉降、构造断裂综合影响的形变区进行反演模型相关参数设置(长度为24km，宽度为7km，走向为66°，倾角为0°，深度为150m)；接下来就是对反演模型进行格网化(一般而言，格网越小计算效率越低，反演精度越高，反之亦然；本试验设置的格网为200m×200m)，并通过格林函数与附加限制条件的线性最小二乘方法计算出各格网的压缩量(如附图3)；在此基础上，结合格林函数正演出因地下水开采引起的地面沉降(如附图2(b))；最后将恢复所得的地面沉降形变从InSAR技术获取的地表形变中去除(如附图2(c))，便得到了构造断裂引起的地表形变(如附图4)。

为了验证本发明的应用效果，对扣除地面沉降形变后的残差形变进行精度统计，如附图2(d)。从附图2(c)和2(d)可以看出，残余形变的均方根误差很小，说明对于受地面沉降、构造断裂综合影响的区域，利用分布式矩形位错模型能够高精度地模拟出因地下水开采引起的地面沉降；附图4所示的构造断裂形变与其他监测结果具有很高的一致性，说明本发明是可行和可靠的。

综上所述，本发明实施例提供的一种基于力源反演的构造断裂与地面沉降形变分解法，有效地分离出不同尺度形变对地表形变的贡献量级，更好地为不同类型地质灾害的成灾机理分析提供基础信息。根据地面沉降、构造断裂形变的不同形变特征与力源机制，设置合理的弹性半空间矩形位错模型对因地下水开采引起的地面沉降形变进行恢复，然后从残余形变中获取构造断裂的形变信息，对于分析与揭示同时受到构造断裂、地下水开采等多源机制影响，构造断裂形变与地面沉降形变综合出现、共同影响的研究区域的地质灾害的成灾模式以及防灾减灾具有指导意义。

以上公开的仅为本发明的几个具体实施例，但是，本发明实施例并非局限于此，任何本领域的技术人员能思之的变化都应落入本发明的保护范围。

Claims (1)

1.一种基于力源反演的构造断裂与地面沉降形变分解法，其特征在于，

首先对研究区的SAR影像进行时序干涉处理，获取对应的年平均形变速率或形变时间序列；

然后选取受地面沉降、构造断裂综合影响的形变区进行反演模型的相关参数设置，进而对反演模型进行格网化，并通过格林函数与附加限制条件的线性最小二乘方法计算出各格网的压缩量；

再采用格林函数正演出因地下水开采引起的地面沉降；

最后将正演所得的地面沉降形变从InSAR技术获取的综合地表形变中扣除，得到构造断裂引起的地表形变；

该基于力源反演的构造断裂与地面沉降形变分解法的具体步骤如下：

步骤1：模型参数的最优化确定，采用分布式矩形位错模型进行地表沉降形变反演；

矩形位错反演模型参数与地表形变观测数据之间满足如下公式：

d＝G(m)+ε (1)

其中，d表示地表形变观测数据资料，m表示反演模型参数，包括：矩形位错模型参考点的坐标(X,Y)、长度、宽度、走向、倾角和深度；G(*)为联系观测数据与模型参数的非线性函数，ε表示观测数据的误差；在具体反演中，参考点的坐标、走向、长度和宽度可通过空间大地测量技术监测的地表形变来确定；深度参数可结合待研究区地下水开采深度给定，也可通过参数的最优化搜索确定；由于地下水开采引起的地面沉降主要表现为地层在垂直方向的压缩，故可将模型的倾向参数设置为0；

步骤2：分布式矩形压缩量的计算

首先对设置的矩形位错模型按照给定的格网大小进行格网化，各子块滑动量与地表形变观测值之间满足线性关系，对应的函数模型为：

d＝GS+ε (2)

式中，d为地表形变观测值向量，S为滑动参数向量，G为格林函数矩阵，其元素为各个矩形位错模型子块单位张量引起的地表点的位移，具体形式如下：

其中，上标代表矩形子块的编号，下标代表地表形变观测值个数；

然后，利用线性最小二乘方法反演各格网内的压缩量S：

S＝(G^TG)^-1G^Td (4)

步骤3：地下水引起的地表形变恢复

对步骤1确定的反演模型参数：矩形位错模型的长度、宽度、走向、倾角和深度，利用公式(2)计算各高相干点的地表形变的模拟值，如公式(5)所示，地表任意一点在垂直方向的位移大小可表示为：

其中，(x,y)表示地表高相干点在局部坐标系中的坐标，L为矩形模型长度的一半，W为矩形模型的宽度，d为矩形模型的深度，R₁＝(x+L)²+y²+d²，R₂＝(x+L)²+(y-W)²+d²，R₃＝(x-L)²+y²+d²，R₄＝(x-L)²+(y-W)²+d²；

步骤4：构造断裂形变的获取

利用步骤3所计算的由地下水开采引起的地面沉降形变量，从空间大地测量技术获取的综合形变场中扣除，所得的残差形变即为构造断裂的形变，如公式(6)所示；

d_fault＝d_obs-u_z (6)。

Images