## 2 Tomography Theory and Method  
### Establishment of Tropospheric Tomography Modeling  
It is well known that the zenith total delay (ZTD) is an average parameter in spatial aspects, and some important information about the local variability of the atmosphere is lost. In fact, ZTD is an integral parameter, and the information about the vertical structure of the atmosphere is lost [20]. The zenith wet delay (ZWD) can be extracted from ZTD by removing the zenith hydrostatic delay (ZHD) calculated using the Saastamoinen model [21]. Therefore, the SWD (Slant wet delay) can be obtained by projecting the ZWD value to the slant path, which is regarded as the input information of GNSS troposphere tomography modeling. According to the distances of GNSS rays calculated in each divided voxel of the tomography area and the unknown wet refractivity to be estimated, the observation equation of troposphere tomography can be expressed as:  
\[SWD=\sum d_{ijk}\cdot x_{ijk}, \tag{1}\]  
where the subscripts of \((i,j,k)\) represent the location of voxels in horizontal and vertical directions. \(d\) is the distance of satellite ray in the voxel \((i,j,k)\), while \(x\) is the unknown parameters to be estimated in the voxel \((i,j,k)\). As mentioned before, some empirical constraints are considered to overcome the rank deficiency problem. For the small tomography area, a certain voxel can be described by the weighted value of its neighbors horizontally [22]. The relative horizontal constraint can be unified expressed as follows:  
\[h_{1}\cdot x_{1}+h_{2}\cdot x_{2}+\cdots+h_{m-1}\cdot x_{m-1}-x_{m}+h_{m+1} \cdot x_{m+1}+\cdots+h_{n}\cdot x_{n}=0, \tag{2}\]  
where \(h\) is the horizontal weighted coefficient while \(n\) is the total number of voxels in the same layer. The vertical constraint is constructed based on the exponential function in different layers as follows [23; 24]:  
\[x_{i}=x_{i-1}\cdot e^{(h_{i-1}-h_{i})/H}, \tag{3}\]where \(x_{i}\) and \(h_{i}\) are the wet refractivity value and the corresponding height of the voxel \(i\), respectively. \(H\) denotes the water vapor scale height, which ranges from 1 to 2 km. Aside from horizontal and vertical constraints, an a priori constraint is also imposed for the location of radiosonde station using the mean value of radiosonde data in the first three days before the tomography epoch. The priori constraint can be described as follows:  
\[x_{q}=x_{0_{q}}^{rs}, \tag{4}\]  
where \(x_{0}^{rs}\) represents the initial wet refractivity value calculated using the radiosonde data in the first three days before the tomography epoch. \(q\) is equal to the number of vertical layers of the tomography area. Due to the data from multi-GNSS (GPS, GLONASS and BDS) being used in our study, the final tomography modeling can thus be expressed as follows:  
\[ \begin{bmatrix}\mathbf{y}_{m1\times 1}^{gps}\\ \mathbf{y}_{m2\times 1}^{gloans}\\ \mathbf{y}_{m3\times 1}^{bds}\\ \mathbf{y}_{m4\times n1}^{bds}\\ \mathbf{y}_{m5\times n1}^{b}\\ \mathbf{y}_{m6\times 1}^{b}\end{bmatrix} = \begin{bmatrix}\mathbf{A}_{m1 \times 1}^{gps}\\ \mathbf{A}_{m2\times 1}^{gloans}\\ \mathbf{A}_{m3\times n1}^{bds}\\ \mathbf{A}_{m4\times n1}^{bds}\\ \mathbf{A}_{m5\times n1}^{bV}\\ \mathbf{A}_{m6\times 1}^{p}\end{bmatrix} \cdot\mathbf{X}_{n1\times 1}+ \begin{bmatrix}\mathbf{\Delta}_{m1\times 1}^{gps}\\ \mathbf{\Delta}_{m2\times 1}^{gloans}\\ \mathbf{\Delta}_{m3\times 1}^{bds}\\ \mathbf{\Delta}_{m4\times 1}^{bds}\\ \mathbf{\Delta}_{m5\times 1}^{V}\\ \mathbf{\Delta}_{m6\times 1}^{p}\end{bmatrix} , \tag{5}\]  
where \(\mathbf{y}^{tp1}\), \(\mathbf{A}^{tp1}\) and \(\mathbf{\Delta}^{tp1}\) (\(tp1=gps,gloans,bds\)) represent the column vectors with a set of SWD observations derived from multi-GNSS (GPS, GLONASS and BDS), the design matrices with the value of the distance of the satellite ray inside a certain voxel and the noise of the multi-GNSS observation equations, respectively.