留言板

尊敬的读者、作者、审稿人, 关于本刊的投稿、审稿、编辑和出版的任何问题, 您可以本页添加留言。我们将尽快给您答复。谢谢您的支持!

姓名
邮箱
手机号码
标题
留言内容
验证码

一种基于多面函数的改进最小二乘配置方法

谢曦霖 许才军 温扬茂 周力璇

谢曦霖, 许才军, 温扬茂, 周力璇. 一种基于多面函数的改进最小二乘配置方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(4): 592-598. doi: 10.13203/j.whugis20150664
引用本文: 谢曦霖, 许才军, 温扬茂, 周力璇. 一种基于多面函数的改进最小二乘配置方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(4): 592-598. doi: 10.13203/j.whugis20150664
XIE Xilin, XU Caijun, WEN Yangmao, ZHOU Lixuan. A Refined Least Squares Collocation Method Based on Multiquadric Function[J]. Geomatics and Information Science of Wuhan University, 2018, 43(4): 592-598. doi: 10.13203/j.whugis20150664
Citation: XIE Xilin, XU Caijun, WEN Yangmao, ZHOU Lixuan. A Refined Least Squares Collocation Method Based on Multiquadric Function[J]. Geomatics and Information Science of Wuhan University, 2018, 43(4): 592-598. doi: 10.13203/j.whugis20150664

一种基于多面函数的改进最小二乘配置方法

doi: 10.13203/j.whugis20150664
基金项目: 

国家高技术研究发展计划 2013AA122501-3

地震行业科研专项基金 201308009

国家自然科学基金 41431069

国家自然科学基金 41574002

国家自然科学基金 41721003

国家自然科学基金 41774011

详细信息
    作者简介:

    谢曦霖, 硕士, 主要从事地球物理大地测量的理论与方法研究。xlxie@whu.edu.cn

    通讯作者: 温扬茂, 博士, 副教授。ymwen@sgg.whu.edu.cn
  • 中图分类号: P223;P231

A Refined Least Squares Collocation Method Based on Multiquadric Function

Funds: 

The National High-Tech Research and Development Program of China (863 Program) 2013AA122501-3

the Special Earthquake Industry Research Project 201308009

the National Natural Science Foundation of China 41431069

the National Natural Science Foundation of China 41574002

the National Natural Science Foundation of China 41721003

the National Natural Science Foundation of China 41774011

More Information
    Author Bio:

    XIE Xilin, master, specializes in the theories and methods of Geophysical Geodesy.xlxie@whu.edu.cn

    Corresponding author: WEN Yangmao, PhD, associate professor.E-mail:ymwen@sgg.whu.edu.cn
图(4) / 表(2)
计量
  • 文章访问数:  1213
  • HTML全文浏览量:  102
  • PDF下载量:  347
  • 被引次数: 0
出版历程
  • 收稿日期:  2016-03-06
  • 刊出日期:  2018-04-05

一种基于多面函数的改进最小二乘配置方法

doi: 10.13203/j.whugis20150664
    基金项目:

    国家高技术研究发展计划 2013AA122501-3

    地震行业科研专项基金 201308009

    国家自然科学基金 41431069

    国家自然科学基金 41574002

    国家自然科学基金 41721003

    国家自然科学基金 41774011

    作者简介:

    谢曦霖, 硕士, 主要从事地球物理大地测量的理论与方法研究。xlxie@whu.edu.cn

    通讯作者: 温扬茂, 博士, 副教授。ymwen@sgg.whu.edu.cn
  • 中图分类号: P223;P231

摘要: 利用最小二乘配置对非平稳空间随机场进行推估时,趋势项数学模型的选择通常无法完整体现非平稳空间随机场的系统性,这将导致经验协方差函数估计出现偏差,最终推估结果可能错误。提出了一种基于多面函数的改进最小二乘配置方法来解决上述问题。该方法引入多面函数拟合区域内的趋势项,通过多次迭代计算得到稳定的待定系数值与协方差函数的参数值,最后综合趋势项与信号项得到最终估值。分别采用了模拟地震垂直形变数据和2009年意大利L’Aquila地震的合成孔径雷达干涉测量(Interferometric SAR,InSAR)与GPS同震位移数据来对该方法进行验证,并将其结果与常规方法进行比较。结果表明,改进方法在外部检核点估值的均方残差要小于多面函数法与常规的最小二乘配置法,且受采样点位的影响最小。

English Abstract

谢曦霖, 许才军, 温扬茂, 周力璇. 一种基于多面函数的改进最小二乘配置方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(4): 592-598. doi: 10.13203/j.whugis20150664
引用本文: 谢曦霖, 许才军, 温扬茂, 周力璇. 一种基于多面函数的改进最小二乘配置方法[J]. 武汉大学学报 ● 信息科学版, 2018, 43(4): 592-598. doi: 10.13203/j.whugis20150664
XIE Xilin, XU Caijun, WEN Yangmao, ZHOU Lixuan. A Refined Least Squares Collocation Method Based on Multiquadric Function[J]. Geomatics and Information Science of Wuhan University, 2018, 43(4): 592-598. doi: 10.13203/j.whugis20150664
Citation: XIE Xilin, XU Caijun, WEN Yangmao, ZHOU Lixuan. A Refined Least Squares Collocation Method Based on Multiquadric Function[J]. Geomatics and Information Science of Wuhan University, 2018, 43(4): 592-598. doi: 10.13203/j.whugis20150664
  • 最小二乘配置(least squares collocation)作为一种广义估计方法,已广泛应用于大地测量以及地球物理中空间随机场的插值问题中,如重力异常的插值[1]以及地壳形变场的估计[2-5]。平稳性是标准最小二乘配置的前提条件[6],然而实际中很多空间随机场并不满足这一假设,即出现非平稳的情况,如一些特殊地理位置的重力异常[7]以及地壳形变场[8-9]。空间随机场的非平稳表示其均值在空间上不一定为常值或(和)协方差是各项异性的且与空间位置相关,即所谓的二阶非平稳[10]

    大地测量中对非平稳场进行最佳线性无偏估计的方法最常用的为趋势去除法(trend removal),即先将数据中的趋势项去除,进行拟合推估后将趋势项加入到信号中得到最终结果。不同学者对于趋势项表达形式的选择有所区别:Tscherning[11]简单地将数据的均值作为趋势项从数据中去除,Goad等[12]使用线性多项式来去除布格重力异常引起的趋势项,Stopar等[13]使用人工神经网络来获取趋势面。趋势去除法目前有一些不足:①不能完全保证去除趋势项后的数据满足0均值的条件;②趋势项去除本身可能会引入系统误差[14]

    多面函数法最早由Hardy[15-16]提出,此后Franke[17]研究了29种不同的插值算法在离散数据插值问题中的应用,结果表明Hardy提出的多面函数法能得到准确的结果。Kansa[18]给出了多面函数效果好的解释。国内学者[19-20]利用多面函数法对地壳垂直形变场拟合做了大量研究。由于多面函数法是一种几何上的逼近,不需要推估点之间的统计信息,得到的是区域内趋势性和系统性的连续场,对于区域内随机的、不规则的变化则无法很好地逼近。

    为了削弱应用最小二乘配置的趋势去除法对非平稳空间随机场进行推估的不足,本文在最小二乘配置的框架下引入多面函数法对趋势项进行拟合,以期更好地描述空间随机场的系统性。由于非平稳空间随机场中包含有趋势项与随机项,本文通过迭代的方法得到信号稳定的协方差函数参数值和多面函数的待定参数值,最后综合区域的趋势项与信号项得到最终的推估值。

    • 由Hardy[15-16]提出的基本方案可知,任意函数f可以写成N个连续可微的基础函数的组合:

      $$ f\left( {x,y} \right) = \sum\limits_{j = 1}^N {{a_j}{g_j}\left( {x,y,{x_j},{y_j}} \right)} $$ (1)

      式中,gj为核函数;αj为待定系数;N为核函数个数。核函数一般可选为:

      $$ {g_j}\left( {x,y,{x_j},{y_j}} \right) = {\left[ {{{\left( {x - {x_j}} \right)}^2} + {{\left( {y - {y_j}} \right)}^2} + \delta } \right]^k} $$ (2)

      式中,(xj, yj)为节点坐标;δ为光滑因子;k决定了核函数的形状。陶本藻等[19]研究表明,在推估中k值可取为0.5,即采用正双曲面型,δ取0。在推估中,将研究区域视为一数学表面,根据多面函数的原理可利用多个核函数的组合来逼近该数学表面。利用已知离散点数据计算得到所有待定系数αj的估值,即可利用式(1)推估得到所有未测点的估值。

      标准的最小二乘配置函数模型为:

      $$ \mathit{\boldsymbol{l}} = \mathit{\boldsymbol{Bx}} + {\mathit{\boldsymbol{G}}_\mathit{\boldsymbol{y}}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} $$ (3)

      式中,B为设计矩阵; l为观测向量;Δ为观测噪声向量;Gy为倾向部分; x为随机信号向量。最小二乘配置的原理是通过已测点数据和已测点信号与未知点信号之间的协方差来推估未知点的值,其核心问题在于已测点信号与未知点信号之间协方差的确定。式(3)的最佳线性无偏解为[6]

      $$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\hat y}} = {\left\{ {{\mathit{\boldsymbol{G}}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}}} \right)}^{ - 1}}\mathit{\boldsymbol{G}}} \right\}^{ - 1}} \cdot \\ \;\;\;\;\;\;{\mathit{\boldsymbol{G}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{l}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{l}}_\mathit{\boldsymbol{x}}}} \right)\\ \mathit{\boldsymbol{\hat x}} = {\mathit{\boldsymbol{l}}_\mathit{\boldsymbol{x}}} + {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{l}} - {\mathit{\boldsymbol{G}}_{\mathit{\boldsymbol{\hat y}}}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{l}}_\mathit{\boldsymbol{x}}}} \right)\\ \mathit{\boldsymbol{\hat x'}} = {\mathit{\boldsymbol{l}}_{\mathit{\boldsymbol{x'}}}} + {\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{x'x}}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{l}} - {\mathit{\boldsymbol{G}}_{\mathit{\boldsymbol{\hat y}}}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{l}}_\mathit{\boldsymbol{x}}}} \right) \end{array} \right. $$ (4)

      式中,lx为信号的均值,通常取lx = 0Dx为已测点信号之间的协方差矩阵;Dx'x为推估点信号与已测点信号之间的协方差阵;$\mathit{\boldsymbol{\hat x}}{\rm{'}}$为待推估的信号值。

      在最小二乘配置的框架下,本文将多面函数法融入到最小二乘配置的数学模型中,即用多面函数法对区域的倾向项进行估计,然后利用拟合推估方法对空间随机变化项进行补偿,综合两种方法的优点以达到更精确的推估效果。改进最小二乘配置的函数模型可以表示为:

      $$ \mathit{\boldsymbol{l}} = \mathit{\boldsymbol{Bx}} + \mathit{\boldsymbol{Qa}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} $$ (5)

      式中,ln×1观测向量;Qn×m核函数阵,m为节点数目;am×1待定系数向量;Δn×1误差向量。与式(4)同理, 得到待定系数a以及信号x的最佳线性无偏解:

      $$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\hat a}} = {\left\{ {{\mathit{\boldsymbol{Q}}^{\rm{T}}}{{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}}} \right)}^{ - 1}}\mathit{\boldsymbol{Q}}} \right\}^{ - 1}} \cdot \\ \;\;\;\;\;\;\;{\mathit{\boldsymbol{Q}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}}} \right)^{ - 1}}\left( {\mathit{\boldsymbol{l}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{l}}_\mathit{\boldsymbol{x}}}} \right)\\ \mathit{\boldsymbol{\hat x}} = {\mathit{\boldsymbol{l}}_\mathit{\boldsymbol{x}}} + {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}}} \right)^{ - 1}} \cdot \\ \;\;\;\;\;\;\;\left( {\mathit{\boldsymbol{l}} - \mathit{\boldsymbol{Q\hat a}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{l}}_\mathit{\boldsymbol{x}}}} \right)\\ \mathit{\boldsymbol{\hat x'}} = {\mathit{\boldsymbol{l}}_{\mathit{\boldsymbol{x'}}}} + {\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{x'x}}}}{\mathit{\boldsymbol{B}}^{\rm{T}}}{\left( {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{B}}^{\rm{T}}} + {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{ \boldsymbol{\varDelta} }}}} \right)^{ - 1}} \cdot \\ \;\;\;\;\;\;\;\left( {\mathit{\boldsymbol{l}} - \mathit{\boldsymbol{Q\hat a}} - \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{l}}_\mathit{\boldsymbol{x}}}} \right) \end{array} \right. $$ (6)

      式中,lx为信号的均值,通常取lx =0;Dx为已测点信号之间的协方差矩阵;Dx'x为推估点信号与已测点信号之间的协方差阵;$\mathit{\boldsymbol{\hat x}}{\rm{'}}$为待推估的信号值。

      由式(4)可知,多面函数待定参数的估计需要已知的信号协方差阵,而信号的协方差阵则需要首先从观测值中将趋势项去除才能统计拟合得到。本文使用陈家军等[21]提出的多次迭代的方法解决这一矛盾,获取稳健的经验协方差函数和多面函数待定系数值。对于地表垂直形变,可以近似认为是各项同性的,用高斯模型来拟合经验协方差值[1]。高斯模型的表达式如下:

      $$ C\left( d \right) = C\left( 0 \right){e^{ - {\gamma ^2}{d^2}}} $$ (7)

      式中,d为两点间的距离;C(0)与γ为待定系数, C(0)表示信号的方差, γ为信号的空间相关尺度。改进方法具体的计算步骤如下:

      1) 初始化数据协方差阵为单位阵(Dx= E),即C(0)i=1,γi=0,i=0;

      2) 利用协方差函数模型参数C(0)iγi得到协方差阵Dx。由于BQ已知,通过式(6)得到待定系数的估值ai

      3) 根据式(5)去除观测数据l中的趋势项Qai,得到信号xi

      4) 根据文献[1]中的方法对信号xi进行实验协方差统计,得到不同距离段的协方差值(d, C(d)),继而利用最小二乘拟合得到协方差函数模型参数C(0)i+1γi+1,更新协方差阵Dx

      5) 将更新协方差阵Dx代入式(6)得到新的待定系数值ai+1

      6) 判断条件:|C(0)i+1-C(0)i| < 10-6且|γi+1-γi| < 10-6,不满足则i=i+1,回到第2步。否则继续;

      7) 根据待定系数ai+1,利用式(6)得到未测点的信号值$\mathit{\boldsymbol{\hat x}}{\rm{'}}$,最后综合未测点的趋势项与信号值得到最终估值:$\mathit{\boldsymbol{\hat l}}' = \mathit{\boldsymbol{\hat x}}'\mathit{ + }\mathit{\boldsymbol{Q}}{\mathit{\boldsymbol{a}}_{i + 1}}$。

    • 为了验证改进方法的有效性,本文首先通过一个模拟的地震滑动数据来比较改进方法与常规方法的推估效果。考虑一个位于均匀的,各项同性的弹性半空间中的倾向断层,其长度为75 km,宽为60 km,走向角为70°,倾角为15°,模拟的滑动分布见图 1

      图  1  模拟的断层滑动分布图

      Figure 1.  Fault Geometry and Slip Distribution Used in the Synthetic Test

      利用Okada公式[22]计算模拟断层滑动引起的覆盖研究区域的374个均匀格网点(格网点XY方向间隔均为4 km)的地表垂直形变,分布模式见图 2。此外,分别按1 km和2 km进行采样得到的结果与4 km采样的结果类似,鉴于篇幅所限,文中主要展示了点位最为稀疏的采样(即4 km)。

      图  2  模拟垂直形变值与数据点位

      Figure 2.  Vertical Surface Displacement and

      从生成的374个已知形变值的点中随机选择200个点作为计算点,另外174个点作为外部检核点。对计算点加入均值为0 mm,标准差为3 mm的高斯误差。多面函数计算中,选取了全部计算点数作为节点数,且均匀覆盖研究区域。为了进行比较,使用趋势项为线性多项式的最小二乘配置与多面函数法进行了计算。实验中对计算点的随机选取共模拟了500次,几种不同方法的外部检核残差统计见表 1。通过对计算结果以及拟合残差均方根(RMS)统计结果(表 1)的比较分析可以发现:多面函数法得到的结果外部检核点的RMS为最大,以线性多项式为趋势的常规最小二乘配置的RMS次之,改进方法的RMS是最小的,表明改进方法得到的结果最为精确。通过比较这几种方法RMS的标准差可知,改进方法受计算点选取的影响最小,得到的结果更加稳定。

      表 1  基于500次随机实验的不同方法外部检核点RMS统计结果/mm

      Table 1.  External Check Points RMS Statistics Based on 500 Random Experiments/mm

      推估方法 最大值 最小值 均值 标准差
      多面函数 17.44 7.51 11.04 1.42
      常规最小二乘配置 12.91 5.38 8.31 1.11
      改进方法 10.04 4.08 6.79 0.88

      为了评估不同噪声水平对推估结果的影响,本文还对模拟数据加入均值为0 mm,标准差分别为5 mm和10 mm的高斯噪声进行试验。当加入标准差为5 mm的高斯噪声后,多面函数法得到RMS均值为13.9 mm,常规最小二乘配置为8.7 mm,改进方法为7.5 mm。当加入标准差为10 mm的高斯噪声后,多面函数为22.8 mm,常规最小二乘配置为10.5 mm,改进方法为10.1 mm。随着噪声方差值的增加,多面函数得到的结果RMS增加最迅速,而改进方法与常规最小二乘配置得到的RMS趋于一致。这是因为随着外部噪声水平的增加,当噪声的数值与信号数值相当时(图 3),无法有效地将信号从噪声中提取出来,因此改进方法无法发挥其优势。

      图  3  噪声均值为0 mm, 标准差为3 mm时利用改进算法得到的信号项

      Figure 3.  Signal Distribution Estimated by Refined Least Squares Collocation with 0 mm Mean and 3 mm Standard Deviation Gaussian Error and Trend Distribution Estimated

    • 2009年4月6日,意大利中部Abruzzo地区的L’Aquila发生Mw 6.3级地震, 致使300人死亡。为了对改进方法的效果进行检验,本文使用了意大利L’Aquila地震InSAR和GPS的同震形变数据来进行计算分析。Envisat卫星Track 129的干涉经降采样得到的1 282个点沿视线向的形变值来自文献[23],56个GPS水平形变值与26个GPS垂直形变值来自文献[24],点位分布见图 4。InSAR点视线向形变中误差为6.8 mm,GPS水平形变中误差为3.9 mm,GPS垂直形变中误差为7.8 mm。由于InSAR数据得到的是沿视线向的形变值,并不能得到垂直方向的形变。通过InSAR成像几何可以得到视线向形变值与东、北和垂直向形变值的转换关系:

      图  4  L’Aquila地震震中区的InSAR点与垂直GPS点分布图

      Figure 4.  Distribution of InSAR and Vertical GPS Observations

      $$ \begin{array}{l} {d_{{\rm{los}}}} = \\ \begin{array}{*{20}{c}} {{{\left[ { - \sin \theta \sin \left( {{\alpha _h} + \frac{{\rm{ \mathsf{ π} }}}{2}} \right) - \sin \theta \cos \left( {{\alpha _h} + \frac{{\rm{ \mathsf{ π} }}}{2}} \right)\cos \theta } \right]}^{\rm{T}}} \cdot }\\ {\left[ {\begin{array}{*{20}{c}} {{d_e}}&{{d_n}}&{{d_u}} \end{array}} \right]} \end{array} \end{array} $$ (8)

      式中,[de dn du]为地面点在东、北和垂直向的形变量;αh为雷达飞行方向的方位角;θ为雷达的侧视角。等式右边第一项表示视线向形变值方向单位矢量在东、北和垂直向的分量。利用GPS东、北向离散形变值插值得到所有InSAR点东、北方向形变值,再利用式(8)即可将雷达视线向形变dlos转换到垂直向形变du

      为了检核推估效果,本文将1 282个InSAR数据点作为已测点,将26个GPS点垂直向形变值作为外部检核。计算中,由于InSAR点数量较多,且空间分布较为接近,导致信号的协方差阵的复共线性严重,从而导致配置解的不稳定。为了保证得到稳定的结果,需要对其计算模型进行正则化修正。由于推估点与计算点是同类型数据,即设计矩阵为单位阵B = E。根据文献[25]的方法可将式(6)改写为:

      $$ \left\{ \begin{array}{l} \mathit{\boldsymbol{\hat a}} = {\left\{ {{\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{K}}^{ - 1}}\mathit{\boldsymbol{Q}}} \right\}^{ - 1}}{\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{K}}^{ - 1}}\mathit{\boldsymbol{l}}\\ \mathit{\boldsymbol{\hat x}} = {\mathit{\boldsymbol{D}}_\mathit{\boldsymbol{x}}}{\mathit{\boldsymbol{K}}^{ - 1}}\left( {\mathit{\boldsymbol{l}} - \mathit{\boldsymbol{Q\hat a}}} \right)\\ \mathit{\boldsymbol{\hat x'}} = {\mathit{\boldsymbol{D}}_{\mathit{\boldsymbol{x'x}}}}{\mathit{\boldsymbol{K}}^{ - 1}}\left( {\mathit{\boldsymbol{l}} - \mathit{\boldsymbol{Q\hat a}}} \right) \end{array} \right. $$ (9)

      式中,K = Dx + DΔ。令:

      $$ \mathit{\boldsymbol{q}} = {\mathit{\boldsymbol{K}}^{ - 1}}\mathit{\boldsymbol{l}} $$ (10)
      $$ \mathit{\boldsymbol{G}} = {\left\{ {{\mathit{\boldsymbol{Q}}^{\rm{T}}}{\mathit{\boldsymbol{K}}^{ - 1}}\mathit{\boldsymbol{Q}}} \right\}^{ - 1}}{\mathit{\boldsymbol{Q}}^{\rm{T}}} $$ (11)

      假设G为稳定的,则只需要得到稳定的q,即可确保可以获得稳定的待估量${\mathit{\boldsymbol{\hat a}}}$。将式(10)改写为:

      $$ \mathit{\boldsymbol{Kq}} = \mathit{\boldsymbol{l}} $$ (12)

      对式(12)引入Tikhonov正则化法,其正则化准则为[26]

      $$ {F_\alpha }\left( \mathit{\boldsymbol{q}} \right) = {\left\| {\mathit{\boldsymbol{Kq}} - \mathit{\boldsymbol{l}}} \right\|^2} + \alpha {\left\| \mathit{\boldsymbol{q}} \right\|^2} = \min $$ (13)

      式中,α为正则化参数;‖qq‖范数。通过式(13)可得q的正则化解:

      $$ {{\mathit{\boldsymbol{\hat q}}}_\alpha } = {\left( {{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{K}} + \alpha \mathit{\boldsymbol{I}}} \right)^{ - 1}}{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{l}} $$ (14)

      通过式(14)可知协方差阵K的逆的稳定形式可以用(KTK+αI)-1 K T表示。考虑G为不稳定的,令:

      $$ \begin{array}{*{20}{c}} {\mathit{\boldsymbol{J}} = {{\left( {{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{K}} + \alpha \mathit{\boldsymbol{I}}} \right)}^{ - 1}}{\mathit{\boldsymbol{K}}^{\rm{T}}}}\\ {\mathit{\boldsymbol{H}} = {\mathit{\boldsymbol{Q}}^{\rm{T}}}\mathit{\boldsymbol{JQ}}} \end{array} $$ (15)

      则式(9)可改写为:

      $$ \mathit{\boldsymbol{\hat a}} = {\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{Q}}^{\rm{T}}}\mathit{\boldsymbol{Jl}} $$ (16)

      同理,要得到稳定的待估量${\mathit{\boldsymbol{\hat a}}}$,则必须得到稳定的H。通过相似的推导可以得到${\mathit{\boldsymbol{\hat a}}}$的正则化解:

      $$ {{\mathit{\boldsymbol{\hat a}}}_{\alpha ,\beta }} = {\left( {{\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{K}} + \beta \mathit{\boldsymbol{I}}} \right)^{ - 1}}{\mathit{\boldsymbol{H}}^{\rm{T}}}{\mathit{\boldsymbol{Q}}^{\rm{T}}}\mathit{\boldsymbol{Jl}} $$ (17)

      式中,β为正则化参数。

      利用上述方法进行正则化修正得到信号的正则化解为:

      $$ \hat x = {\mathit{\boldsymbol{D}}_x}\mathit{\boldsymbol{J}}\left( {\mathit{\boldsymbol{l}} - \mathit{\boldsymbol{Q}}\hat a} \right) $$ (18)
      $$ \hat x' = {\mathit{\boldsymbol{D}}_{x'x}}\mathit{\boldsymbol{J}}\left( {\mathit{\boldsymbol{l}} - \mathit{\boldsymbol{Q}}\hat a} \right) $$ (19)

      上文中正则化参数的估计可以采用Hansen[27]提出的L曲线法求得,其核心是定位L曲线上曲率最大的那个点的正则化参数值。改进方法的推估流程与§2中的流程一致,只是用式(17)~式(19)替代式(6)。

      为了对计算结果进行比较分析,文中还同时使用了正则化修正的多面函数和以线性多项式为趋势项的最小二乘配置来计算。不同方法的外部检核残差统计结果(表 2)表明, 改进方法相对于另外两种常规方法,各项指标都是最优的。实测数据计算得到了与模拟数据计算相同的结论,这表明将多面函数引入后能够更真实地反映出该区域系统性的形变,将系统部分去除得到的随机场也更加平稳,因此能够通过推估的方法得到准确的信号值。但是由于InSAR数据处理中,视线向形变值残留了轨道误差、DEM误差以及大气误差等的影响,而且GPS的垂直向精度也不是很高,这会导致两者之间残差RMS较大。

      表 2  各种方法的GPS外部检核点残差统计结果/mm

      Table 2.  External GPS Check Points Residuals Statistics/mm

      推估方法 残差最小值 残差最大值 残差均值 RMS
      多面函数 -40.7 39.1 -0.17 23.4
      最小二乘配置 -39.3 30.3 0.02 17.2
      改进方法 -38.6 25.1 -0.01 14.2
    • 在利用最小二乘配置对非平稳空间场进行推估时,趋势去除法是最为广泛的策略,但是在实际应用中由于趋势项数学模型选择的不准确导致无法完全去除区域的趋势。本文在分析多面函数法的优良数学特性后,提出了一种改进的最小二乘配置方法来解决上述问题。该方法利用多面函数法来准确拟合区域的倾向项,利用迭代的方法获得稳定的待定系数值与协方差函数系数值,最后综合趋势项与信号项得到了待估点的估值。

      本文首先通过对一组模拟数据进行计算验证了改进方法的有效性。改进方法得到的结果要优于多面函数法以及以线性多项式为趋势的最小二乘配置法。此外,还利用该方法对2009年意大利L’Aquila地震的同震形变数据进行计算。针对InSAR数据量多,且分布密集的特点,引入了正则化方法对改进方法解的表达式进行修正以解决病态对结果带来的影响。文中使用了多面函数与常规最小二乘配置法作为对比,通过外部检核GPS点残差的统计结果表明,改进方法同样能够得到最佳的推估值。

      本文所使用的协方差函数模型只考虑了各项同性的情况,如果区域内信号有明显的各向异性,那么在进行推估时需要引入各向异性的协方差函数模型。由于改进方法是基于多面函数的,那么在利用多面函数进行推估时遇到的问题,改进方法也同样会遇到,即节点数量以及位置的选择、平滑因子的选择和核函数形状的选择。至于用多面函数法进行趋势拟合后对信号的影响,有待于更进一步研究。

参考文献 (27)

目录

    /

    返回文章
    返回