# A New GNSS-Derived Water Vapor Tomography Method Based on Optimized Voxel for Large GNSS Network
## 2 Data and Methods
### Conventional Solution for WV Tomography  
As shown in Figure 1, owing to the unsatisfactory geometric distribution of satellite constellations and ground receivers, inevitably, some voxels are not passed by any signal. This demonstrates that Equation (7), is ill-posed. Moreover, the problem becomes more significant for sparser receivers ([PERSON], 2016). Thus, some constraints are required to stabilize Equation (7). The WV density in a certain voxel can be determined by its surrounding voxels ([PERSON], 2016; [PERSON], 2016). Based on the above theory, the general method for WV tomographic equation solution introduces the following three types of constraints.  
(1) Horizontal Constraint  
These constraints are used to determine the WV value in null-ray passed voxels according to their neighbors in the same level [18], which can be described:  
\[H\times X\ =\ 0, \tag{8}\]  
where \(H\) is the constraint matrix for the horizontal direction.  
(2) Vertical Constraint  
The relation between the WV densities of two vertically adjacent voxels can be expressed by the vertical constraint expressed as:  
\[V\times X\ =\ 0, \tag{9}\]  
where \(V\) stands for vertical constraint matrix.  
(3) Boundary Constraint  
The tropospheric WV density rapidly decreases with the increase in height. In this case, we assume that the WV density in the top-layer voxels is zero. Accordingly, in addition to the horizontal and vertical constraints, we add the following top-level boundary constraints,  
\[B\times X\ =\ 0. \tag{10}\]  
Considering the above three types of constraint equations, the coefficient matrix \(A\) and result matrix \(d\) in Equation (7) (i.e., original tomographic equation) are modified to  
\[A_{new}\ =\ \left[\begin{array}{c}A_{ori}\\ H\\ V\\ B\end{array}\right]\ d\ =\ \left[\begin{array}{c}d_{ori}\\ 0\\ 0\\ 0\end{array}\right]. \tag{11}\]  
However, even with the expanded coefficient matrix \(A\) and \(d\) (Equation (11)), an ill-posed problem exists in Equation (7) owing to the limited effective observations in one epoch. To overcome this  
Figure 1: Principle of water vapor (WV) tomography with Global Navigation Satellite System (GNSS) observations.  
problem, we multiply the rows of matrix \(A_{ori}\), which is composed of original GNSS observations, by setting a \"tomography window.\" The tomographic window is based on the assumption that the WV content is stationary in a period of time and can significantly expand matrix \(A_{ori}\) through multi-epoch observations.  
After the above process steps, the solution of the updated Equation (7), i.e., the value of \(x\), expresses the WV density in each voxel of the tomographic region. However, considering that the coefficient matrix composed of observations and constraints is a large-scale sparse matrix with the corresponding inverse problem, it is challenging to obtain \(x\). Two main methods can be used to solve the tomographic equations, including non-iterative and iterative algorithms [17]. The latter exhibits a higher performance than that of the former. Among the iterative methods, the multiplicative algebraic reconstruction technique algorithm [17] and Kalman filtering (KF) technique [34; 35] have been widely used. The KF technique includes an advanced algorithm, simple implementation, and good performances in numerous fields. Thus, we use KF to solve the tomographic equations. Besides, references [36; 37] describe the KF for tomographic use in detail.  
Overall, the general tomographic approach introduces auxiliary constraints (horizontal, vertical, and boundary constraints) into the original tomographic Equation (7) to solve its ill-posed problem and obtain its solutions through an iterative algorithm (such as the KF in this case).