大地测量学科发展的总趋势是向地球科学纵深发展，深入到其他地学学科的交叉领域，其主要任务是监测和研究地球动力学现象，研究地球本体的各种物理场，认识与探索地球内部的各种物理过程并揭示其规律。以空间大地测量为标志的现代大地测量技术，不论在测量的空间尺度上还是已达到的精度水平方面，已经有能力监测地球动力学过程产生的运动状态和物理场的微变化^[1]。利用GPS监测点复测资料所得的位移或速度研究区域地壳运动与应变是大地测量深入地球科学研究最常用的方法之一^[2-7]，在数据处理过程中基本都采用最小二乘方法(least squares, LS)，没有考虑到反演模型中系数矩阵误差的影响。文献[8]采用总体最小二乘方法研究地壳应变参数的反演，考虑了系数矩阵的误差影响以及系数矩阵中含有非随机元素的情况，但是在公式推导时对公式进行了线性化处理，引入了模型误差，因此得到的参数结果存在偏差。目前关于地壳应变参数反演模型中系数矩阵误差对模型参数的具体影响情况，还没有明确的定论。

部分变量误差(partial-EIV，PEIV)模型作为变量误差(errors-in-variables，EIV)模型一般化的表达形式，既考虑了函数模型中观测向量的误差，又顾及了系数矩阵含有误差的情况，同时对系数矩阵的结构及组成无特殊限制，更适合用于EIV模型问题的求解。目前，国内外研究并使用PEIV模型的文献并不多。文献[9]在PEIV模型的基础上，推导出了加权总体最小二乘(weighted total least squares, WTLS)更一般化的解，该算法涵盖了一般WTLS算法需要特殊处理的各种情况，同时大大减少了模型的待估量，算法的简单形式使后续估计值的精度评定更加便利；文献[10]提出了基于PEIV模型的可靠性理论，并分析了模型中观测向量和系数矩阵可靠性的基本特点。本文以PEIV模型为基础，采用加权总体最小二乘算法进行地壳应变参数反演，进一步揭示GPS站点坐标误差对地壳应变模型参数估计的影响。同时本文算法可以很好地解决反演模型的系数矩阵中存在的随机元素、非随机元素和重复元素的情况，算法的推导公式中也避免了文献[8]中存在的问题，结合推导得到的最小二乘改正项公式，对地壳反演模型中坐标点误差对反演参数求解的影响进行了分析。通过模拟数据和川滇地区的实际数据，分析了系数矩阵误差对地壳应变参数反演的影响，展示了更多反演过程的细节，对深入理解地壳应变反演机理有一定的价值。

1 由坐标位移反演地壳应变参数的PEIV模型算法

在假定介质均匀且邻近测点间应变均匀的情况下，可以建立邻近点间相对形变量与地壳应变量的线性关系。设第j点两期纵坐标位移u_j，横坐标位移v_j，网中共有m点，由坐标位移反演应变参数的公式(二维)为^[11]：

式中，x_j、y_j(j=1, 2, …, m)为第j点近似坐标；u、v为网的平移量；w为其转动量；ε_x、ε_y、ε_xy为应变状态参数。

将式(1) 改写成矩阵形式：

令y=[u₁ v₁ … u_m v_m]^T，y∈R^2m×1；β= [u v ε_x ε_y ε_xy w]^T, β∈R^6×1; A∈R^2m×6

所以，式(2) 可以表示为：

基于PEIV模型的坐标位移反演应变参数的思想是：不仅考虑坐标位移的误差，即y的误差对应变参数的影响，以及系数矩阵A误差的影响，同时考虑系数矩阵的元素并非全都是随机的，而部分是由0或者1等常数组成的，所以有^[9]:

式中，y∈R^2m×1表示观测向量；ε∈R^2m×1表示观测向量y的随机误差；h∈R^12m×1主要由系数矩阵A中非随机元素组成；B∈R^12m×t是固定矩阵，其阶数取决于系数矩阵A中随机元素的数目(此处假设为t)；I_2m ∈R^2m×2m为单位阵；A∈R^t×1是系数矩阵A中随机元素组成的列向量(若随机元素重复出现，则只提取一次)，其真值用a表示；ε_a∈R^t×1是A中包含的随机误差；A·β=(β^T⊗I_n)vec(A)，其中A表示系数矩阵A的真值，vec(·)表示矩阵拉直运算。

地壳应变参数反演的PEIV模型中各矩阵的具体形式表示为：

其中⊗表示直积运算。

平差的随机模型为：

式中，D表示方差-协方差矩阵；σ²为单位权方差；Q_εε=P₁^-1为观测值的协因数阵；Q_εaεa=P₂^-1为系数矩阵中随机元素的协因数阵。

通过观测值y和系数矩阵随机元素A估计PEIV模型中a和β的加权总体最小二乘平差准则为^[9]：

当初始权阵P₁和P₂已知时，可以采用下面的步骤进行未知参数求解^[9]

第1步  初始值设定。给定观测向量y和系数矩阵A，提取矩阵A中随机元素组成列向量a，根据系数矩阵结构，构造常数向量h和固定矩阵B，计算参数初始值$ {{\mathit{\boldsymbol{\hat \beta }}}^{\left( 0 \right)}} = {\left( {{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}A} \right)^{ - 1}}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{y}} $；

第2步  迭代过程，利用公式$ \widehat {\mathit{\boldsymbol{\overline a}} }\mathit{\boldsymbol{ = }}{\left( {\mathit{\boldsymbol{S}}_\beta ^{\rm{T}}{\mathit{\boldsymbol{P}}_1}{\mathit{\boldsymbol{S}}_\beta } + {\mathit{\boldsymbol{P}}_2}} \right)^{ - 1}}\left( {{\mathit{\boldsymbol{P}}_2}\mathit{\boldsymbol{a}} - \mathit{\boldsymbol{S}}_\beta ^{\rm{T}}{\mathit{\boldsymbol{P}}_1}\left( {\sum\limits_{i = 1}^m {{\mathit{\boldsymbol{h}}_\mathit{\boldsymbol{i}}}{{\mathit{\boldsymbol{\hat \beta }}}_\mathit{\boldsymbol{i}}}} } \right) + \mathit{\boldsymbol{S}}_\beta ^{\rm{T}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{y}}} \right) $计算$ \widehat {\mathit{\boldsymbol{\overline a}} } $，式中$ {\mathit{\boldsymbol{S}}_\beta } = \sum\limits_{i = 1}^m {{\mathit{\boldsymbol{B}}_\mathit{i}}{{\mathit{\hat \beta }}_\mathit{i}}} $；利用公式$ \mathit{\boldsymbol{\hat \beta }} = {\left( {{{\mathit{\boldsymbol{\hat A}}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{\hat A}}} \right)^{ - 1}}{{\mathit{\boldsymbol{\hat A}}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{y}} $计算$ {{\mathit{\boldsymbol{\hat \beta }}}^{\left( {i + 1} \right)}} $，式中$ \mathit{\boldsymbol{\hat A}}\mathit{ = }{\rm{ve}}{{\rm{c}}^{ - 1}}\left( {\mathit{\boldsymbol{h}} + \mathit{\boldsymbol{B}}\widehat {\mathit{\boldsymbol{\overline a}} }} \right) $，其中vec^-1(·)表示拉直运算的逆运算。

第3步  迭代终止。$ \left\| {{{\mathit{\boldsymbol{\hat \beta }}}^{\left( {i + 1} \right)}} - {{\mathit{\boldsymbol{\hat \beta }}}^{\left( i \right)}}} \right\| \le {\delta _0} $，其中δ₀表示设定的阈值。

2 系数矩阵误差对WTLS估计值的影响

对基于PEIV模型的地壳应变参数加权总体最小二乘解的表达式进行相应调整可以得到^[10]：

式中，$ {{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{WTLS}}}} $表示加权总体最小二乘估值；$ {{{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{LS}}}}} $表示加权最小二乘估值；$ {{{\mathit{\boldsymbol{\hat E}}}_A}} $表示通过WTLS估计得到的系数矩阵误差的残差；$ {{{\mathit{\boldsymbol{\hat e}}}_y}} $表示通过WLS估计得到的观测向量的残差。

将式(7) 展开可以得到：

令$ {\mathit{\boldsymbol{N}}_1} = {\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{A}}, {\mathit{\boldsymbol{N}}_2} = {\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}{{\mathit{\boldsymbol{\hat E}}}_A} + {{\mathit{\boldsymbol{\hat E}}}^{\rm{T}}}_A{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{A}} + {{\mathit{\boldsymbol{\hat E}}}^{\rm{T}}}_A{\mathit{\boldsymbol{P}}_1}{{\mathit{\hat E}}_A} $，根据分块矩阵求逆公式，式(8) 可以改写为：

对上式进行移项调整可以得到：

令$ \begin{array}{l} {\mathit{\boldsymbol{H}}_1} = \mathit{\boldsymbol{N}}_1^{ - 1}{{\mathit{\boldsymbol{\hat E}}}^{\rm{T}}}_A{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{A}} - \mathit{\boldsymbol{N}}_1^{ - 1}{\mathit{\boldsymbol{N}}_2}{\left( {{\mathit{\boldsymbol{I}}_6} + \mathit{\boldsymbol{N}}_1^{ - 1}{\mathit{\boldsymbol{N}}_2}} \right)^{ - 1}} - \mathit{\boldsymbol{N}}_1^{ - 1}{\mathit{\boldsymbol{N}}_2}{\left( {{\mathit{\boldsymbol{I}}_6} + \mathit{\boldsymbol{N}}_1^{ - 1}{\mathit{\boldsymbol{N}}_2}} \right)^{ - 1}}\\ \;\;\;\;\;\;\;\;\mathit{\boldsymbol{N}}_1^{ - 1}{{\mathit{\boldsymbol{\hat E}}}^{\rm{T}}}_A{\mathit{\boldsymbol{P}}_1}\mathit{\boldsymbol{A}}, \;\;{\mathit{\boldsymbol{H}}_2} = \mathit{\boldsymbol{N}}_1^{ - 1}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}{{\mathit{\boldsymbol{\hat e}}}_y} + \mathit{\boldsymbol{N}}_1^{ - 1}{{\mathit{\boldsymbol{\hat E}}}^{\rm{T}}}_A{\mathit{\boldsymbol{P}}_1}{{\mathit{\boldsymbol{\hat e}}}_y} - \\ \;\;\;\;\;\;\;\;\mathit{\boldsymbol{N}}_1^{ - 1}{\mathit{\boldsymbol{N}}_2}{\left( {{\mathit{\boldsymbol{I}}_6} + \mathit{\boldsymbol{N}}_1^{ - 1}{\mathit{\boldsymbol{N}}_2}} \right)^{ - 1}}\mathit{\boldsymbol{N}}_1^{ - 1}{\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{P}}_1}{{\mathit{\boldsymbol{\hat e}}}_y} - \\ \;\;\;\;\;\;\;\;\mathit{\boldsymbol{N}}_1^{ - 1}{\mathit{\boldsymbol{N}}_2}{\left( {{\mathit{\boldsymbol{I}}_6} + \mathit{\boldsymbol{N}}_1^{ - 1}{\mathit{\boldsymbol{N}}_2}} \right)^{ - 1}}\mathit{\boldsymbol{N}}_1^{ - 1}{{\mathit{\boldsymbol{\hat E}}}^{\rm{T}}}_A{\mathit{\boldsymbol{P}}_1}{{\mathit{\boldsymbol{\hat e}}}_y} \end{array} $。

式(9) 可以改写为：

式中，b表示WTLS估值和WLS估值之间的改正项，由于采用最小二乘求解地壳应变参数反演模型得到的是一个近似解，在求解过程中未考虑函数模型中系数矩阵误差的影响，且系数矩阵误差对WLS估计值的影响未知，因此解算结果的优劣无法判断。采用总体最小二乘进行地壳应变参数的求解，得到的是一个严密解，此时在WTLS估值和WLS估值之间需要添加一个改正项。在式(10) 的基础上进行模拟算例和实际算例的定量分析，可以对最小二乘准则和总体最小二乘准则之间的差异性有一个初步的了解，同时可以客观研究系数矩阵误差对地壳应变参数WTLS解的影响。

3 算例分析

3.1   模拟实验

对3个方案(见表 1)进行模拟实验，以方案1为例，在边长为100 m²的正方形区域，沿x和y轴每隔10 m进行取样，共得到121个坐标点，如图 1。给定该区域应变参数ε_x=200×10^-8，ε_y=-80×10^-8，ε_xy=1 000×10^-8，ω=1 100×10^-8(模拟算例中未考虑平移量)，将坐标数据代入式(2) 得到各点对应的位移量真值。分别给坐标点和位移量加入服从正态分布的随机误差，此时函数模型式(4) 中，观测向量y由加入误差后的位移量组成，其中D(ε)=σ²Q，σ²=1 cm²，Q=diag(0.1, 0.1, …, 0.1)；而系数矩阵中的随机元素A由含有误差的坐标组成，且满足D(ε_a)=σ²Q_εaεa，Q_εaεa=diag(0.3, 0.3, …，0.3)。采用加权最小二乘(WLS)、基于PEIV模型的加权总体最小二乘法(PEIV-WTLS)两种方法求解表 1中3种方案的应变参数，同时计算相应的改正项b用于衡量WLS与PEIV-WTLS法解算结果之间的差距。表 2为方案1、方案2、方案3的应变参数计算结果。图 2给出了利用基于PEIV模型的加权总体最小二乘法(PEIV-WTLS)求解得到的各方案观测向量和系数矩阵随机元素残差值。在算例模拟过程中，为了使模拟的位移量符合实际情况，在方案2和方案3中，当模拟区域扩大时，相应的应变参数需要进行调整。

表  1  各方案的情景设置

Table  1.  The Scene of the Three Schemes

	方案1	方案2	方案3
区域	100 m×100 m	200 m×200 m	1 000 m×1 000 m
采样间隔	沿x、y轴每隔10 m	沿x、y轴每隔20 m	沿x、y轴每隔100 m
坐标点总数(采样)	121个	121个	121个
模拟应变参数真值	εx=200×10-8 εy=-80×10-8 εxy=1 000×10-8 ω=1 100×10-8	εx=200×10-9 εy=-80×10-9 εxy=1 000×10-9 ω=1 100×10-9	εx=200×10-10 εy=-80×10-10 εxy=1 000×10-10 ω=1 100×10-10

下载: 导出CSV 

| 显示表格

图  1  方案1中地壳应变参数反演的模拟站点

Figure  1.  Simulated Network Stations of the Inversion of Crustal Strain in Case 1

下载: 全尺寸图片 幻灯片

表  2  各方案参数计算结果

Table  2.  Results of Parameters in Each Case

		εx	εy	εxy	w
方案1	真值	200×10-8	-80×10-8	1 000×10-8	1 100×10-8
	WLS	189×10-8	-86×10-8	1 006×10-8	1 097×10-8
	PEIV-WTLS	189×10-8	-86×10-8	1 006×10-8	1 097×10-8
	b	-1.5×10-18	2.6×10-17	-1.5×10-17	-1.7×10-17
方案2	真值	200×10-9	-80×10-9	1000×10-9	1100×10-9
	WLS	195×10-9	-81×10-9	1005×10-9	1100×10-9
	PEIV-WTLS	195×10-9	-81×10-9	1005×10-9	1100×10-9
	b	-3.2×10-21	8.9×10-21	-4.8×10-21	-7.1×10-21
方案3	真值	200×10-10	-80×10-10	1 000×10-10	1 100×10-10
	WLS	202×10-10	-83×10-10	999×10-10	1 103×10-10
	PEIV-WTLS	202×10-10	-83×10-10	999×10-10	1 103×10-10
	b	-2.9×10-23	9.8×10-23	-1.0×10-22	-1.2×10-22

下载: 导出CSV 

| 显示表格

图  2  各方案观测向量和系数矩阵随机元素残差值(PEIV-WTLS)

Figure  2.  Residuals of Observation Vector and Coefficient Matrix in Each Case(PEIV-WTLS)

下载: 全尺寸图片 幻灯片

从表 2中可以看出，采用WLS和PEIV-WTLS这两种算法求解得到的应变参数估值基本相同，无法看出两者之间的区别，因此通过分析表格中的改正项b对这种现象进行探索。从式(10) 中可以看出，影响改正项b大小的主要因素来自地壳应变参数的加权最小二乘估值$ {{\mathit{\boldsymbol{\hat \beta }}}_{{\rm{WLS}}}} $，以及矩阵H₁和H₂中的元素(主要为系数矩阵的残差$ {{\mathit{\boldsymbol{\hat E}}}_A} $和WLS估计得到的观测向量的残差$ {{\mathit{\boldsymbol{\hat e}}}_y} $)。在地壳应变参数反演模型中，系数矩阵的随机元素是由观测站点的坐标组成的，站点坐标所含误差的量级一般为10^-3~10^-2，观测向量元素由位移量组成，残差$ {{\mathit{\boldsymbol{\hat e}}}_y} $的量级一般在10^-3~10^-2，应变参数估值的量级一般与观测站点坐标的量级呈负相关。从表 1中可以看出，随着模拟区域扩大，站点坐标值增大，应变参数随之减小。比较方案1和方案2可以看出，当模拟区域扩大4倍，即坐标点观测值量级扩展到10²时，相应的改正项减小10^-4个量级；比较方案2和方案3，当模拟区域扩大25倍，即坐标点观测值量级扩展到10³时，相应的改正项减小10^-2个量级，产生这一现象的主要原因为随着坐标点观测值量级的增大，求解得到的应变参数变小，从式(3) 中可以看出，系数矩阵中随机元素的误差是通过与应变参数的乘积形式对函数模型产生影响的，因此当应变参数减小的同时也会削弱系数矩阵中随机元素对模型解算的影响。也可以解释为，由于反演模型系数矩阵基本是由监测点坐标组成的，数值一般在10⁶量级，坐标观测值中含有的误差被坐标值本身淹没，导致计算结果与WLS结果基本相同。在许多研究坐标转换的文献中也存在相同的问题^[12]。从图 2可以看出，采用严密的PEIV-WTLS算法可以求解得到观测向量的残差值，却无法恢复系数矩阵中随机元素的残差值。区域性形变研究是一个长期的过程，任何微小的变化经过长时间的积累都会达到一个可观的量级，因此在监测地壳形变和运动时，再微小的变化也需要考虑在内，地壳应变参数反演模型中系数矩阵的误差随着观测区域的扩大会逐渐减小其对整个模型参数估计的贡献，却不能忽略其微小的影响，需要采用严密的方法将其考虑在内。

3.2   川滇地区地壳应变参数反演

在实际情况中，GPS监测站点的分布一般比较分散，且数量较少，站点质量存在差异。这也给区域性地壳形变研究带来了更多的不确定性。本文采用GPS速度场数据对川滇地区8个块体的应变率参数进行反演，以此来探索在实际情况中是否也存在上述3个模拟方案中相同的情况。川滇地区主要包括北纬20°~34°、东经96°~106°的中国大陆西南区域的青藏高原东南部，分布于其中的大量活动断层将该区域分成了许多块体，主要的断层有红河断层、小江断层、腾冲-景洪断层、鲜水河断层、安宁河断层、丽江断层、金沙江断层和龙门山断层等，这些断层所分割成的块体有马尔康块体(Ⅰ)、川西北块体(Ⅱ1)、滇中块体(Ⅱ2)、保山块体(Ⅲ1)、景谷块体(Ⅲ2)、勐腊块体(Ⅲ3)、羌塘块体(Ⅳ)和华南块体(Ⅴ)^[13]，各块体与断层之间的空间关系见图 3。

图  3  川滇地区GPS监测点分布图^[14]

Figure  3.  The Distribution of GPS Stations in the Region of Sichuan and Yunnan^[14]

下载: 全尺寸图片 幻灯片

本文所用的GPS水平速度场的数据来自于文献[15]，即1998~2004年中国地壳运动观测网络的GPS水平速度场及其方差-协方差矩阵数据，主要来源于中国地壳运动观测网络。在应变率参数反演过程中，首先将每个GPS速度场数据对应的测站的经、纬度坐标经过高斯投影转换为高斯平面直角坐标(区域坐标)，在坐标转换过程中，选取每个块体中心单独投影^[16]。由于GPS测站的经纬度坐标的测定必然含有误差，而且在坐标转换中也会引入其他误差，而本文采用的地壳应变参数反演的PEIV模型求解方法可以顾及这些坐标误差。利用式(10) 计算各个块体的改正项b，具体计算结果见表 3。此处需要明确，当采用的观测数据为位移时，计算得到的是地壳应变参数; 当采用的观测数据为速度时，计算得到的是应变率参数。图 4给出了利用基于PEIV模型的加权总体最小二乘法求解得到的各个块体反演过程中观测向量和系数矩阵随机元素残差值，图中实线表示观测向量残差，虚线表示系数矩阵残差。

表  3  各块体应变率参数计算改正项b结果

Table  3.  The Calculation Results of Each Block Strain Rate Parameter

改正项	u/a	v/a	εx/a	εy/a	εxy/a	ω/a
bⅠ	-5.6×10-16	-2.2×10-15	-4.4×10-22	8.7×10-24	1.9×10-23	6.9×10-23
bⅡ1	2.2×10-15	-1.2×10-15	3.1×10-23	-1.6×10-23	-3.4×10-23	-5.0×10-23
bⅡ2	2.2×10-15	7.5×10-18	1.2×10-24	-4.5×10-25	-2.7×10-24	-4.3×10-24
bⅢ1	-2.7×10-15	9.4×10-15	1.6×10-20	-1.1×10-21	-2.7×10-21	-2.5×10-21
bⅢ2	8.6×10-16	2.1×10-15	1.7×10-23	-2.4×10-23	-2.8×10-24	-1.2×10-24
bⅢ3	-3.4×10-13	5.4×10-13	5.3×10-22	2.2×10-22	3.0×10-22	1.8×10-22
bⅣ	1.8×10-13	1.8×10-13	-1.3×10-20	-4.8×10-22	-6.9×10-21	-6.4×10-21
bⅤ	-5.6×10-16	-3.7×10-16	-1.9×10-24	7.1×10-25	5.4×10-26	1.2×10-24

下载: 导出CSV 

| 显示表格

图  4  川滇地区地壳应变参数反演过程中观测向量和系数矩阵随机元素残差值

Figure  4.  Residuals of Observation Vector and Coefficient Matrix in the Region of Sichuan and Yunnan

下载: 全尺寸图片 幻灯片

1) 在利用GPS监测点复测资料所得的位移、速度反演应变或应变率参数时，用最小二乘方法未考虑由监测点坐标组成的反演模型系数矩阵的误差影响，而PEIV-WTLS算法考虑了观测向量的误差，同时考虑了反演模型系数矩阵的扰动对反演应变参数的影响，从理论上，PEIV模型算法更加合理。

2) 结合图 3和表 3可以看出，保山块体(Ⅲ1) 区域内的GPS监测点分布较为集中，且覆盖区域较小，因此相比于其他块体，改正项b的量级更大些；滇中块体(Ⅱ2) 和羌塘块体(Ⅳ)区域内的GPS监测点覆盖范围较广，因此相比于其他块体，改正项b的量级更小些，这也与上述模拟方案中得到的结论相一致。

3)表 3中的改正项b主要用于衡量WLS与PEIV-WTLS算法解算结果之间的差距。通过比较不同块体地壳应变率参数反演改正项b的大小，可以更直观地分析影响地壳应变率参数反演过程的因素。地壳应变率参数反演模型中系数矩阵的随机元素由监测点坐标组成，含有随机误差，若不考虑其影响，则会导致解算结果存在偏差。

4) 从图 4中可以看出，利用PEIV-WTLS算法可以估算出各块体相应区域中GPS监测站点位移中包含的随机误差，而无法估算出GPS监测站点坐标含有的随机误差，即系数矩阵中的随机误差，这与模拟方案中的结果相同，主要是由于GPS监测站点坐标值量级过大，相应的应变参数偏小，较小的应变参数抑制了系数矩阵中的随机误差对模型估计的贡献度。

5) 区域性地壳应力应变研究需要极高的测量精度，因此在进行地壳应变率参数求解时必须采用理论严密的算法。虽然系数矩阵中随机元素的误差相对于随机元素本身量级过小，但在一定程度上也会影响模型估计过程，因此必须将其考虑在内。

4 结语

本文将以PEIV模型为基础的加权总体最小二乘算法应用于地壳应变参数反演，该算法考虑了以往一直被忽略的系数矩阵误差的影响，得到了合理的参数解。通过改正项b定量分析系数矩阵误差对参数估计的影响。通过对模拟数据和川滇地区的实际数据进行处理，可以得知系数矩阵误差对地壳应变参数反演的影响主要受到GPS站点坐标值量级以及应变参数量级的牵制。同时存在加权总体最小二乘算法无法恢复系数矩阵随机元素误差这一问题，对于这一问题的解决方法还有待进一步研究。