InSAR已被证实为一种强大的空间对地观测技术, 可以用来研究许多种地球物理现象^{[1, 2, 3, 4]}。然而, 如同其他微波遥感一样, 当雷达信号2次穿越大气层时, 会因大气水汽含量在时间和空间上的不确定性而产生不同程度的传播延迟, 从而影响了差分InSAR的监测精度^{[5, 6, 7, 8]}。因此, 大气延迟成为影响差分InSAR高精度获取缓慢微小地表形变的重要因素。作为一种新型的大气探测手段, 20世纪90年代发展起来的地基GPS可以实时、连续、全天候、高精度地获取大气可降水汽值, 且精度可以达到1~2 mm ^{[9, 10]}。随着越来越多的区域连续运行GPS网的建立, 为利用GPS水汽观测值建立区域水汽延迟改正模型, 从而对SAR影像中的大气延迟进行改正提供了可能性。但是, GPS在区域上是以离散点形式分布的, 即使相距最近的区域性GPS连续跟踪站间距一般也有几km到几十km不等, 远远大于SAR影像干涉像元的实际大小。因此, 在利用GPS观测数据改正InSAR影像中的大气延迟时, 必须要对GPS水汽观测值进行空间插值, 以建立与干涉图像元的空间分辨率相当的水汽延迟改正模型。以往的研究主要是用逆距离权(inverse distance weighted, IDW)倒数和克里格(Kriging)等插值方法直接对GPS水汽值进行空间插值^[11], 但这些方法并未考虑到地形对大气延迟的影响, 效果并不理想。

本文提出了一种改进的克里格插值(IKriging)方法。该方法在分别考虑大气在高程和距离上变异特性的基础上, 按照双因素对插值结果的影响, 建立一个综合的样本变异模型, 并根据该变异模型建立已知数据与待插值点之间的协方差矩阵, 求解已知数据对待插点的贡献权值, 进而确定待插值点值。

Kriging算法是法国地质统计学家Matheron在1963年提出的^[12], 用于研究空间分布既有随机性又有结构性的空间变量, 并可以对这些变量进行最优无偏估计。该法是以研究变量的空间相关性和依赖性, 或空间结构和变异性为基础^[13]。在二阶平稳假设条件下, 变量的空间变异函数γ (h)可以表示为

γ (h)= 12E[Z(x)-Z(x+h)]², (1)

式中: 变异函数γ (h)是一个仅与2个采样点之间距离h有关、而与采样点位置x无关的函数; Z(x)表示空间任意采样点; E表示求取期望值。对于任意待估计点的估计值Z'(x₀)均可以通过待估测点搜索范围内n个观测样本Z(x_j)(j=1, 2, …, n)的线性组合, 即

Z'(x₀)= ∑j=1nλ _jZ(x_j)(2)

求解。式中λ _j为权重系数。

由此可见, 内插估计值的精度取决于权重的求解。根据最优无偏估计原则, λ _j的确定应满足^[13]

E[Z'(x₀)-Z(x₀)]=0, (3)

E[Z'(x₀)-Z(X₀)]²=min, (4)

式中Z(x₀)为其真值。

根据拉格朗日定理, λ _j可由式

∑i, j=1nλjC¯ij+μ=C¯i0∑j=1nλj=1(5)

求解^[14]。式中: μ 为拉格朗日参数; C¯ij表示i, j 2个样本点间的变异函数值(i, j=1, 2, …, n ); C¯i0表示第i个样本点与待估点间的变异函数值。由式(5)求出各样本点的权重值之后, 再根据式(2)即可求得待估点值。

从以上分析可知, 应用Kriging方法进行插值时, 考虑的仅是一维的距离因素。但在某些情况下, 影响待估点值的因素不仅与距离有关, 还与其他因素有关, 如本文研究的大气延迟影响, 其各待估点的值不仅与距离有关, 还与高程有关, 如果仅用一维Kriging插值, 显然不能得到较为合理的估计值。为此, 本文提出了改进的Kriging(improved Kriging, IKriging)插值方法, 它可以考虑多维因素对插值结果的影响。以考虑2个影响因素的情况为例, IKriging算法如下:

已知样本数据为n个GPS站点三维坐标(x_i, y_i, z_i)以及可降水汽值p_i(i=1, 2, …, n)。令待求位置(x₀, y₀, z₀)上的可降水汽值为p₀。如果i和j分别代表任意2个已知GPS站点, 那么各站点之间的距离l_ij和高差绝对值z_ij可通过

l_ij= (xi-xj)2+(yi-yj)2(6)

z_ij= zi-zj(7)

求解。这样, 即可根据各样本站点上的可降水汽值p_i和i, j 2个站点间的距离l_ij, 由式(1)分析出p_i在距离上的变异函数γ (l); 同样, 由各样本站点上的可降水汽值p_i和高差绝对值z_ij分析出p_i在高程上的变异函数γ (z)。在假定距离和高程变化对水汽插值具有等权影响的条件下, GPS站点上水汽变异函数可写为

γ (l, z)= 12[γ (l)+γ (z)]。 (8)

将式(5)中的 C¯ij和 C¯i0分别用式(8)计算的样本间变异函数值矩阵表示, 则式(5)可以写作

∑i, j=1nλjγij(l, z)+μ=γi0(l, z)∑j=1nλj=1, (9)

用矩阵的形式可表示为

γ11(l, z)…γ1n(l, z)1………γn1(l, z)…γnn(l, z)11…10· λ1…λnμ= γ10(l, z)…γn0(l, z)1。(10)

由于式(10)等号左边第一项和等号右边项均为已知, 因此很容易求解出各已知数据被分配的权值λ _j。根据式(2), 将各权值λ _i与相应的水汽值p_i相乘求和, 即可得到待插值点(x₀, y₀, z₀)上的可将水汽值p₀, 即

p₀= ∑j=1nλ _jp_i。 (11)

因此, 当考虑m个影响插值结果的因素、并按照其对结果的影响程度设置权重因子α _k(k=1, 2, …, m)时, IKriging算法的变异函数就变为

γ(h)=∑k=1mαkγk(h)γk(h)=12E[Z(xk)-Z(xk+h)]2∑k=1mαk=1, (12)

式中: γ _k(h)为各因素各自的变异函数; γ (h)表示样本最终变异函数。

考虑到式(12), IKriging算法求解权重系数λ _j的式(9)就变为

(13)

同样, 由式(13)求出各样本点的权值λ _j之后, 再按照式(2)即可得到待估点值。

本文在研究Kriging方法原理的基础上, 提出了改进的IKriging方法, 并将其应用于GPS水汽的空间插值, 通过研究得到如下结论:

1)与Kriging及IDW插值方法相比, IKriging插值方法具有明显的优势, 不仅考虑了插值点与样本点之间的距离, 还考虑了插值点周围的样本点之间的空间关系, 同时还能顾及多因素对插值结果的影响。

2)通过对香港地区GPS水汽空间插值的实验表明, IKriging插值比Kriging及IDW插值更接近真值; IKriging插值方法误差小且分布集中。

3)本研究中将距离和高程变化对GPS水汽插值的影响作等权考虑, 对数值精度有一定影响, 今后将通过研究合理确定2个因素的影响权重, 进一步提高模型的插值精度。