Combining CHAMP and Swarm Satellite Data to Invert the Lithospheric Magnetic Field in the Tibetan Plateau

CHAMP and Swarm satellite magnetic data are combined to establish the lithospheric magnetic field over the Tibetan Plateau at satellite altitude by using zonal revised spherical cap harmonic analysis (R-SCHA). These data are integrated with geological structures data to analyze the relationship between magnetic anomaly signals and large-scale geological tectonic over the Tibetan Plateau and to explore the active tectonic region based on the angle of the magnetic anomaly. Results show that the model fitting error is small for a layer 250–500 km high, and the RMSE of the horizontal and radial geomagnetic components is better than 0.3 nT. The proposed model can accurately describe medium- to long-scale lithospheric magnetic anomalies. Analysis indicates that a negative magnetic anomaly in the Tibetan Plateau significantly differs with a positive magnetic anomaly in the surrounding area, and the boundary of the positive and negative regions is generally consistent with the geological tectonic boundary in the plateau region. Significant differences exist between the basement structures of the hinterland of the plateau and the surrounding area. The magnetic anomaly in the Central and Western Tibetan Plateau shows an east–west trend, which is identical to the direction of the geological structures. The magnetic anomaly in the eastern part is arc-shaped and extends along the northeast direction. Its direction is significantly different from the trend of the geological structures. The strongest negative anomaly is located in the Himalaya block, with a central strength of up to −9 nT at a height of 300 km. The presence of a strong negative anomaly implies that the Curie isotherm in this area is relatively shallow and deep geological tectonic activity may exist.


Introduction
The lithospheric magnetic field accounts for approximately 4% of the total energy of the earth's magnetic field [1]. It originates from the crust and upper mantle of magnetic rocks and reflects a 3D spatial distribution of rock magnetism. Given its different rock magnetic environments, temperature pressure conditions, magnetic carriers, and tectonic evolutions, the lithospheric magnetic field carries abundant and complex geological structures and tectonic movement information [2,3]. Therefore, analyzing the magnetic distribution frame in the depth of the crust and the top of the upper mantle using a high-resolution and high-accuracy lithospheric magnetic field model can reveal important information regarding regional tectonics and their evolution, the division of rock masses, and the heat flow of the lithosphere [4][5][6][7][8].
direction. Common regional modeling methods mainly include the equivalent source method [25], rectangular harmonic analysis [26,27], the spherical cap harmonic/revised spherical cap harmonic analysis (R-SCHA) method [28][29][30][31], and the 3D Taylor/Legendre polynomial method [32][33][34]. These methods have various advantages and disadvantages. Using CHAMP and Swarm satellite magnetic data, the current study utilizes Sq local time conditions, Dst/Kp geomagnetic index conditions, Hamming along-track high-pass filtering, and the hierarchical gridding method, and combines EMM2015, MF7 and other auxiliary models to thoroughly process satellite magnetic data. The present paper adopts zonal revised spherical cap harmonic analysis to inverse the lithospheric magnetic field over the Tibetan Plateau and analyzes the relationship between magnetic anomalies and geological structures to explore the active region of geological structures.

Mathematical Model
The model is established in a source-free conical domain denoted by Ω, which is bounded by a finite cone semi-angle of θ 0 . The domain lower and upper boundaries along the radial r are defined with two spherical surfaces of radii a and b, respectively, as shown in Figure 1. The boundary of cone Ω was defined as follows: ∂Ω = ∂ θ 0 Ω ∪ ∂ a Ω ∪ ∂ b Ω. ∂ θ 0 Ω : {θ = θ 0 , a ≤ r ≤ b}, ∂ a Ω : {θ ≤ θ 0 , r = a} and ∂ b Ω : {θ ≤ θ 0 , r = b} represent the lateral boundary, the lower and upper boundary of Ω, respectively [18,30]. Suppose that the geomagnetic field meets the irrotational and solenoidal requirements in the closed region Ω. The geomagnetic field scalar potential meets the Laplace equation ∇ 2 V = 0 under the condition of boundary constraint. By using the separation of variables method, the solution of Laplace's equation can be written as a sum of infinite series in a cone coordinate system [30]: ∑ m≥0 a r n k +1 G i,m n k cos(mϕ) + H i,m n k sin(mϕ) · P m n k (cos θ) r a n k G e,m n k cos(mϕ) + H e,m n k sin(mϕ) · P m n k (cos θ) R p (r) G m p cos(mϕ) + H m p sin(mϕ) · K m p (cos θ) R 0 G m 0 cos(mϕ) + H m 0 sin(mϕ) · P m 0 (θ) (1) The corresponding expressions for the geomagnetic component are as follows [18,31]: (n k + 1) a r n k +2 P m n k G i,m n k cos(mϕ) + H i,m n k sin(mϕ) n k r a n k −1 P m n k G e,m n k cos(mϕ) + H e,m n k sin(mϕ) dr K m p G m p cos(mϕ) + H m p sin(mϕ) (2) In Equation (2), X, Y, Z represent the north, east, and vertical downward components of the geomagnetic field at an observation station, respectively. a denotes the mean radius of the earth, which is taken to be 6371.2 km. ϕ, θ, r represent the longitude, geocentric colatitude, and geocentric distance, respectively, in the cone coordinate system. P m n k (cos θ) is the associated Schmidt quasi-normalization Legendre functions of real degree n k and integer order m. K m p (cos θ) is the Mehler functions. P m 0 (cos θ) is the particular form of the K m p (cos θ) function when n ≤ m and n = 0. R p (r) are the radial functions.
is the particular form of the R p (r) function when n = 0 [35]. (G i,m n k , H i,m n k , G e,m n k , H e,m n k ) are the Legendre coefficients, and (G m p , H m p , G m 0 , H m 0 ) are the Mehler coefficients. k, p are integers that correspond to the sequence number of roots n k , n p Superscripts i, e of the Legendre coefficient refer to the components of the internal and external source field. Achieving the separation of the internal and external field sources is more difficult when the R-SCHA algorithm is used with the aid of the boundary conditions of cone Ω compared with spherical harmonic analysis. Therefore, careful selection and model correction are conducted on the satellite data to exclude the interference of external source fields as much as possible. The component of the external source field introduced in the algorithm represents neither the ionosphere nor the magnetic field. Instead, it adopts two groups of Legendre coefficients to fit the distribution form of the geomagnetic field in cone Ω. According to radial function form (r/a) corresponding to the Legendre coefficient (G e,m n k , H e,m n k ) in the formula (2), we can infer that (G e,m n k , H e,m n k ) mainly describes the contribution of the geomagnetic field in the middle and upper space altitude layers and several long wavelength components of the geomagnetic field in the earth's surface space. After testing, Torta et al. found that achieving the best fit is difficult using the model if the component of the so-called external source field was removed from the algorithm, but the independent component of the internal source field was only used even if the observation data was distributed in the lower bottom boundary of the spherical cap. The fitting effect improves significantly if the component of the external source field is added to the Legendre function [36,37].

Inversion Method
The observation vector of the field L = {l 1 , l 2 , · · · l n } and parameter vector to be estimated (revised spherical cap harmonic coefficient) X = {x 1 , x 2 , · · · , x k }, n >> k, satisfy the following linear relation: In Equation (3), ∆ represents the observation noise. Its stochastic characteristic is E(∆) = 0, D(∆) = σ 2 0 P L −1 where P L represents the observation vector weight matrix and σ 0 denotes the mean square error of unit weight. A target function is constructed according to the estimation criterion of the least square principle of compensation to provide a unique stable resolution for Equation (3): Parameter vectorX that meets the minimization of Equation (4) is the solution of the linear model by Equation (3). Ω(X) is called a stable function. It functions to transfer ill-posed problems to well-posed problems. α represents the regularization parameter (smoothing factor). It functions to balance the weights of two items of target function J α (L, X). · 2 P L represents the weighted Euclidean two-norm [38]. In practice, stable functional Ω(X) typically uses the Tikhonov regularization method with the following form [39][40][41][42]: Ω(X) = X 2 P X = X T P X X.
In Equation (5), P X represents the parameter weight matrix. For the above stable functional, the minimum value of target function J α (L, X) is calculated as follows: If ∂J α (L, X)/∂X = 0, then a parametric solution can be obtained as follows: Based on the abovementioned analysis, the selection of reasonable parameter weight function P X is the key to solving ill-posed problems. Given that the Legendre and Mehler functions are mutually orthogonal, parameter weight function P X can be set as a diagonal matrix. The elements in the diagonal line of matrix P X can be determined by geomagnetic energy per coefficient. Besides the additional constraint or a stable functional Ω(X), regularization parameter α also plays an important role in ill-posed problems [43]. α balances the influence of data fitting degree and regularized error for parameter estimation [43]. The optimal regularization parameter should minimize the comprehensive error caused by both factors. When α is excessively small, the inversion model tends to fit the observation data effectively and with a small constraint effect, which may induce over-fitting of high-frequent noise in the observation data. As a result, the solution of target function loses physical meaning because of shock. If α is excessively large, then the inversion model will be inclined to constraints, thereby leading to excessive smoothing of the model parameters. As a result, high-frequency information is lost, and the fitting effect of observation data is poor.
To estimate the optimal regularization parameter a, Thébault et al. used CM4 synthetic data to calculate the correlation coefficient between the R-SCHA prediction and the expected magnetic field and in the data gap layer at a height of 100~300 km. After testing, set α = 1e −6 was identified for the internal Legender functions. Set α = 9e −7 was determined for the Mehler functions. The R-SCHA algorithm can provide the best balance between the fit and the smoothness of the edges, as well as the reliability of the data gap layer [31]. By using multiple simulation data, Thébault and Gaya-Piqué (2008) found that the model can also provide a stable space continuation in equal weight without considering regularization conditions. From a side view, when the truncation level of a model is reasonably designed, the R-SCHA algorithm is stable for local regional modeling without obvious shock effects [44].

Coordinate Transformation
Before model inversion, the geographical coordinates (ϕ, θ, r) S and geomagnetic component (N, E, U) S of the observation data must be transformed from the spherical coordinate reference frame to the cone coordinate reference frame [45], as shown in Figure 1. inversion model tends to fit the observation data effectively and with a small constraint effect, which may induce over-fitting of high-frequent noise in the observation data. As a result, the solution of target function loses physical meaning because of shock. If  is excessively large, then the inversion model will be inclined to constraints, thereby leading to excessive smoothing of the model parameters. As a result, high-frequency information is lost, and the fitting effect of observation data is poor.
To estimate the optimal regularization parameter a , Thébault et al. used CM4 synthetic data to calculate the correlation coefficient between the R-SCHA prediction and the expected magnetic field and in the data gap layer at a height of 100~300 km. After testing, set  =9e  was determined for the Mehler functions. The R-SCHA algorithm can provide the best balance between the fit and the smoothness of the edges, as well as the reliability of the data gap layer [31]. By using multiple simulation data, Thébault and Gaya-Piqué (2008) found that the model can also provide a stable space continuation in equal weight without considering regularization conditions. From a side view, when the truncation level of a model is reasonably designed, the R-SCHA algorithm is stable for local regional modeling without obvious shock effects [44].

Coordinate Transformation
Before model inversion, the geographical coordinates ( , , ) S r   and geomagnetic component ( , , ) S N E U of the observation data must be transformed from the spherical coordinate reference frame to the cone coordinate reference frame [45], as shown in Figure 1.

Transformation of Geographical Coordinates
The geographical coordinate of the observation data (ϕ, θ, r) S is transformed from the spherical coordinate system to the cone coordinate system. The transformation is conducted in three steps: (1) In the spherical coordinate reference frame, Equation (8) is used to transform the geographical coordinate (ϕ, θ, r) S to the earth's core space rectangular coordinate system (X, Y, Z) S ; (2) Equation (9) is utilized to transform (X, Y, Z) S to the cone space rectangular coordinate In Equation (9), (ϕ 0 , θ 0 ) refers to the longitude and geocentric colatitudes of the cone pole; (3) In the cone coordinate reference frame, Equation (10) is used to transform (X, Y, Z) C to the cone geographical coordinate system (ϕ, θ, r) C :

Transformation of Geomagnetic Observation Component
The transformation of geomagnetic observation component (N, E, U) S is conducted through three steps: (1) In the spherical coordinate reference frame, Equation (11) is used to transform the geomagnetic observation component (N, E, U) S in station center coordinate system to the geocentric space In Equation (11), (B 0 , L 0 ) represent the geocentric latitude and longitude in the spherical coordinate system, respectively.
(2) Similarly, Equation (12) is used to transform (X, Y, Z) S−VAL into the cone space rectangular coordinate system (X, Y, Z) C−VAL , namely: (3) In the cone coordinate reference frame, Equation (13) is used to transform (X, Y, Z) C−VAL to geomagnetic observation component (N, E, U) C in station center coordinate system. (13) (Note: in Equation (13), (N, E, U) C correspond to the geomagnetic observation component (X, Y, Z) in the spherical cap harmonic model in Equation (2)).

Data Source
The satellite data used in the present paper include: (1) CHAMP vector observation data from 2008.1~2010.9. The data in this period correspond to low solar activity years with a corresponding satellite orbital altitude of 250~340 km, as shown in Figure 2. Therefore, the data are less affected by the external magnetic field. This condition is conducive to the exploration of middle and small-scale magnetic anomaly signals. (2) Given that the Swarm C satellite orbit is excessively high to be sensitive to middle and small-scale magnetic anomaly signals, the vector observation data of Swarm A and B satellites from 2014.1~2015.12 with an orbital altitude of 450~510 km are selected. Satellite circular half-orbit separation and data are subsampled every 5s before data processing. Track separation is conducted along-track filtering corrections of observations. At satellite orbit altitude, sampling with 5 s intervals corresponds to a spacing of approximately 35 km. The sampling can weaken the reciprocal mixing of small-scale magnetic field signals and effectively control the anisotropic error to be introduced in the final model.
X Y Z in the spherical cap harmonic model in Equation (2)).

Data Source
The satellite data used in the present paper include: (1) CHAMP vector observation data from 2008.1~2010.9. The data in this period correspond to low solar activity years with a corresponding satellite orbital altitude of 250~340 km, as shown in Figure 2. Therefore, the data are less affected by the external magnetic field. This condition is conducive to the exploration of middle and small-scale magnetic anomaly signals. (2) Given that the Swarm C satellite orbit is excessively high to be sensitive to middle and small-scale magnetic anomaly signals, the vector observation data of Swarm A and B satellites from 2014.1~2015.12 with an orbital altitude of 450~510 km are selected. Satellite circular half-orbit separation and data are subsampled every 5s before data processing. Track separation is conducted along-track filtering corrections of observations. At satellite orbit altitude, sampling with 5 s intervals corresponds to a spacing of approximately 35 km. The sampling can weaken the reciprocal mixing of small-scale magnetic field signals and effectively control the anisotropic error to be introduced in the final model.

Data Screening
To decrease the interference of the external source field in establishing lithospheric magnetic field model, according to the satellite data screening criteria proposed by Maus and Sabaka et al. [24,46], the data are selected in local time 22:00~03:00 to minimize the contributions from the ionospheric Sq field. The s D t and Kp index conditions in Equation (14) are utilized to select the data corresponding to a quiet period of geomagnetic activity to weaken the effects of the equatorial ring current and magnetospheric convection intensity:

Data Screening
To decrease the interference of the external source field in establishing lithospheric magnetic field model, according to the satellite data screening criteria proposed by Maus and Sabaka et al. [24,46], the data are selected in local time 22:00~03:00 to minimize the contributions from the ionospheric Sq field. The Dst and Kp index conditions in Equation (14) are utilized to select the data corresponding to a quiet period of geomagnetic activity to weaken the effects of the equatorial ring current and magnetospheric convection intensity: |Dst| ≤ 20nT and d|Dst| ≤ 10nT Kp ≤ 2 and d|Kp|≤ 2 (14) d|Dst| and d|Kp| represent changes over the previous three hours. A selection condition is selected using quality flags indicators, such as spacecraft maneuvers and the on/off status of the star camera, to judge the running state of a star camera. Finally, short-arc tracks without data beyond 240 epochs are eliminated to fully guarantee the quality of the observed data. Figure 3 shows the data preprocessing procedure of extracting a lithospheric magnetic field signal based on satellite magnetic survey data.
Sensors 2017, 17, 238 9 of 18 camera, to judge the running state of a star camera. Finally, short-arc tracks without data beyond 240 epochs are eliminated to fully guarantee the quality of the observed data. Figure 3 shows the data preprocessing procedure of extracting a lithospheric magnetic field signal based on satellite magnetic survey data.

Main Magnetic Field and External Magnetic Field Correction
The non-lithospheric magnetic field must be corrected for the selected data to capture lithospheric magnetic field signal [2]. The corrected magnetic fields cover the main magnetic field, magnetosphere, and ionospheric magnetic field, as well their induced counterparts. Considering the weakness of the ionospheric magnetic field at night, corrections are not conducted using the model. To correct the magnetospheric sources, the long wavelength signals are separated by means of along-track filtering in Section 3.4. To correct the main magnetic field, the POMME-10, EMM2015, and MCO_SHA_2D models were tested [47]. Their mutual differences are small. Therefore, the EMM2015 model (1-15 degrees) is selected as the earth's main magnetic field and its secular variation. The model is removed from the satellite observation value to obtain the so-called residual field. Figure 4 describes the radial components of CHAMP and Swarm satellite observations after deleting the main magnetic field. The MF7 lithospheric magnetic field model is eliminated to comprehensively analyze the along-track filtering data in Section 3.4. Figure 4 shows that before making along-track filtering corrections, the residual magnetic field of CHAMP and Swarm observations range from −10~10 nT and 0~30 nT, respectively.

Main Magnetic Field and External Magnetic Field Correction
The non-lithospheric magnetic field must be corrected for the selected data to capture lithospheric magnetic field signal [2]. The corrected magnetic fields cover the main magnetic field, magnetosphere, and ionospheric magnetic field, as well their induced counterparts. Considering the weakness of the ionospheric magnetic field at night, corrections are not conducted using the model. To correct the magnetospheric sources, the long wavelength signals are separated by means of along-track filtering in Section 3.4. To correct the main magnetic field, the POMME-10, EMM2015, and MCO_SHA_2D models were tested [47]. Their mutual differences are small. Therefore, the EMM2015 model (1-15 degrees) is selected as the earth's main magnetic field and its secular variation. The model is removed from the satellite observation value to obtain the so-called residual field. Figure 4 describes the radial components of CHAMP and Swarm satellite observations after deleting the main magnetic field. The MF7 lithospheric magnetic field model is eliminated to comprehensively analyze the along-track filtering data in Section 3.4. Figure 4 shows that before making along-track filtering corrections, the residual magnetic field of CHAMP and Swarm observations range from −10~10 nT and 0~30 nT, respectively.

Along-Track Filtering
After data screening and main magnetic field correction, several non-lithospheric signals that were not modeled remain in the observations. These residual signals are dominated by long-wave components of magnetospheric contributions [20,48]. Another important source of residuals is the ionospheric disturbance magnetic field, and the sources of secondary importance are the induced magnetic field and the secular variation error from the main magnetic field. Along-track high-pass filtering of fast Fourier transform is adopted in the current study to eliminate the large-scale residual signals. The challenge in employing this method is the selection of the cut-off wavelength. The magnetic field contributions from wavelength signals greater than 1200 km are reduced with high pass filtering by using a Hamming window function (filter of order 128). In the study, a lithospheric field (MF7) is removed before performing along-track filtering for the residual signals in the actual filtering process to lower the level of spectral leakage [49]. Finally, a residual analysis is conducted for the values after filtering, and magnetic anomaly values greater than ±25 nT are eliminated. At satellite altitude, the intensity of lithospheric magnetic anomaly is generally less than 25 nT [50].

Data Hierarchical Gridding
Hierarchical gridding of the data is applied to avoid the shock effect of the model in the horizontal and vertical directions caused by the uneven data distributions. In the horizontal direction, the grid resolution of data is set to 0.5° × 0.5°. In the vertical direction, the layered interval of data is set to 10 km. A layer of insufficient data is merged into an adjacent layer. For the data distributed in the same grid unit, the data nearest the average value of the grid is selected as the magnetic anomaly value of the grid to ensure the equilibrium of the data in the horizontal and vertical directions.

Related Parameter Setting
To achieve the best fit between the selected spherical cap and the modeling region, the pole of the spherical cap is generally set at the center of the region or nearby. The geographic latitude and

Along-Track Filtering
After data screening and main magnetic field correction, several non-lithospheric signals that were not modeled remain in the observations. These residual signals are dominated by long-wave components of magnetospheric contributions [20,48]. Another important source of residuals is the ionospheric disturbance magnetic field, and the sources of secondary importance are the induced magnetic field and the secular variation error from the main magnetic field. Along-track high-pass filtering of fast Fourier transform is adopted in the current study to eliminate the large-scale residual signals. The challenge in employing this method is the selection of the cut-off wavelength. The magnetic field contributions from wavelength signals greater than 1200 km are reduced with high pass filtering by using a Hamming window function (filter of order 128). In the study, a lithospheric field (MF7) is removed before performing along-track filtering for the residual signals in the actual filtering process to lower the level of spectral leakage [49]. Finally, a residual analysis is conducted for the values after filtering, and magnetic anomaly values greater than ±25 nT are eliminated. At satellite altitude, the intensity of lithospheric magnetic anomaly is generally less than 25 nT [50].

Data Hierarchical Gridding
Hierarchical gridding of the data is applied to avoid the shock effect of the model in the horizontal and vertical directions caused by the uneven data distributions. In the horizontal direction, the grid resolution of data is set to 0.5 • × 0.5 • . In the vertical direction, the layered interval of data is set to 10 km. A layer of insufficient data is merged into an adjacent layer. For the data distributed in the same grid unit, the data nearest the average value of the grid is selected as the magnetic anomaly value of the grid to ensure the equilibrium of the data in the horizontal and vertical directions.

Related Parameter Setting
To achieve the best fit between the selected spherical cap and the modeling region, the pole of the spherical cap is generally set at the center of the region or nearby. The geographic latitude and longitude of the Tibetan Plateau are 26.0 • N ∼ 40.0 • N, 73.0 • E ∼ 105.0 • E, respectively. The difference of longitude is approximately twice the difference of latitude. Under a certain level of truncation, partition modeling of a small spherical cap can be used to improve the accuracy of the model in the whole region and generate good spatial resolution [51]. Therefore, the Tibetan Plateau was divided into two subregions, and each subregion was independently modeled. The pole of the spherical cap in Region 1 is located at (33 • N, 81 • E), and the pole of spherical cap in Region 2 is located at (33 • N, 97 • E). To suppress the boundary effect, the half aperture of the spherical cap was amplified of 1 • , and θ 1 = θ 2 = 10 • was set. The largest overlap between adjacent spherical caps is approximately 7.8 • . Thus, all the observed data are located within the range of the spherical cap.
The mean altitude of the Tibetan Plateau is greater than 4 km. The height of the lowest orbit of the selected CHAMP satellite data is approximately 250 km. The height of the highest orbit of Swarm satellite data is approximately 510 km. Therefore, the geocentric distances of the lower and upper boundaries of the cone are defined as a = Re + 240 km and b = Re + 520 km, where Re = 6371.2 km refers to the earth's mean radius. The related parameters are listed in Table 1. In Table 1, K i max , K e max represent the maximum truncation degree of the Legendre function in the internal and external source fields, respectively. P max denotes the maximum truncation degree of the Mehler function. Comprehensive consideration of the data size in each layer and the magnetic sensors do not easily capture the weak middle-and small-scale lithospheric magnetic field signals with lengths less than 200 km at the CHAMP orbital altitude. Therefore, K i max = 15, K e max = 10 are selected to correspond to the minimum wavelengths of magnetic fields in the horizontal direction, which are approximately 200 and 250 km, respectively. The Mehler function mainly provides a stable space continuation in the radial direction to prevent the shock effect of the model in the data gap layer (340~450 km). After testing, P max = 5 is set, and this can balance the data resolution in radials and the fitting in each height layer.

Spherical Cap Splicing
An overlapping area of approximately 7.8 • exists between adjacent spherical caps. In the overlapping area, certain differences exist in the geomagnetic components fitted by R-SCHA models in two subregions. These differences involve splicing problems between adjacent spherical caps. The present paper initially utilizes the R-SCHA model to calculate the fitting value in grid spots in overlapped areas and analyze the mutual difference. If the mutual difference is greater than the given threshold value (i.e., 0.5 nT), then the mean fitting value of each subregion is regarded as the observed value. In addition, iterative inversion is conducted to gradually obtain fitted geomagnetic components in overlapped area grid spots in two subregions that are approximately equal. Finally, its mean value is taken as the final fitting result. Figure 5 shows the splicing flow of the fitting result of adjacent spherical caps.

Error Analysis
Considering that the inversion model adopts data sets of different heights, the intensity of the lithospheric magnetic field varies significantly at different height layers. Therefore, errors are fitted at CHAMP and Swarm satellite altitudes separately to analyze the consistency between the R-SCHA model fitting value and the original observation data. The model inversion coefficients in Regions 1 and 2 are used to calculate the R-SCHA model inversion values in two regions and to seamlessly splice the overlapped region. The difference between the model inversion value and the original observation over the Tibetan Plateau and the root mean square error (RMSE) are calculated. Table 2 lists the statistical results of the errors of X, Y, and Z geomagnetic elements. Table 2. Fitting errors of X, Y, and Z components on the height layer of CHAMP and Swarm.  Table 2 shows that model fitting error decreases as altitude increases. The model fitting error in the CHAMP satellite altitude is slightly larger than that of Swarm satellite. However, the RMSE of all geomagnetic components are better than 0.3 nT. Figure 6 shows the difference values ( Z  ) between the R-SCHA and CHAMP satellite at orbital altitude of 330 ± 5 km (a) and Swarm satellite at 470 ± 5 km (b). Figure 6 shows that the maximum errors are less than 0.8 and 0.6 nT at the height layers, respectively. This comparative analysis indicates that the inversion model is consistent with the original observation. The R-SCHA algorithm can accurately describe the distribution characteristics of the lithospheric magnetic field in the Tibetan Plateau. Figure 6 shows a slight shock in the boundary regions of two independent spherical caps. The shock may be related to the large extension degree of the Mehler basis function. On the edge of the spherical cap, the value of the Mehler basis function increases constantly with extension degree. This trend affects the fitting accuracy of the edge data. The lack of boundary data and uneven distribution also strongly influence the fitting effect.

Error Analysis
Considering that the inversion model adopts data sets of different heights, the intensity of the lithospheric magnetic field varies significantly at different height layers. Therefore, errors are fitted at CHAMP and Swarm satellite altitudes separately to analyze the consistency between the R-SCHA model fitting value and the original observation data. The model inversion coefficients in Regions 1 and 2 are used to calculate the R-SCHA model inversion values in two regions and to seamlessly splice the overlapped region. The difference between the model inversion value and the original observation over the Tibetan Plateau and the root mean square error (RMSE) are calculated. Table 2 lists the statistical results of the errors of X, Y, and Z geomagnetic elements.  Table 2 shows that model fitting error decreases as altitude increases. The model fitting error in the CHAMP satellite altitude is slightly larger than that of Swarm satellite. However, the RMSE of all geomagnetic components are better than 0.3 nT. Figure 6 shows the difference values (∆Z) between the R-SCHA and CHAMP satellite at orbital altitude of 330 ± 5 km (a) and Swarm satellite at 470 ± 5 km (b). Figure 6 shows that the maximum errors are less than 0.8 and 0.6 nT at the height layers, respectively. This comparative analysis indicates that the inversion model is consistent with the original observation. The R-SCHA algorithm can accurately describe the distribution characteristics of the lithospheric magnetic field in the Tibetan Plateau. Figure 6 shows a slight shock in the boundary regions of two independent spherical caps. The shock may be related to the large extension degree of the Mehler basis function. On the edge of the spherical cap, the value of the Mehler basis function increases constantly with extension degree. This trend affects the fitting accuracy of the edge data. The lack of boundary data and uneven distribution also strongly influence the fitting effect. To further verify the stability of the model in the radial direction based on the model coefficient, the error between inversion model and EMM2015 (16°~720°) at an altitude range of 250~500 km every 25 km is computed based on a 0.5° × 0.5° grid resolution. Given that the absolute error of the model decreases as height increases, RMSE (in %) is used to show the statistical result in Figure 7b. Figure 7 does not indicate an obvious increase of errors (in %) from 350 km to 450 km, although the data lack constraints at these intermediate altitudes. The R-SCHA model can provide a stable spatial continuation in the radial direction with a reasonable setting of the truncation degree.

Geological Structures Analysis of Magnetic Anomalies in the Tibetan Plateau
The geological structures of the Tibetan Plateau are complex and have a unique and deep tectonic background. The boundary between the northern edge of plateau and the Tarim Basin is the Altun Tagh fault zone. The northeast part of the Tibetan Plateau is the Qaidam Basin. The plateau mainly includes the Songpan-Ganzi block, the Qiangtang-Chengdu block, the Lhasa block, and the Himalaya block from north to south. These blocks are bounded by the LungmuCo-Jinshajiang suture zone, the Bangong-Nujiang suture zone, and the Yarlung-Zangbo suture zone [52,53] (Figure 8).
To analyze the relationship between the magnetic anomalies and geological structures of the Tibetan Plateau, a contour map of the lithospheric magnetic field in the Tibetan Plateau at a height of 300 km based on 0.5° × 0.5° grid resolution is first drawn (Figure 9). The red region in the figure represents the positive anomaly, and the blue region denotes the negative anomaly. To further verify the stability of the model in the radial direction based on the model coefficient, the error between inversion model and EMM2015 (16 •~7 20 • ) at an altitude range of 250~500 km every 25 km is computed based on a 0.5 • × 0.5 • grid resolution. Given that the absolute error of the model decreases as height increases, RMSE (in %) is used to show the statistical result in Figure 7b. Figure 7 does not indicate an obvious increase of errors (in %) from 350 km to 450 km, although the data lack constraints at these intermediate altitudes. The R-SCHA model can provide a stable spatial continuation in the radial direction with a reasonable setting of the truncation degree. To further verify the stability of the model in the radial direction based on the model coefficient, the error between inversion model and EMM2015 (16°~720°) at an altitude range of 250~500 km every 25 km is computed based on a 0.5° × 0.5° grid resolution. Given that the absolute error of the model decreases as height increases, RMSE (in %) is used to show the statistical result in Figure 7b. Figure 7 does not indicate an obvious increase of errors (in %) from 350 km to 450 km, although the data lack constraints at these intermediate altitudes. The R-SCHA model can provide a stable spatial continuation in the radial direction with a reasonable setting of the truncation degree.

Geological Structures Analysis of Magnetic Anomalies in the Tibetan Plateau
The geological structures of the Tibetan Plateau are complex and have a unique and deep tectonic background. The boundary between the northern edge of plateau and the Tarim Basin is the Altun Tagh fault zone. The northeast part of the Tibetan Plateau is the Qaidam Basin. The plateau mainly includes the Songpan-Ganzi block, the Qiangtang-Chengdu block, the Lhasa block, and the Himalaya block from north to south. These blocks are bounded by the LungmuCo-Jinshajiang suture zone, the Bangong-Nujiang suture zone, and the Yarlung-Zangbo suture zone [52,53] (Figure 8).
To analyze the relationship between the magnetic anomalies and geological structures of the Tibetan Plateau, a contour map of the lithospheric magnetic field in the Tibetan Plateau at a height of 300 km based on 0.5° × 0.5° grid resolution is first drawn (Figure 9). The red region in the figure represents the positive anomaly, and the blue region denotes the negative anomaly.

Geological Structures Analysis of Magnetic Anomalies in the Tibetan Plateau
The geological structures of the Tibetan Plateau are complex and have a unique and deep tectonic background. The boundary between the northern edge of plateau and the Tarim Basin is the Altun Tagh fault zone. The northeast part of the Tibetan Plateau is the Qaidam Basin. The plateau mainly includes the Songpan-Ganzi block, the Qiangtang-Chengdu block, the Lhasa block, and the Himalaya block from north to south. These blocks are bounded by the LungmuCo-Jinshajiang suture zone, the Bangong-Nujiang suture zone, and the Yarlung-Zangbo suture zone [52,53] (Figure 8).  The large-background lithospheric magnetic field in the hinterland of the Tibetan Plateau shows the negative anomalies that formed a sharp contrast with positive anomalies observed in the Tarim Basin located in the north, the Sichuan Basin located in the east, and the Indian subcontinent located in the southwest. This demonstrates that an immense difference exists in the basement tectonics of the Plateau hinterland and the surrounding blocks. The midwestern part of the weak magnetic anomaly region is approximately distributed in a strip shape along the east-west direction. This distribution is basically consistent with the orientation of the regional geological structures, and its boundary basically agrees with the boundary of regional tectonic. However, in the east of the plateau, the magnetic anomaly gradually extends along the northeast direction, thereby exhibiting an arc shape. Thus, the magnetic anomaly is inconsistent with the east-south distribution of the geological structures zone. This phenomenon indicates that the anomaly in the east of the plateau may be caused by the magnetic materials in the middle and shallow layers of the earth's crust.
The satellite magnetic anomaly map shows that the magnetic anomaly intensity and its distribution in all secondary structures zones in the Tibetan Plateau from the north to the south differ. A regional clump by positive anomaly is distributed in the middle part of the Tarim Basin. The western region adjacent to the Arkin fault zone shows a linear positive anomaly, but the eastern regions adjacent to the Qilian Mountain fault zone exhibit a negative anomaly. However, the positive and negative anomalies are weak. The Songpan-Ganzi block and the Qiangtang-Chengdu To analyze the relationship between the magnetic anomalies and geological structures of the Tibetan Plateau, a contour map of the lithospheric magnetic field in the Tibetan Plateau at a height of 300 km based on 0.5 • × 0.5 • grid resolution is first drawn (Figure 9). The red region in the figure represents the positive anomaly, and the blue region denotes the negative anomaly.  The large-background lithospheric magnetic field in the hinterland of the Tibetan Plateau shows the negative anomalies that formed a sharp contrast with positive anomalies observed in the Tarim Basin located in the north, the Sichuan Basin located in the east, and the Indian subcontinent located in the southwest. This demonstrates that an immense difference exists in the basement tectonics of the Plateau hinterland and the surrounding blocks. The midwestern part of the weak magnetic anomaly region is approximately distributed in a strip shape along the east-west direction. This distribution is basically consistent with the orientation of the regional geological structures, and its boundary basically agrees with the boundary of regional tectonic. However, in the east of the plateau, the magnetic anomaly gradually extends along the northeast direction, thereby exhibiting an arc shape. Thus, the magnetic anomaly is inconsistent with the east-south distribution of the geological structures zone. This phenomenon indicates that the anomaly in the east of the plateau may be caused by the magnetic materials in the middle and shallow layers of the earth's crust.
The satellite magnetic anomaly map shows that the magnetic anomaly intensity and its distribution in all secondary structures zones in the Tibetan Plateau from the north to the south differ. A regional clump by positive anomaly is distributed in the middle part of the Tarim Basin. The western region adjacent to the Arkin fault zone shows a linear positive anomaly, but the eastern regions adjacent to the Qilian Mountain fault zone exhibit a negative anomaly. However, the positive and negative anomalies are weak. The Songpan-Ganzi block and the Qiangtang-Chengdu The large-background lithospheric magnetic field in the hinterland of the Tibetan Plateau shows the negative anomalies that formed a sharp contrast with positive anomalies observed in the Tarim Basin located in the north, the Sichuan Basin located in the east, and the Indian subcontinent located in the southwest. This demonstrates that an immense difference exists in the basement tectonics of the Plateau hinterland and the surrounding blocks. The midwestern part of the weak magnetic anomaly region is approximately distributed in a strip shape along the east-west direction. This distribution is basically consistent with the orientation of the regional geological structures, and its boundary basically agrees with the boundary of regional tectonic. However, in the east of the plateau, the magnetic anomaly gradually extends along the northeast direction, thereby exhibiting an arc shape. Thus, the magnetic anomaly is inconsistent with the east-south distribution of the geological structures zone. This phenomenon indicates that the anomaly in the east of the plateau may be caused by the magnetic materials in the middle and shallow layers of the earth's crust.
The satellite magnetic anomaly map shows that the magnetic anomaly intensity and its distribution in all secondary structures zones in the Tibetan Plateau from the north to the south differ. A regional clump by positive anomaly is distributed in the middle part of the Tarim Basin. The western region adjacent to the Arkin fault zone shows a linear positive anomaly, but the eastern regions adjacent to the Qilian Mountain fault zone exhibit a negative anomaly. However, the positive and negative anomalies are weak. The Songpan-Ganzi block and the Qiangtang-Chengdu block are mainly manifested in a stable negative anomaly that approximates an east-west distribution. The eastern regions near the thrust nappe belt of the Longmen Mountains shows a northeast-southwest negative anomaly. Weak negative anomalies of these blocks imply that the Curie isothermal surface in these blocks is deeper than that in the south of the plateau. Magnetic anomalies in the Lhasa (Gangdise belt) and Himalaya blocks approximately show an east-west orientation along the west section of the Himalaya and across the Himalaya mountains, and this orientation continuously extends to the east region of Lhasa with a distinct negative anomaly. The strongest negative anomaly is located in the Himalaya region. The center of this anomaly is located approximately at 83.5 • E and 29 • N and has a center strength that reaches −9 nT. The negative anomaly may be a result of the strong collision between the blocks. Thrust nappes on the boundary of the block generate a friction phenomenon, thereby promoting constant increases in the temperature of deep crustal materials and the upward migration of hot materials. These effects lead to the Curie isothermal surface uplift. This explanation can be proved by the high heat flow activities in the Tibetan Plateau [54]. Therefore, a strong negative anomaly may imply that deep geological structures in this region are highly active.

Conclusions
Zonal revised of spherical cap harmonic analysis(R-SCHA) was used to model the lithospheric magnetic field of the Tibetan Plateau at satellite altitude. The accuracy and stability of this model were tested and evaluated in the horizontal and radial directions. Furthermore, the relationship between magnetic anomalies and geological structures was analyzed.
One of the key factors for modeling is refining satellite magnetic data. A complete data preprocessing procedure was designed, which involved data screening, model correction, track filtering, and hierarchical gridding. Reliability is another important factor of a model. Thus, the model must not generate significant shock effects in the horizontal and vertical directions but still achieve a certain spatial resolution. Considering that the difference between the latitude and longitude over the Tibetan Plateau is large, a small cap partition modeling can be used in a particular truncation level to provide good spatial resolution while ensuring accuracy. Therefore, the Tibetan Plateau was divided into two independent parts and the half-angle of the cap (approximately 1 • ) was moderately amplified to eliminate the influence of the boundary effect. Meanwhile, hierarchical gridding was conducted for satellite data distributed at a height of 250~510 km to guarantee the balance of each layer of data. The small semi angle of cap improved the fitting accuracy of the model throughout the region, but the difficulty of this method lies at the joint of each subregion model. The present paper used iterative inversion to create a method for seamless splicing between adjacent spherical caps. Moreover, model error and reliability analyses played important roles in the modeling process. The difference between the model value and the satellite observation value were initially calculated before the RMS of the model was quantitatively analyzed at the heights of the CHAMP and Swarm satellites. The fitting error of the model is small, and the RMS of each geomagnetic component is better than 0.3 nT. Moreover, the stability of the model in horizontal and vertical directions was examined by introducing the EMM2015 (16~720 order) high-accuracy lithospheric magnetic field model. The model does not have any significant shock effects in the horizontal direction. The model can provide a stable space continuation in the radial direction. The analysis indicates that the R-SCHA algorithm can accurately describe the middle-to long-scale lithospheric magnetic field in the Tibetan Plateau at the satellite altitude layer.
The separation of internal and external field sources remains a problem in the modeling process. The R-SCHA algorithm does not consider the joint estimation of the external source field. Instead, it facilitates the separation of an external field source through the model and filtering method in the data preprocessing steps. This method may introduce new errors in the process of model fitting. Therefore, Combining Swarm constellation differential gradient measurement is necessary to effectively suppress the interference of the long wavelength signal that remains in the external field and the main magnetic field. However, a map of a high-resolution lithospheric magnetic field over the Tibetan Plateau cannot be developed because of the lack of support for ground magnetic survey and aeromagnetic data. Although a good relationship between the distribution of magnetic anomalies and deep crustal tectonic structures at 300 km satellite altitude was identified, the model fails to extract abundant shallow geological structures information. Future research will be conducted on a downward continuation algorithm for geomagnetic anomaly data and multi-scale modeling will be integrated to obtain a high-resolution and high-accuracy lithospheric magnetic field. This endeavor will provide technical support for the study on the crustal structures of the Tibetan Plateau and related dynamic mechanisms.