-
最小二乘配置(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 = 0;Dx为已测点信号之间的协方差矩阵;Dx'x为推估点信号与已测点信号之间的协方差阵;$\mathit{\boldsymbol{\hat x}}{\rm{'}}$为待推估的信号值。
在最小二乘配置的框架下,本文将多面函数法融入到最小二乘配置的数学模型中,即用多面函数法对区域的倾向项进行估计,然后利用拟合推估方法对空间随机变化项进行补偿,综合两种方法的优点以达到更精确的推估效果。改进最小二乘配置的函数模型可以表示为:
$$ \mathit{\boldsymbol{l}} = \mathit{\boldsymbol{Bx}} + \mathit{\boldsymbol{Qa}} + \mathit{\boldsymbol{ \boldsymbol{\varDelta} }} $$ (5) 式中,l为n×1观测向量;Q为n×m核函数阵,m为节点数目;a为m×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。由于B、Q已知,通过式(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。
利用Okada公式[22]计算模拟断层滑动引起的覆盖研究区域的374个均匀格网点(格网点X与Y方向间隔均为4 km)的地表垂直形变,分布模式见图 2。此外,分别按1 km和2 km进行采样得到的结果与4 km采样的结果类似,鉴于篇幅所限,文中主要展示了点位最为稀疏的采样(即4 km)。
从生成的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),无法有效地将信号从噪声中提取出来,因此改进方法无法发挥其优势。
-
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成像几何可以得到视线向形变值与东、北和垂直向形变值的转换关系:
$$ \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) 式中,α为正则化参数;‖q为q‖范数。通过式(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点残差的统计结果表明,改进方法同样能够得到最佳的推估值。
本文所使用的协方差函数模型只考虑了各项同性的情况,如果区域内信号有明显的各向异性,那么在进行推估时需要引入各向异性的协方差函数模型。由于改进方法是基于多面函数的,那么在利用多面函数进行推估时遇到的问题,改进方法也同样会遇到,即节点数量以及位置的选择、平滑因子的选择和核函数形状的选择。至于用多面函数法进行趋势拟合后对信号的影响,有待于更进一步研究。
-
摘要: 利用最小二乘配置对非平稳空间随机场进行推估时,趋势项数学模型的选择通常无法完整体现非平稳空间随机场的系统性,这将导致经验协方差函数估计出现偏差,最终推估结果可能错误。提出了一种基于多面函数的改进最小二乘配置方法来解决上述问题。该方法引入多面函数拟合区域内的趋势项,通过多次迭代计算得到稳定的待定系数值与协方差函数的参数值,最后综合趋势项与信号项得到最终估值。分别采用了模拟地震垂直形变数据和2009年意大利L’Aquila地震的合成孔径雷达干涉测量(Interferometric SAR,InSAR)与GPS同震位移数据来对该方法进行验证,并将其结果与常规方法进行比较。结果表明,改进方法在外部检核点估值的均方残差要小于多面函数法与常规的最小二乘配置法,且受采样点位的影响最小。Abstract: Trend removal is the most common approach in conventional least squares collocation (LSC) to deal with nonstationarity. Due to the inaccuracy of the trend model, conventional LSC can barely eliminate the drift component of the data that results in estimation bias of the local covariance function and error of the interpolation values. Here we present a refined LSC method to estimate the drift component of the field and compensate the residualfrom LSC. The refined method employs the mulitiquadric function to approach the drift and the collocation to estimate signals. We apply the refined method to a synthetic data set and coseismic displacement data from the 2009 L'Aquila, Italy earthquake, and compare the results of refined method with conventional methods. The statistical results of residuals indicate that the refined method can achieve a more accurate result than conventional methods and is affected less by the observation distribution.
-
Key words:
- multiquadric function /
- least squares collocation /
- spatial interpolation /
- ill-condition
-
表 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 表 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 -
[1] Moritz H. Least Squares Collocation[J].Reviews of Geophysics and Space Physics, 1978(16):421-430 doi: 10.1029/RG016i003p00421/full [2] 吴啸龙, 杨志强, 党永超.基于球面最小二乘配置的福建省地壳水平形变研究[J].武汉大学学报·信息科学版, 2015, 40(3):401-405 http://ch.whu.edu.cn/CN/abstract/abstract3218.shtml Wu Xiaolong, Yang Zhiqiang, Dang Yongchao. Study on Horizontal Crustal Deformation in Fujian Province Using Spherical Least Square Collocation[J].Geomatics and Information Science of Wuhan University, 2015, 40(3):401-405 http://ch.whu.edu.cn/CN/abstract/abstract3218.shtml [3] 陶本藻, 姚宜斌.基于多面核函数配置型模型的参数估计[J].武汉大学学报·信息科学版, 2003, 28(5):547-550 http://ch.whu.edu.cn/CN/abstract/abstract4778.shtml Tao Benzao, Yao Yibin. Parameter Estimation Based on Multi-quadric Collocation Model[J]. Geo-matics and Information Science of Wuhan University, 2003, 28(5):547-550 http://ch.whu.edu.cn/CN/abstract/abstract4778.shtml [4] Kahle H G, Cocard M, Peter Y, et al. The GPS Strain Field in the Aegean Sea and Western Anatolia[J].Geophysical Research Letters, 1999, 26(16):2513-2516 doi: 10.1029/1999GL900403 [5] Kahle H G, Cocard M, Peter Y, et al. GPS-Derived Strain Rate Field Within the Boundary Zones of the Eurasian, African, and Arabian Plates[J].Journal of Geophysical Research, 2000, 105(B10):23353-23370 doi: 10.1029/2000JB900238 [6] Moritz H. Advanced Physical Geodesy[M]. Tunbridge Wells:Abacus, 1980 [7] Lambeck K. The Perth Basin:A Possible Framework for its Formation and Evolution[J].Exploration Geophysics, 1987, 18(1-2):124-128 doi: 10.1071/EG987124 [8] 柴洪洲, 崔岳, 明锋.最小二乘配置方法确定中国大陆主要块体运动模型[J].测绘学报, 2009, 38(1):61-65 https://www.wenkuxiazai.com/doc/e81c2c748e9951e79b892774.html Chai Hongzhou, Cui Yue, Ming Feng. The Determination of Chinese Mainland Crustal Movement Model Using Least-Squares Collocation[J].ActaGeodaetica et Cartographica Sinica, 2009, 38(1):61-65 https://www.wenkuxiazai.com/doc/e81c2c748e9951e79b892774.html [9] 陈光保, 陈永奇, 何秀凤.基于改进最小二乘配置的地壳垂直形变分析[J].大地测量与地球动力学, 2010, 30(4):128-132 http://www.cnki.com.cn/Article/CJFDTotal-CHTB201301008.htm Chen Guangbao, Chen Yongqi, He Xiufeng. Analysis of Vertical Crust Deformation by Using Improved Least Square Collocation[J].Journal of Geodesy and Geodynamics, 2010, 30(4):128-132 http://www.cnki.com.cn/Article/CJFDTotal-CHTB201301008.htm [10] Armstrong M. Basic Linear Geostatistics[M]. Berlin:Springer, 1998 [11] Tscherning C C. Geoid Determination by 3D Least-Squares Collocation[M]//Sanso F, Sideris M G. Geoid Determination; Theory and Methods, Lecture Notes in Earth System Sciences 110. Berlin Heidelberg:Springer-Verlag, 2013:311-329 [12] Goad C C, Tscherning C C, Chin M M. Gravity Empirical Covariance Values for the Continental United States[J]. Journal of Geophysical Research, 1984, 89(B9):7962-7968 doi: 10.1029/JB089iB09p07962 [13] Stopar B, Kuhar A T M, Turk G. GPS-Derived Geoid Using Artificial Neural Network and Least Squares Collocation[J]. Survey Review, 2006, 38(300):513-524 doi: 10.1179/sre.2006.38.300.513 [14] Cressie N. Statistics for Spatial Data[M]. New York:Wiley Inc., 1993 [15] Hardy R L. Multiquadric Equations of Topography and Other Irregular Surfaces[J]. Journal of Geophysical Research, 1971, 76(8):1905-1915 doi: 10.1029/JB076i008p01905 [16] Hardy R L. Research Results in the Application of Multiquadric Equations to Surveying and Mapping Problems[J]. Surveying and Mapping, 1975, 35:321-332 https://trid.trb.org/view/66619 [17] Franke R. Scattered Data Interpolation-Tests of Some Methods[J]. Mathematics of Computation, 1982, 38(157):181-200 http://www.ams.org/journals/mcom/1982-38-157/S0025-5718-1982-0637296-4/home.html [18] Kansa E J. Multiquadrics-A Scattered Data Approximation Scheme with Applications to Computational Fluid-Dynamics-I Surface Approximations and Partial Derivative Estimates[J]. Computers & Mathematics with Applications, 1990, 19(8-9):127-145 https://www.sciencedirect.com/science/article/pii/089812219090270T [19] 武艳强, 江在森, 杨国华, 等.利用多面函数整体求解GPS应变场的方法及应用[J].武汉大学学报·信息科学版, 2009, 34(9):1085-1089 http://ch.whu.edu.cn/CN/abstract/abstract1372.shtml Wu Yanqiang, Jiang Zaisen, Yang Guohua, et al. Application and Method of GPS Strain Calculating in Whole Mode Using Multi-surface Function[J].Geomatics and Information Science of Wuhan University, 2009, 34(9):1085-1089 http://ch.whu.edu.cn/CN/abstract/abstract1372.shtml [20] 刘经南, 施闯, 姚宜斌, 等.多面函数拟合法及其在建立中国地壳平面运动速度场模型中的应用研究[J].武汉大学学报·信息科学版, 2001, 26(6):500-508 http://ch.whu.edu.cn/CN/abstract/abstract5221.shtml Liu Jingnan, Shi Chuang, Yao Yibin, et al. Hardy Function Interpolation and its Applications to the Establishment of Crustal Movement Speed Field[J].Geomatics and Information Science of Wuhan University, 2001, 26(6):500-508 http://ch.whu.edu.cn/CN/abstract/abstract5221.shtml [21] 陈家军, 郭乔羽, 王红旗, 等.应用协同-泛克立格法估计地下水水位[J].水利学报, 2000, 31(7):7-13 https://www.wenkuxiazai.com/doc/ac9a6ec0360cba1aa911da05.html Chen Jiajun, Guo Qiaoyu, Wang Hongqi, et al. The Application of Co-universal Kriging Method for Estimation of Groundwater Table[J]. Journal of Hydraulic Engineering, 2000, 31(7):7-13 https://www.wenkuxiazai.com/doc/ac9a6ec0360cba1aa911da05.html [22] Okada Y. Surface Deformation due to Shear and tensile Faults in a Half-Space[J].International Journal of Rock Mechanics and Mining Sciences &Geomechanics Abstracts, 1985, 23(4):1135-1154 http://bssa.geoscienceworld.org/content/75/4/1135.short [23] 温扬茂, 何平, 许才军, 等.联合Envisat和ALOS卫星影像确定L'Aquila地震震源机制[J].地球物理学报, 2012, 55(1):53-65 Wen Yangmao, He Ping, Xu Caijun, et al. Source Parameters of the 2009 L'Aquila Earthquake, Italy from Envisat and ALOS Satellite SAR Images[J]. Chinese Journal of Geophysics, 2012, 55(1):53-65 [24] Cheloni D, D'Agostino N, D'Anastasio E, et al. Coseismic and Initial Post-Seismic Slip of the 2009 Mw 6.3 L'Aquila Earthquake, Italy, from GPS Measurements[J].Geophysical Journal International, 2010, 181(3):1539-1546 doi: 10.1111/j.1365-246X.2010.04584.x [25] 黄谟涛, 欧阳永忠, 翟国君, 等.融合多源重力数据的Tikhonov正则化配置法[J].海洋测绘, 2013, 33(3):6-12 http://d.wanfangdata.com.cn/Periodical_hych201303002.aspx Huang Motao, Ouyang Yongzhong, Zhai Guojun, et al. Tikhonov Regularization Collocation for Multi-source Gravity Data Fusion Processing[J]. Hydrographic Surveying and Charting, 2013, 33(3):6-12 http://d.wanfangdata.com.cn/Periodical_hych201303002.aspx [26] Tikhonov A N, Arsenin V Y. Solutions of Ill-posed Problems[M]. New York:Wiley Inc., 1977 [27] Hansen P C. Analysis of Discrete Ill-posed Problems by Means of the L-curve[J]. Siam Review, 1992, 34(4):561-580 doi: 10.1137/1034115 -