Study of the Correction Method of Secular Variation in the Main Geomagnetic Field by the Field Seismogeomagnetic Survey

: The correction of the secular variation (SV) of the main geomagnetic field is a key link of field seismogeomagnetic data processing, and the current method relies on the observatory data for the relevant technical processing. To optimize the data products and obtain more accurate and reliable seismomagnetic information, this study adopted a new technical idea, which uses the repeated survey data from field stations to obtain the SV of the main geomagnetic field over the survey area by the weighted least - squares method, and compared the results with those of the current methods. The results were as follows: 1. The SV results of the main geomagnetic field produced by the new method are closer to those of the International Geomagnetic Reference Field (IGRF)_SV model. The mean square error (MSE) of the difference of the three elements F, D, and I between the new method and the IGRF_SV model is 10.7%, 47.0%, and 14.5% of that of the original method, respectively. 2. By applying the new SV correction method, more stable and reasonable variations in Earth’s crustal magnetic field can be obtained. The average amplitude of the Earth’s crustal magnetic field variation in the three elements F, D, and I is 28.5%, 55.4%, and 34.4 of the original results, the MSE is 59.1%, 56.5%, and 40.3% of the original results, and the mean gradient is 93.6%, 91.9%, and 97.0%, respectively. 3. In the processed results of the new method, the seismomagnetic information is clearly optimized, and the location of the epicenter is more consistent with the 0 value line of the Earth’s crustal magnetic field. The processed results of the new method are significantly better than those of the original method and have a higher application value.


Introduction
The field seismogeomagnetic survey is a repeated measurement taken periodically at field sites using high-precision geomagnetic measuring instruments to monitor anomalous changes in the geomagnetic field and to serve for earthquake prediction.In mainland China, this survey has a history of more than 30 years.The China Earthquake Administration started to implement large-scale geomagnetic field vector surveys in 2009 and gradually built a network of field stations that basically covers the whole of mainland China, with a total of more than 1300 field stations, which plays an important role in earthquake prediction (Gu,2017).In the presence of abundant observation data, how to obtain effective seismomagnetic information through scientific and reasonable technical ways is the focus of data processing.
The geomagnetic field consists of the main magnetic field, the Earth's crustal magnetic field (or lithospheric magnetic field), and the varying magnetic field.The purpose of the field seismogeomagnetic survey is to obtain the variation information of the Earth's crustal magnetic field to extract the seismomagnetic information, which requires the removal of the varying magnetic field and secular variation (SV) of the main geomagnetic field from the measured data, resulting in two main data processing steps of diurnal variation reduction and SV of the main geomagnetic field correction (Chen et al.,2017).The Earth's crustal magnetic field is much less variable than the other two fields, and the seismomagnetic information is even more complex and variable (Ding et al.,2011), which puts forward very demanding technical requirements for data processing.At present, some research results have been obtained on diurnal variation reduction (Roden R B et al.,1964;Sander E L,1982;Whitmarsh R B & Jones M T,2010;Pan et al.,2020).Studies of SV of the main geomagnetic field (hereafter referred to as SV) mainly focused on the establishment methods of global and regional models and the characterization of spatial and temporal variations (Alldredge L R,1981;Malin S R C,1985;Korte M & Haak V;Heikki NEVANLINNA,2010;Finlay C C et al.,2016;Rother M et al.,2021).Field seismogeomagnetic surveying requires high accuracy of the data products, and more accurate, reliable, and targeted SV correction methods are needed.
The currently used SV correction method is as follows: after performing the natural orthogonal decomposition (NOC) on the data from 31 basic geomagnetic observatories in China, the top six principal components of each observatory are spatially interpolated to calculate the SV values at the field stations to complete the SV correction (Chen,2011).The method relies on the daily observation data from each observatory, and the processed results may be problematic in some regions or periods.This targeted study adopted a new technical idea to construct a new SV model based on repeated survey data from high-density field stations, and it quantitatively compared and evaluated the processed results of the two methods using historical data.

Data Information
The North-South Seismic Belt (NSSB) in mainland China is a seismically active area and a key monitoring area for field seismogeomagnetic measurements where systematic geomagnetic vector observations have been conducted since 2009.The locations of the field stations are buried in the ground with artificial markers, which are fixed during each period of operation.The operation is conducted for one period each year, which is from March to May.Six sets of data are collected at each field station, and the elements are F, D, and I.The results of the new and old SV correction methods were compared using the two-period data from 186 field stations in the area in 2016 and 2017.
The purpose of seismogeomagnetic monitoring is to obtain seismomagnetic information.From July 2017 to December 2018, the seismic activity within the survey area was very active, and a total of eight earthquakes of magnitude 5.0 or greater occurred, including the Jiuzhaigou earthquake of magnitude 7.0 on August 8, 2017 [Table 1].Using these earthquakes, the effects of different SV correction methods on the seismomagnetic information were examined, and the application value of the new method was tested.The locations of the relevant observatories, field stations, and earthquakes are shown in Figure 1.

Weighted Least-Squares Method
After the field geomagnetic survey data are generalized by diurnal variation reduction, the observation data from each field station are uniformly corrected to the same time.The difference in the diurnal variation reduction data in the survey area is obtained by subtracting the diurnal variation reduction data from the two periods, which contains the components of SV and the Earth's crustal magnetic field variation.These two components have different temporal and spatial variation characteristics, and SV is the main component, with a smooth and regular spatial distribution pattern (Xu,2009).An appropriate modeling method can be used to extract the stable main data component in the above difference, which can be used as the SV of the survey area during the two periods.
The Taylor polynomial is an important geomagnetic field modeling method that can effectively describe the spatial distribution pattern of regional SV (An et al.,1982).It uses the least-squares method to fit the data, i.e., the sum of squares of the residual of each data point is minimized.During the fitting, the weight of each data point is the same.The survey results may be disturbed by factors, such as measurement errors, environmental changes, and even seismomagnetic information, resulting in anomalies in the data from some field stations.These anomalies usually have large residuals, which clearly affect the fitting, cause errors in SV results, and directly interfere with the subsequent extraction of the seismomagnetic information.
To suppress the influence of these anomalies, the weighted least-squares method (Schaffrin B & Wieser A,2008;Vahid, Mahboub,2011) was used for the fitting.In other words, the conventional least-squares method was first used to calculate the residual of each field station, and then a different weight was assigned to each field station according to the principle of larger residuals, smaller weights.Thus, the weighted least-squares method was used for fitting.An iterative algorithm was used to repeat the fitting until the fitting results were stabilized, and then the final fitting coefficients were obtained.The main technical aspects are expressed as follows: The weighted least-squares method is implemented by an iterative algorithm.The conventional least-squares method is used first.The spatial distribution of SV in the region is simple and regular, so the order of polynomials is two.Let there be m points in the survey area.For the i th point (i=1, 2, ..., m), the longitude is   , the latitude is   , and the measured value is   .The vector Z of dependent variables, the matrix A of independent variables, the vector  of coefficients, and the residual vector  are as follows: (1) Then, the matrix form of the linear equations to be solved is Using the conventional least-squares method, the coefficient vector is The weighted least-squares method is based on the conventional least-squares method, and according to the residual  of the previous fit, different weights are assigned to field stations.Let the weight of the i th point be   (i=1,2,…,m), and the weight vector W is Then, when using the weighted least-squares method, the vector of coefficients is

Weight Function
Before applying the weighted least-squares method, a weight function needs to be determined.The weight function uses the residual  of the previous fit as the independent variable, and the dimensions and changes of elements F, D, and I are different.To unify the calculation results, the  of the three elements should first be standardized(Van Eck N J P & Waltman L,2014).Let the number of measurement points be m, the set of residuals of the previous fit of a certain element be { 1 ,  2 , ⋯ ,   , ⋯ ,   }, the mean value be ̅ , and the standard deviation (SD) be   .Then, the standardized residual   ′ of the i th point is The normal distribution is a probability distribution that exists widely in various scientific fields.The probability density function of the standard normal distribution was selected (Martarelli G C & Zanini A,1975) as the weight function; then, the weight   of the point corresponding for   ′ is The weight function (Figure 2) shows that the weights of the measurement points converge rapidly with increasing residuals.When the absolute values   ′ are set as 1, 2, and 3 (indicating that the deviation of the residuals from the mean value of the survey area is 1-, 2-, and 3-fold the SD, respectively), the corresponding weights are 60.7%, 13.5%, and 1.1% of the baseline weight (the weight at   ′ =0).The effect of anomalies is greatly suppressed.The weighted least-squares method was used to perform an SV correction on the 2016 and 2017 difference in diurnal variation reduction data in the NSSB, and the results of Earth's crustal magnetic field variation in the survey area were recalculated.When other technical aspects were identical, the results were compared with the original processed results in terms of the morphological characteristics of Earth's crustal magnetic field changes and the display of seismomagnetic information (Figure 3).The comparison of the new and old maps shows that the spatial distribution patterns of the two results for Earth's crustal magnetic field change have obvious differences.Compared with the original results, the distributions of positive and negative anomalous regions become balanced and scattered in the new method results.

Morphological Characteristics
To quantitatively describe the morphological characteristics of the map, a statistical analysis of the mean absolute value, MSE, and mean gradient of the data from two maps was conducted [Table 2], which represent the magnitude, dispersion, and uniformity of the spatial distribution of the Earth's crustal magnetic field changes, respectively.According to [Table 2], the mean absolute value of the results of the new method for the three elements F, D, and I is 28.5%, 55.4%, and 34.4 of the original results; the MSE is 59.1%, 56.5%, and 40.3% of the original results; and the mean gradient is 93.6%, 91.9%, and 97.0% of the original results, respectively.Compared with the original results, the average amplitude and dispersion of the new method are significantly reduced, and the spatial distribution of uniformity is slightly increased, indicating that the magnetic field variation in the survey area obtained by the new SV correction method has a smaller amplitude and is more stable and uniform.The Earth's crustal magnetic field originates from rock magnetism, and its spatial structure is complex and stable in time (Xu et al.,2008) , which means that the results of the new method are closer to the theoretical pattern of Earth's crustal magnetic field variation and more reasonable.

Seismomagnetic Information
The purpose of seismogeomagnetic monitoring is to capture seismomagnetic information for the prediction of potential epicenter locations.Can the new SV correction method help to capture the seismomagnetic information?The existing experience shows that the 0 value line and the high-gradient region of the Earth's crustal magnetic field variation are the seismically active areas (Gu et al.,2010;Feng,2019).Based on this empirical understanding, the contour lines near the epicenter were analyzed.
The Jiuzhaigou earthquake of magnitude 7.0 on August 8, 2017, (the largest red dot in Figure 3) was taken as an example.At the epicenter, compared with the original results, the density of contour lines of the three elements F, D, and I of the new results is not changed significantly.However, the agreement between the epicenter and the 0 value of the three elements is greatly improved, and the location characteristics become clear.The other seven earthquakes of magnitude 5 or greater (only six small dots are shown in Figure 3 because of the overlapping of two earthquakes) also show similar results.The above results show that the new SV correction method can effectively strengthen the display of seismomagnetic information and facilitate the extraction of seismomagnetic information.

Traceability and Evaluation
The two results obtained for the Earth's crustal magnetic field variation based on the new and old SV correction methods differ greatly, indicating that the two methods yield significantly different SV results, and the rationality must be evaluated.The results of the International Geomagnetic Reference Field (IGRF)_SV model for the survey area were calculated as a benchmark for comparative evaluation.Meanwhile, the inconsistency of SV results obtained based on data from different sources in the same area indicates that there may be problems with one data source.The daily relative observations produced by the observatories need to be corrected for monthly and annual baseline values before they become the final data products.For time-sensitive reasons, the final data cannot be used when correcting the SV for field survey data, and the data used were only provisionally corrected for baseline values.That is, these data were not fully quality controlled and may affect the SV results.To trace the problems, the SV results of the survey area were recalculated using the final data (data corrected by monthly and annual baselines) of the relevant observatories for comparison.The SV of the survey area obtained by the four methods is shown in Figure 4.
The IGRF_SV results (denoted as SV_IGRF) of the survey area have a simple and smooth shape.In contrast, the SV results of the original method (denoted as SV_old) have clearly different morphological characteristics, i.e., the contour lines of F and D elements show severe distortion; the contour lines of element I are basically the same as those of SV_IGRF, but they are clearly too dense, and the spatial variation rate is too large.The SV results using the final observatory data (denoted as SV_observatory) are significantly more consistent with SV_IGRF, but there are still some local distortions, and the location and shape of the distortions are not consistent with SV_old.The results of related studies have shown regional features in SV (Xu et al.2005;Peqini K et al.,2018).However, in the two SV results using observatory data, the positions of distortions are inconsistent, and the morphology is significantly different.These distortions can't be found in the other two SV results, which indicates that these local distortions are not real and reliable and are due to the data distortions of some observatories.In contrast, the SV results of the new method (denoted as SV_new) are most similar to SV_IGRF.To quantify the agreement between the three SV results and SV_IGRF, the correlation coefficient and the MSE of the difference between SV_IGRF and the other three SV results were calculated.The correlation coefficients show that SV_observatory has the highest linear correlation with SV_IGRF, followed by SV_new and then SV_old; the MSE shows that SV_new has the highest agreement with SV_IGRF, and it is significantly better than SV_observatory and SV_old.In general, SV_new is closest to SV_IGR.The above statistical results show that the SV results obtained using observatory data may have distortion problems in some periods and regions.The use of observatory data cannot avoid the problems of instrument stability and environmental disturbances (Li et al.,2019), which have a hard-to-quantify impact on the accuracy of the observatory data.This problem cannot be completely avoided even when using the final data products from the observatories.Field seismogeomagnetic surveying requires extremely high accuracy of processed results, and failure to adequately check the quality of observatory data before performing the SV correction may have severe consequences.IGRF_SV is a global data model for the SV correction of field seismogeomagnetic surveys, but it is too simple and has insufficient accuracy.In addition, the difference between the mean value of IGRF_SV and those of the other three SVs was calculated [Table 4], and the mean values of all three elements of IGRF_SV are higher than those of the other three SVs.There is a systematic bias that cannot be ignored.Although the technical conditions of field surveys are inferior to those of geomagnetic observatories, the number and spatial density of field stations are far better than those of geomagnetic observatories, and the repeated survey data from many field stations can be used to obtain more ideal SV results in the survey area by the weighted least-squares method.Verhoef J and Williams C A used the repeated survey data from the intersection points of marine magnetic survey lines to extract SV information (Verhoef J & Williams C A,1993), and better results were obtained.The accuracy of field seismogeomagnetic surveys and the density of field stations are much better than those of marine magnetic surveys, so the accuracy of processed results should be better.This method extracts the SV components from the field survey data, which is adaptive and avoids the systematic bias that may be brought by using data from other sources.Compared with the existing SV correction method, this method is simple, fast, and reliable, optimizes the seismomagnetic information, and has a strong practical value.
Earthquake prediction is a worldwide scientific challenge, and seismogeomagnetic monitoring also faces many problems.Only by using practical work as a basis and constantly seeking technological progress in the process of finding and solving problems can we gradually approach our goal.Declarations

Figure
Figure 2 Weight function curve 2 Calculation Examples

Figure 3
Figure 3 Earth's crustal magnetic field variation in the survey area (The left panel is the results of the original method, and the right panel is the results of the new method; row 1 of the panels corresponds to one geomagnetic element, from top to bottom, F, D, and I)

Figure 4
Figure 4 SV results of the survey area by the four methods

Table 2
Table of the results of Earth's crustal magnetic field variation

Table 3
Comparison results between IGRF_SV and the other three SVs in terms of agreement

Table 4
Systematic difference between SV_IGRF and the other three SV results