多核函数法是Hardy在1971年提出的一种数值逼近方法。其基本思想是,任何数学表面和任何不规则的圆滑表面,总可以通过一系列有规则的数学表面的合成,以任意精度逼近。其数学表达式为

f(x,y)=∑^m_j=1c_js_j(x,y,x_j,y_j).(1)

式中,s_j(x,y,x_j,y_j)为核函数,(x_j,y_j)为核点的坐标,c_j为待定系数。

核函数是一个非常简单的曲面形态,通常为“钵”形或“钟”形曲面函数。虽然哪一种形态的核函数更优尚无非常明确的定论,但在后来的应用研究中,该法已被广泛地用于地壳形变研究(黄立人,刘天奎,1989,杨国华,黄立人,1990)、DEM内插(周光文,黄莜容,1996)和地球重力场元内插(陶本藻,1992),并被证明是非常成功和有效的。然而,这种应用主要是侧重于数值逼近的应用。我们发现该方法若运用得当还具有滤波的功能。

事实上,多核函数在具体应用中的效果取决于如何选择,对于不同的描述对象可能会有不同的选择。所以,在实际应用中主要取决于使用者对多核函数及描述对象的理解。一般来说主要考虑以下三个方面:① 核函数的选择; ② 核函数中的参数选择; ③ 核函数的核点选择。由于我们所描述的对象是较大尺度的地壳应变场,而地壳介质由于其各向非同性、活动构造带的存在及应力大小与方向的空间变化等,使得应变场具有广谱性。兼顾数值逼近与滤波,本研究选用的核函数为“钵”形函数,具体为

s_j(x,y,x_j,y_j)=d^a_j+b.(2)

式中,d_j=((x-x_j)²+(y-y_j)²)^1/2。选用如此简单的核函数是因为使用该函数不存在量纲因子,使用者易于掌握,效果比较好。式中a的大小在核点确定之后影响波形,a越大波形的“峰”、“谷”性越突出。一般情况下,它的取值范围以1～2.5为宜,由于观测值带有误差,当测点不均匀时不宜取较大值。b以前被称为滤波因子,经验表明这种称呼不准确,因为它的大小对波形几乎无影响,但不可以为0。经过大量的试算发现,当测站分布不均匀或运动结果误差较大时,取值较大则影响数值逼近的正确性。为了稳妥起见,本文中取a=1.1、b=1。事实上,控制滤波效果或滤波成份,在宏观上由核点的空间分布所决定; 一般来说,核点间隔越大,得到的就是较低频段的结果,反之则为较高频段的结果。由于较大尺度应变场的描述是在球面上,故将(1)、(2)式合并为下式

〖JB<4{〗f(λ,φ)=∑c_js_j(λ,φ,λ_j,φ_j)=S^TC^T,

s_j(λ,φ,λ_j,φ_j)=d^1.1_j+1.

(3)

式中,d_j为球面上两点间的大地线长度,S^T=(s₁,s₂,…,s_nx),C^T=(c₁,c₂,…,c_nx)。

在GNSS计算中通常以ITRF作为参考框架获取速度场。现假定(λ_j,φ_j)为测站j的位置坐标,v_e(λ_j,φ_j)、v_n(λ_j,φ_j)为其相应的东向和北向运动结果,在用多核函数进行数值逼近时应将所有测站位置作为核函数的核点位置,因此依据(3)式,对于任一方向均可列出n个方程(测站的个数),故可求解n个未知数(C)。这时可获得东、北向运动的解析式f_e(λ,φ)、f_n(λ,φ)。由此即可计算研究区任意位置的水平运动结果。即

f(λ,φ)=f_e(λ,φ)i+f_n(λ,φ)j.(4)

同理,利用测站东、北向运动的误差也可获得相应误差的解析式:

m(λ,φ)=m_e(λ,φ)i+m_n(λ,φ)j.(5)

在球面坐标系非直角坐标系,若将平面应变计算式或由运动与应变关系式移植到球面上则必然存在偏差。根据石耀霖和朱守彪(2006)、刘序俨等(2007)所给出的球面应变计算式,只要基准恰当即可获取球面应变的无偏差结果。因此,若从参考基准的角度来看,球面上应变的正确计算似乎与运动场计算的参考基准有关。这一结论与经典应变分析所得到的与起算基准无关的结论相悖,故此该算式不被广泛地认同或使用者不知如何正确使用该模型。事实上,这是一种假象,球面应变计算式包含的非微分项应变计算恰恰是对偏差所进行的修正。对于如何正确使用该算式,杨国华等(2010)也给出了描述,只要利用ITRF参考框架所得到的水平运动(位移)场即可。现假定东西向应变为ε_e(λ,φ)、南北向应变为ε_n(λ,φ)、它们之间的剪应变为ε_en(λ,φ),在(4)式基础上和现行球面坐标系统下球面应变算式则为

{ε_e(λ,φ)=1/(Rcosφ)(f_e(λ,φ))/(λ)-(f_n(λ,φ))/Rtanφ,

ε_n(λ,φ)=1/R(f_n(λ,φ))/(φ),

ε_en(λ,φ)=1/2[1/(Rcosφ)(f_n(λ,φ))/(λ)+

1/R(f_e(λ,φ))/(φ)+(f_e(λ,φ))/Rtanφ].(6)

式中R为地球的平均半径。

在椭球面上的应变算式则为

{ε_e(λ,φ)=1/(R_ncosφ)(f_e(λ,φ))/(λ)-(f_n(λ,φ))/(R_n)tanφ,

ε_n(λ,φ)=1/(R_m)(f_n(λ,φ))/(φ),

ε_en(λ,φ)=1/2[1/(R_ncosφ)(f_n(λ,φ))/(λ)+

1/(R_m)(f_e(λ,φ))/(φ)+(f_e(λ,φ))/(R_n)tanφ].(7)

式中,R_n为卯酉圈曲率半径,R_m为子午圈曲率半径,λ为经度,φ为纬度。

上述数值逼近所得的水平运动场的解析结果,不仅包含运动信息,也包含误差干扰,所以由此得到的应变结果也必然不纯正。因此,进行滤波与信息分离是必要的。由于多核函数也具有滤波的功能,故本文利用多核函数法进行空间滤波。在进行滤波计算时,滤波效果或滤波层次的控制则取决于核点的选用与分布。从性质上看,核点的间隔越小,获得的应变信息所包含的频谱带宽越宽,可涵盖从低频至高频的范围,反之带宽越窄,也越趋于低频。在进行滤波计算时,由于应变场是由上述3个参数共同描述,而其他应变参数是通过这3个参数计算得到,所以为了保持应变场描述的协调性和彼此之间的相关性不被破坏,进行该3参数的滤波计算时应采用相同的滤波计算准则。因此,应变参数的滤波函数均表述为

{ε_e(λ,φ)=∑a_is_i(λ,φ,λ_i,φ_i)=S^TA^T,

ε_en(λ,φ)=∑b_is_i(λ,φ,λ_i,φ_i)=S^TB^T,

ε_n(λ,φ)=∑c_is_i(λ,φ,λ_i,φ_i)=S^TC^T.

(8)

式中,A^T=(a₁,…,a_nx), B^T=(b₁,…,b_nx), C^T=(c₁,…,c_nx), 它们均为待定系数; S^T=(s₁,…,s_nx), s_i的具体形式如(3)所示,(λ_i,φ_i)为核点坐标。

在实际计算时,首先以相同的网格化形式获取各类应变数据,并以此作为应变观测值。然后,根据最小二乘法求解上述待定系数,即对每一类应变通过误差方程

V=DX+L,(9)

建立相应的法方程

NX+W=0,(10)

并由此求解(待定系数)为

X=-N^-1W.(11)

其中,N=D^TPL,W=D^TPD,L为应变观测值的矩阵,P为应变观测值的权矩阵(p_i=1/m²_i,m_i在下一节给出具体算式)。

求得X后,(8)式则变为经滤波后应变计算的函数解析式,据此可获得其他应变的函数解析式,即最大主应变为

ε_max(λ,φ)=(ε_e(λ,φ)+ε_n(λ,φ))/2+

((4ε²_en(λ,φ)+[ε_e(λ,φ)-ε_n(λ,φ)]²)^1/2)/2,(12)

最小主应变为

ε_min(λ,φ)=(ε_e(λ,φ)+ε_n(λ,φ))/2-

((4ε²_en(λ,φ)+[ε_e(λ,φ)-ε_n(λ,φ)]²)^1/2)/2,(13)

最大剪应变为

γ_max(λ,φ)=(ε_max(λ,φ)-ε_min(λ,φ))/2

=((4ε²_en(λ,φ)+[ε_e(λ,φ)-ε_n(λ,φ)]²)^1/2)/2,(14)

面应变为

Δ(λ,φ)=ε_max(λ,φ)+ε_min(λ,φ)

=ε_e(λ,φ)+ε_n(λ,φ),(15)

最大主应变方向为

θ(λ,φ)=arctan((ε_max(λ,φ)-ε_e(λ,φ))/(ε_en(λ,φ))).(16)

上述有关章节只给出了无偏应变的算式((6)～(8)式),但没有给出相应的误差估计。该误差估计不仅可使我们了解无偏应变的误差,更重要的是该误差是此后滤波计算的定权依据,以及有关应变误差传递计算所需要的相关系数。故此我们给出具体算式,为了简单起见做了适当简化,即误差传播时只考虑微分项,因为修正项的误差较微分项的误差小数十倍; 此外,由于协方差矩阵太大,通常给研究者所提供的GNSS计算结果中不包含不同测站速度间的协误差,所以以下计算式不再考虑它们之间的相关性问题。因此,依(6)式并兼顾(5)式及实际计算时的网格化特点,各个应变的误差为

{m_εe(λ,φ)=1/(Δs_λ)(m²_e(λ,φ)+m²_e(λ+Δλ,φ))^1/2,

m_εn(λ,φ)=1/(Δs_φ)(m²_n(λ,φ)+m²_n(λ,φ+Δφ))^1/2,

m_εen(λ,φ)=1/2(1/(Δs_λ)(m²_n(λ,φ)+m²_n(λ+Δλ,φ))^1/2+

1/(Δs_φ)(m²_e(λ,φ)+m²_e(λ,φ+Δφ))^1/2).(17)

其中Δs_λ、Δs_φ、Δλ、Δφ分别为应变网格化计算时经、纬方向上的距离步长和经纬步长。各应变之间的协方差为

{m_εe,n(λ,φ)=0,

m_εe,en(λ,φ)=1/(2Δs_φΔs_λ)m²_e(λ,φ),

m_εn,en(λ,φ)=1/(2Δs_φΔs_λ)m²_n(λ,φ).(18)

其相应的相关系数为

{R_e,n(λ,φ)=0,

R_e,en(λ,φ)=(m_εe,en(λ,φ))/(m_εn(λ,φ)m_εen(λ,φ)),

R_n,en(λ,φ)=(m_εn,en(λ,φ))/(m_εn(λ,φ)m_εen(λ,φ)).(19)

进行滤波计算后,由(10)式利用测站位置的残差可计算相应的应变单位权方差估值,东向、剪切和北向分别为

m²_0e=(V^T_eP_eV_e)/(n-n_x),(20)

m²_0en=(V^T_enP_enV_en)/(n-n_x),(21)

m²_0n=(V^T_nP_nV_n)/(n-n_x).(22)

式中,n_x为核点的个数,n为测站个数,且大于n_x,否则无法给出精度评定,也达不到滤波的目的。(20)～(22)式中的残差均是测站点位置的相应残差。根据误差传播定律,可求各应变参量的误差。令

Q_j=N^-1_j.(23)

式中,j=ee,enen,nn,因此东西向、剪切和南北向应变方差分别为

m²_ee(λ,φ)=m²_0eS^TQ_eeS,(24)

m²_enen(λ,φ)=m²_0enS^TQ_enenS,(25)

m²_nn(λ,φ)=m²_0nS^TQ_nnS,(26)

东西向应变与剪应变协方差为

m_e,en(λ,φ)=m_ee(λ,φ)·m_enen(λ,φ)·R_e,en(λ,φ),(27)

南北向应变与剪应变协方差为

m_n,en(λ,φ)=m_nn(λ,φ)·m_enen(λ,φ)·R_n,en(λ,φ),(28)

最大主应变方差为

m²_εmax(λ,φ)=f ²₁(λ,φ)m²_ee(λ,φ)+f ²₂(λ,φ)m²_enen(λ,φ)+f ²₃(λ,φ)m²_nn(λ,φ)+2f₁(λ,φ)f₂(λ,φ)_me,en(λ,φ)+2f₂(λ,φ)f₃(λ,φ)m_en,n(λ,φ),(29)

最小主应变方差为

m²_rmin(λ,φ)=f ²₃(λ,φ)m²_ee(λ,φ)+f ²₂(λ,φ)m²_enen(λ,φ)+f ²₁(λ,φ)m²_nn(λ,φ)-2f₃(λ,φ)f₂(λ,φ)_me,en(λ,φ)-2f₂(λ,φ)f₁(λ,φ)m_en,n(λ,φ),(30)

最大剪切应变方差为

m²_rmax(λ,φ)=[h²₁(λ,φ)m²_ee(λ,φ)+h²₂(λ,φ)m²_enen(λ,φ)+h²₃(λ,φ)m²_nn(λ,φ)+2h₁(λ,φ)h₂(λ,φ)_me,en(λ,φ)+2h₂(λ,φ)h₃(λ,φ)m_en,n(λ,φ)]/2,(31)

面应变方差为

m²_Δ(λ,φ)=m²_ee(λ,φ)+m²_nn(λ,φ).(32)

式中,

f₁(λ,φ)=(2γ_max(λ,φ)+ε_e(λ,φ)-ε_n(λ,φ))/(4γ_max(λ,φ)),(33)

f₂(λ,φ)=(ε_en(λ,φ))/(γ_max(λ,φ)),(34)

f₃(λ,φ)=(2γ_max(λ,φ)+ε_n(λ,φ)-ε_e(λ,φ))/(4γ_max(λ,φ)),(35)

h₁(λ,φ)=(ε_e(λ,φ)-ε_n(λ,φ))/(4γ_max(λ,φ)),(36)

h₂(λ,φ)=(ε_en(λ,φ))/(γ_max(λ,φ)),(37)

h₃(λ,φ)=(ε_n(λ,φ)-ε_e(λ,φ))/(4γ_max(λ,φ)).(38)