地基GNSS技术已广泛应用于对流层水汽探测中，目前能获得的测站上空大气可降水量参数的值精度基本在2~3 mm^[1-3]。随着GNSS大气可降水量反演技术的日益成熟，很多国家和地区陆续建立了区域气象监测网，这些分布密集的测站提供的高时空分辨率的大气可降水量信息加深了我们对水汽时空变化的了解。

在此基础上，接收机至卫星方向上的水汽延迟信息包含了大气中水汽分布的各向异性部分，能够提供比大气可降水量更为丰富的信息。Flores等最先将利用局域网获得的斜路径湿延迟参数作为观测值，采用层析技术反演水汽的三维分布，为了克服观测方程的不适定性，采用附加水平约束方程及先验信息的方法实现了水汽参数的解算^[4]，目前该方法已经应用到许多国家和地区的三维水汽探测中^[5-9]。随着地基GNSS水汽探测技术的日益完善，目前我国的许多省、市均实现了准实时的大气可降水量解算；与此同时，各地建设完成及在建的测站的数目在逐渐增多，拟利用密集的GNSS区域网逐步实现三维水汽的准实时解算，实现为气象业务服务的目的。另外，随着欧洲的Galileo、俄罗斯的GLONASS及我国的Compass的建设日趋完善，空中可视卫星的数目将显著增多。从GNSS水汽层析的角度来说(观测方程为)，测站数目和可视卫星数目的增多意味着观测值数目m增大；同时，测站数目的增多意味着测站间距减小，而测站间距的减小将使水汽层析的空间分辨率得到提高，这直接导致需要解算的未知参数n的数目增大。观测方程矩阵A维数的日益增大，对水汽层析技术的数据处理方法提出了更高的要求。

代数重构算法(algebraic reconstruction technique, ART)采用迭代的方法来重构图像，它对观测值逐个进行处理，不需要直接对矩阵进行运算，其优点主要体现在解算效率和稳定性上，目前该方法已经成功地应用到电离层电子密度的层析反演中^[10]。文献[11]通过调整松弛参数向量，提高了该方法在电离层层析反演中的收敛速度；在此基础上，文献[12]对该方法在GNSS水汽层析中的应用进行了深入的研究，比较了几种不同的重构方法的优劣；王维等对联合迭代重构算法在水汽层析中的应用进行了分析^[13]。本文基于鄂东局域GNSS气象监测网，对传统的加法代数重构算法和乘法代数重构算法在水汽层析解算中的可行性进行了研究，并分析了观测值误差对解算结果的影响。

1 GNSS水汽层析

在获得了测站至卫星方向上的湿延迟参数(slant wet delay, SWD)的基础上，将对流层层析区域划分为许多小网格，并假设在一定时间段内该网格内的湿折射率参数是不随时间变化而变化的。则对某一射线j来说，该路径上的湿延迟参数SWD可表示为^[4]：

式中，n为划分的总网格数；x_i表示网格i的湿折射率参数；a_j，i表示射线j穿过网格i的长度。

将GNSS局域网内每个接收机到卫星方向上的斜路径湿延迟参数均表示为式(1)所示的求和形式(不考虑那些从层析区域侧面穿出的网格)，则水汽层析的观测方程可表示为：

式中，A为系数阵，其第ｊ行、第i列的元素为a_j，i；X为未知参数向量，即所有网格的大气湿折射率参数N_w；B为一定时间段内所有的SWD参数组成的观测值向量；m、n分别为观测值和未知参数的个数。

由于水汽层析的观测方程(2)的不适定性，通常采用附加水平和垂直约束条件的方法^[4-8]，采用最小二乘解算得到合理的层析解：

式中，H₁X=0表示水平约束方程；H₂X=C表示利用参数间的相关性或先验信息构建的垂直约束方程。

2 代数重构算法

代数重构算法是一种常用的图像重构算法，它通过迭代的方式对初始估计值重复进行修正，直到满足某种特定的条件为止。这种重构方法每次仅对一个观测值进行处理，每次迭代只涉及到简单的运算，因此十分节省计算机内存，同时具有较高的稳定度，特别适合涉及大维度矩阵的运算。下面简单介绍几种常用的代数重构算法^[12]。

1)加法代数重构算法：

式中，k表示第k次迭代；x^k为第k次迭代的结果；Aⁱ表示矩阵A的第i行；〈·〉表示向量内积；为松弛参数向量；γ₀为松弛因子。

2)乘法代数重构算法(multiplicative ART，MART)：

式中，x_j^k表示第j个网格的参数值x_j第k次迭代的结果；a_i，j为矩阵A第i行、第j列上的元素值；其他参数的含义与式(4)相同。MART方法的迭代步骤与ART方法基本相同，但其迭代改进以乘法的形式进行，其优点是可以保证反演结果为正值。

3)同时迭代重构算法(simultaneous iterations reconstruction technique，SIRT)：

式中，参数的含义与式(4)和式(5)相同。

从式(4)、式(5)可以看出，ART和MART的重构结果在一定程度上取决于所采用的观测值的顺序，而SIRT则避开了这种影响。但研究结果显示，使用不同次序的观测值所引起的影响相比其他参数(如初始场X₀的选取和迭代次数等)的影响是可以忽略的^[12]。在GNSS水汽参数的迭代重构中，通常采用式(7)中的δ、σ参数作为收敛条件，在迭代收敛的情况下，当其小于一定的阈值时，就把迭代的结果作为最终的重构结果，一般情况下通过100~200次的迭代便可收敛。

3 代数重构算法在水汽层析中的应用

利用模拟试验对代数重构算法在水汽层析解算中的可行性进行分析。截至2011年，武汉暴雨研究所长江中游暴雨外场观测试验基地的GNSS水汽监测网已建成58个测站，采样时间间隔为30 s，全天24 h连续观测，其中位于鄂东区域的22个测站分布较为均匀，测站平均间距在30 km左右，因此，该区域被选定为层析的试验区域。层析区域测站的分布和水平网格划分情况如图 1所示。在大地高方向上，层析区域设定为0~10 km，分辨率为1 km。因此，在该层析试验中，共需解算的湿折射率参数为7×4×10=280个。

利用2012年5月份某天无线电探空仪的观测结果模拟出每个网格的湿折射率值，结合卫星的实际位置和层析区域的网格划分情况，模拟出接收机至卫星连线方向上的SWD(卫星截止高度角取为5°)，在此基础上加入一定的观测误差，生成用于水汽层析的观测值SWD。关于迭代初值X₀，我们选取了利用武汉多年的无线电探空数据(探空站位置如图 1中黑色圆圈所示)获得的5月份水汽垂直分布的统计结果^[14]。

图  1  层析区域站点分布及水平网格划分

Figure  1.  Distribution of GNSS Sites and Horizontal Grids Defined in Tomography Region

下载: 全尺寸图片 幻灯片

3.1   松弛向量的选取

代数重构算法的基本原理是基于〈Aⁱ, x^k〉与实际观测值SWD的差值(或比值)，按照一定的方式计算出改正值，用于对第k次迭代结果x^k进行修正，直到满足某一条件退出迭代为止。松弛参数λ则用来调整改正值占差值(或比值)的比重。下面将利用ART方法和MART方法进行水汽湿折射率参数的重构。为了进一步了解松弛参数向量λ在迭代中所起的作用，将ART方法中的λ调整为λ=一并进行解算，该迭代方法记为IART。

图 2给出了利用ART、IART、MART等三种方法得到的重构结果与真值的差值序列，图中x轴表示网格号。根据图 1给出的网格划分情况，第一层网格号为1~28，第二层为29~56，以此类推；由于有的网格没有射线穿过，因此其迭代结果维持在迭代初值未发生变化，在图 2中没有给出其误差值。从图 2可以看出，利用三种方法得到的水汽湿折射率参数的重构结果在第一层中相对误差的绝对值均在10%左右；在较低层外的其他层，IART方法均可较好地重构水汽湿折射率场，而采用ART方法得到的重构结果则在较高层具有较大的误差。从ART和IART的松弛参数向量λ的比较可以看出，ART方法为每个网格确定的改正值，与射线穿过网格的长度成正比，即对某一条特定的射线来说，它穿过网格的长度越大，则对该网格参数的改正值就越大；但是，对流层大气湿折射率一般随着高度的增加呈指数递减，因此，利用射线穿过各个网格的长度来确定改正值的大小明显是不合理的，而IART方法考虑了大气湿折射率的这一变化规律，利用该方法可明显降低高层网格的重构误差。另外，MART方法采用乘法的形式进行迭代修正，松弛参数向量位于指数的位置上，其重构结果与IART方法的结果相当。为了更加直观地显示出重构结果，图 3给出了某网格上垂直方向上三种迭代方法的重构值与真值的比对结果。

图  2  利用几种代数重构算法重构的大气湿折射率N_w的误差

Figure  2.  Reconstruction Error of the Wet Refractivity Based on ART, IART and MART

下载: 全尺寸图片 幻灯片

图  3  重构的大气湿折射率与网格真值的比较

Figure  3.  Comparison Between the Reconstructed Field of Wet Refractivity and the Truth Value

下载: 全尺寸图片 幻灯片

表 1为三种代数重构算法与附加约束条件的层析解算方法的误差的统计结果。从表 1中可以看出，利用IART方法和MART方法得到的重构结果的精度相当，均明显优于利用ART方法得到的结果，但略低于附加约束条件的层析解算方法的解算精度，这可能是由于附加约束条件的层析解算方法采用了网格间水汽参数的平滑约束，增加了观测信息。

表  1  不同解算方法重构误差的均值、均方根值

Table  1.  Mean and RMS Value of the Reconstructed Error Based on Different Reconstruction Technique

方法	松弛因子	均值/(mm\5km-1)			均方根值/(mm\5km-1)
		所有网格	有射线穿过的网格		所有网格	有射线穿过的网格
ART	约0.2	-0.50	0.39		8.40	7.47
IART	约0.1	-1.16	-0.26		5.64	4.16
MART	约0.15	-1.28	-0.32		5.78	4.24
附加约束条件的层析解算方法	/	-0.99	-0.25		4.74	3.18

下载: 导出CSV 

| 显示表格

3.2   观测值误差对重构结果的影响分析

GNSS水汽层析的重构结果的质量受多种因素的影响，如层析算法、先验约束信息及观测值的质量等。为了分析在不同的算法中观测值的误差对最终的层析结果的影响，我们在模拟的观测值真值的基础上，分别加入了均值为0、标准差为1~50 mm的观测噪声，并利用IART、MART及附加约束条件的层析解算方法进行解算。

图 4为观测值误差与层析结果误差的对比关系图，从图 4中可以看出，对于IART和MART来说，在观测值误差标准差小于20 mm时，重构结果的误差较为稳定，随着观测值误差的增大，层析结果的误差基本上也随之增大，但MART方法得到的层析结果的误差增加的速度明显低于IART方法；而附加约束条件的层析解算方法的结果误差除了在开始时略微增大外，随着观测值误差的增大，其精度较为稳定(这可能是由于该方法中采用了额外的约束信息的缘故，在实际的解算中应该合理选取观测权阵，否则容易造成约束过紧而产生虚假的重构湿折射率场^[14])。在实际的业务工作中，考虑到代数重构算法的解算稳定性，当采用该方法进行水汽层析解算时，应尽量获取较高精度的观测值信息(如采用精度较高的映射函数VMF1，合理处理观测残差以得到其中的对流层延迟残差等)，同时，采用MART方法可获得相对精度较高的层析结果。

图  4  不同观测值误差下重构的大气湿折射率误差的均方根值

Figure  4.  RMS of the Reconstructed Wet Refractivity Field at Different Observation Error Levels

下载: 全尺寸图片 幻灯片

4 结语

本文基于鄂东局域GNSS气象监测网，考虑到测站数目及空中可视卫星数目增多的现状，对代数重构算法在GNSS水汽层析解算中的应用进行了研究，得出如下结论。

传统的ART方法在GNSS水汽层析解算中，会出现较大的重构误差，尤其在较高层上，几乎无法获得较为准确的水汽参数解；针对对流层大气湿折射率随高度的指数分布情况，本文调整了松弛参数向量。结果表明，利用IART方法可以较好地重构水汽的垂直分布；而利用MART方法得到的水汽层析结果和IART方法的重构精度相当。

与附加约束条件的层析解算方法相比，IART和MART方法更易受到观测值误差的影响。利用IART和MART方法得到的重构结果在观测值误差较小的情况下较为稳定，在误差标准差大于20 mm的情况下，随着观测值误差的增加而增大，但MART方法的结果优于IART方法的结果。在实际中考虑采用代数重构算法的情况下，应提高水汽层析观测值的精度以便得到较好的层析结果。