An hourly inversion method for ionospheric equivalent current system using CM model

Abstract

The invention discloses an hourly inversion method of an ionospheric equivalent current system by introducing a CM model. The method collects and records the original data of the geomagnetic field vector of global geomagnetic stations, eliminates the background changes such as the long-term term and the trend term in the original data, and then performs principal component analysis on the single station data to obtain the geomagnetic diurnal variation field component as the geomagnetic diurnal variation field basic data for hourly inversion of the ionospheric equivalent current; and then uses the CM magnetic field model to calculate the model magnetic field data of the ionospheric current magnetic field component of the global grid point. The measured data and the model data are further combined for gridding interpolation to obtain the global geomagnetic diurnal variation field data for hourly inversion of the ionospheric equivalent current; finally, based on the spherical harmonic analysis method, the global geomagnetic diurnal variation field gridding data is separated into internal and external source fields to obtain the external source field corresponding to the ionospheric equivalent current magnetic potential, and then according to the relationship between the magnetic field and the current given by the Biot-Saffar law, the hourly ionospheric equivalent current is obtained by inversion calculation.

Description

Technical Field

[0001] The invention belongs to the field of digital combing and research of earth's changing magnetic field and ionospheric current data, and relates to an hourly inversion method of an ionospheric equivalent current system by introducing a CM model. Background Art

[0002] The inversion algorithm of the ionospheric equivalent current system is based on the geomagnetic field potential field theory proposed by mathematician Gauss in 1839. With the help of the Gauss-Schmidt method, the basic magnetic field of the earth can be divided into two parts: internal and external sources. The basic idea of ​​the algorithm is: assuming that the current that generates the magnetic field is located on a two-dimensional spherical shell at a certain height (such as the E region of the ionosphere), first use the forward method (spherical harmonic analysis theory) to establish the magnetic field expression generated by the known current system, and then solve the Gaussian coefficient based on the global geomagnetic field observation data, and then use the relationship between the magnetic field and the current to derive the current. The external source field corresponds to the ionospheric equivalent current system. Based on the potential field theory, the spherical harmonic analysis method realizes a concise mathematical expression of the spatial distribution of geophysical parameters, and therefore has been widely used in fields such as geomagnetism and gravity.

[0003] The geomagnetic diurnal variation field (abbreviated as Sq) is the data basis of the ionospheric equivalent current inversion algorithm, so its data quality and the number and distribution of stations are important factors affecting the inversion results. Since the geomagnetic field record is a superposition of multi-source magnetic field signals inside and outside the earth, how to effectively extract the geomagnetic diurnal variation field is a key link in the inversion calculation. After eliminating the long-term changes in the geomagnetic field record, the average calculation can effectively eliminate short-term magnetic disturbances. Therefore, the most commonly used ionospheric equivalent current inversion methods are monthly inversion and daily inversion, that is, the daily geomagnetic diurnal variation field is extracted by monthly averaging or Fourier transform, and its data is used as the data input of the inversion algorithm, thereby realizing the monthly / daily inversion of the ionospheric equivalent current system.

[0004] However, for the inversion of the ionospheric equivalent current system with higher time precision, such as hourly inversion, it is necessary to extract the geomagnetic diurnal variation field of global stations at the same time, and then re-grid interpolate the "instantaneous" value data of the diurnal variation field of all global stations at a single time, and then calculate the global "diurnal variation field" data. Obviously, due to the limitations of topography and economic conditions of various countries, the spatial distribution of global geomagnetic stations has obvious defects, with dense stations in Europe and the United States and no stations in the ocean area. The uneven spatial distribution of global geomagnetic stations has greatly limited the time precision of the ionospheric equivalent current inversion, which seriously restricts further understanding of the instantaneous spatial distribution and dynamic changes of the ionospheric equivalent current.

[0005] The principal component analysis method has the characteristics of "natural generation" from the data without artificially specifying the basis function. In recent years, it has been widely used in the analysis of ionosphere and geomagnetic data, and can improve the extraction quality of geomagnetic diurnal variation field. In addition, overcoming the defect of uneven distribution of global stations and obtaining relatively ideal global grid data of geomagnetic diurnal variation field is the key to achieving the stability of hourly inversion of ionospheric equivalent current. The CM (Comprehensive Modeling) comprehensive magnetic field model, as a new geomagnetic field modeling method jointly developed by NASA and the Danish Space Research Institute in the 1990s, has the ability to calculate the global basic magnetic field and changing magnetic field. Among them, the ionospheric field and its induced field module unit can calculate the geomagnetic diurnal variation field generated by ionospheric current and its induced current. Therefore, the introduction of CM model values ​​as a supplement to ground-based station data can effectively make up for the insufficiency of global station spatial distribution to meet the data needs of hourly inversion of ionospheric equivalent current. Summary of the invention

[0006] In order to solve the above technical problems, the purpose of the present invention is to provide an hourly inversion method for the ionospheric equivalent current system by introducing the CM model. First, the principal component analysis method is used to extract the single station geomagnetic diurnal variation field, and then the model data is calculated according to the CM geomagnetic field model as the site data supplement, and the measured data and the CM model data are combined for grid processing and calculation to obtain the "instantaneous (hourly)" global geomagnetic diurnal variation field data, and then the spherical harmonic analysis algorithm is used to realize the hourly inversion of the ionospheric equivalent current.

[0007] The hourly inversion method of an ionospheric equivalent current system introducing a CM model described in the present invention comprises the following steps:

[0008] S1 collects and records the original data of geomagnetic field vectors from global geomagnetic stations;

[0009] S2 eliminates background changes such as long-term items and trend items in the original data. Specifically, the daily data series is first converted from the world time coordinates to the local time coordinates according to the conversion relationship. Then, the long-term items are calculated by the average of the data from 23:00 to 01:00 local time on the same day of the station and deducted from the data. The gradient values ​​of the data at 00:00 and 23:00 local time are further used to eliminate the trend items of the data at each local time. Then the daily data series is converted from the local time coordinates to the world time coordinate system;

[0010] S3 processes the single station data to obtain the geomagnetic diurnal variation field component, which is used as the basic data of the geomagnetic diurnal variation field for hourly inversion of the ionospheric equivalent current;

[0011] S4 uses the CM magnetic field model, inputs the operation parameters, and calculates the model magnetic field data of the ionospheric current magnetic field components at the global grid points; combines the measured data and the model data for grid interpolation to obtain the global geomagnetic diurnal field grid data of the ionospheric equivalent current hourly inversion;

[0012] S5, based on the spherical harmonic analysis method, separates the internal and external source fields from the global geomagnetic diurnal field grid data, and obtains the external source field corresponding to the ionospheric equivalent current magnetic potential;

[0013] S6 uses the relationship between magnetic field and current given by the Biot-Saffar law to invert and calculate the hourly ionospheric equivalent current.

[0014] Furthermore, the conversion relationship in step S2 is local time = universal time + geographic longitude / 15.

[0015] Furthermore, the method for processing single data in step S3 is a principal component analysis method.

[0016] Furthermore, the CM magnetic field model operation parameters in step S4 are: time, longitude and latitude coordinates and geomagnetic activity index Dst.

[0017] Furthermore, the analysis method in step S5 specifically includes the following steps:

[0018] S51. Calculate the grid data of the geomagnetic daily variation field, and finally obtain the network points evenly distributed in longitude and latitude on a global scale;

[0019] S52. Solve the Laplace equation for the geomagnetic field components in the spherical coordinate system to obtain the harmonic superposition of a series of sine and cosine functions containing associated Legendre coefficients. Establish an equivalent relationship between the observed data Fourier coefficients and the spherical harmonic expansion harmonic coefficients, that is, HA = B;

[0020] S53. Establish the expression of Gaussian coefficient through the relationship between the Fourier harmonic component sine and cosine coefficients and the associated Legendre coefficients, and solve for the Gaussian coefficient;

[0021] Furthermore, the data calculation method in step S51 is:

[0022] S511 expands the grid data of the geomagnetic daily variation field by Fourier series according to the longitude direction to obtain the Fourier coefficients corresponding to each latitude;

[0023] S522. All Fourier coefficients are interpolated in the latitude direction to make them continuously smooth within a certain latitude range, thereby obtaining network points that are evenly distributed in longitude and latitude on a global scale.

[0024] Furthermore, in the equivalent relationship HA=B in step S52, H is a matrix composed of basis functions, and A is The associated Legendre coefficient matrix is ​​composed of the triangular coefficient matrix, and B is the matrix composed of triangular coefficients. It can be obtained that A = H\B, that is, the associated Legendre coefficient is obtained

[0025] Furthermore, the method for solving the Gaussian coefficient in step S53 is: least square method.

[0026] Furthermore, the interpolation method in step S522 is specifically: interpolating in the latitude direction using a cubic spline function.

[0027] Furthermore, in step S522, the latitude range is -90° to 90°.

[0028] Compared with the prior art, the present invention has the following beneficial effects:

[0029] The present invention introduces the CM model as supplementary data for areas lacking ground-based station observations, and then performs spherical harmonic analysis on the global station grid data at an instant, which can obtain the ionospheric equivalent current inversion effect with high time resolution, and the inversion has good stability and high accuracy. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 is a flow chart of the method of the present invention;

[0031] Figure 2 This is the flow chart of the spherical harmonic analysis algorithm;

[0032] Figure 3 Schematic diagram of extracting the geomagnetic diurnal variation field Sq by principal component analysis in the embodiment;

[0033] Figure 4 Schematic diagram of the correlation between the number of principal component analysis data samples and Sq in the embodiment;

[0034] Figure 5 Schematic diagram of the correlation between the quality of principal component analysis data samples and Sq in the embodiment;

[0035] Figure 6 Schematic diagram of the relationship between the number and quality of principal component analysis data samples and the first-order principal component in the embodiment Figure 7 Schematic diagram of the spatial distribution of global measured stations and CM model data grid points in the embodiment;

[0036] Figure 8 Schematic diagram of comparison of hourly inversion results of ionospheric equivalent current of measured data or / and CM model data in the embodiment;

[0037] Fig. 9 This is a schematic diagram comparing the spatial distribution and changes of the ionospheric equivalent current hour by hour before and after the introduction of the CM model data in the embodiment; DETAILED DESCRIPTION

[0038] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0039] like Figure 1 As shown:

[0040] First, the original data x(t), y(t), and z(t) of the geomagnetic field vectors of global geomagnetic stations are collected and recorded.

[0041] Next, the long-term changes and trend term changes of the single raw data x(t), y(t), and z(t) are eliminated to obtain dx(t), dy(t), and dz(t), where t is the hourly sampling rate.

[0042] Perform principal component analysis on the original dx(t), dy(t), and dz(t) of each station, that is, expand each time series into the sum of multiple components: Each time series is a sample, the sample length is n, m samples participate in the covariance matrix and introduce the orthogonal condition, and the least squares method is used to calculate the eigenvector Φ of each order component k (t), amplitude A k (s) and principal component F k (s, t). The first-order principal component obtained by calculating the daily data of the geomagnetic field is the single-station geomagnetic daily variation field ddx(t), ddy(t), ddz(t);

[0043] Based on the CM model calculation software and model coefficients developed and provided by Sabaka et al. of the Danish Space Research Institute and NASA's Goddard Space Flight Center, the CM magnetic field model is called with time and longitude and latitude coordinates as input parameters to calculate the value of the geomagnetic diurnal variation field at any point on the ground. Combined with the spatial distribution of global ground-based stations, a reasonable grid density is selected, and the CM model is called to calculate the model data of the geomagnetic diurnal variation field at the grid point; then the measured data and model data of global ground-based stations are combined, and the data are gridded and interpolated using the biharmonic spline interpolation method, thereby obtaining the global gridded geomagnetic diurnal variation field data X(θ,φ), Y(θ,φ), and Z(θ,φ);

[0044] The spherical harmonic analysis algorithm is used to separate the internal and external source fields, where the external source field corresponds to the ionospheric equivalent current system. The gridded data X(θ, φ), Y(θ, φ), and Z(θ, φ) of the global geomagnetic diurnal variation field at a certain time h, where θ is the geographic longitude and φ is the geographic latitude, are used as the algorithm input data. The specific steps are as follows Figure 2As shown:

[0045] First, the geomagnetic diurnal field grid data X(θ,φ), Y(θ,φ), and Z(θ,φ) are expanded into a series of harmonic superpositions of sine and cosine functions by Fourier series according to longitude, as follows: The Fourier coefficients Xa(m), Xb(m), Ya(m), Yb(m), Za(m), and Zb(m) are obtained for each latitude, where a represents the cosine term and b represents the sin term.

[0046] Second, all Fourier coefficients are interpolated in the latitude direction using the cubic spline interpolation method, thereby calculating the network point coefficient values ​​that are evenly distributed in the latitude range of -90° to 90°.

[0047] Third, the three components of the geomagnetic field X, Y, and Z are solved in the spherical coordinate system to obtain the harmonic superposition of a series of sine and cosine functions containing associated Legendre coefficients, and its expanded form is: and Among them, θ is the geographic co-latitude and λ is the geographic longitude. is an nth order and mth degree associated Legendre function. and Defined as the basis functions of X, Y, and Z, the basis functions can be obtained by solving the Legendre polynomials. The equivalent relationship is established by the Fourier coefficients of the observed data and the harmonic coefficients of the spherical harmonic expansion, that is, HA = B. H is the matrix composed of the basis functions, and A is The associated Legendre coefficient matrix is ​​composed of the triangular coefficient matrix, and B is the matrix composed of triangular coefficients. It can be obtained that A = H\B, that is, the associated Legendre coefficient is obtained

[0048] Fourth, according to Gaussian theory, the relationship between the Gaussian coefficient of the internal and external source fields and the associated Legendre coefficient is: The Gaussian coefficient of the internal source field is obtained by solving and the Gaussian coefficient of the external field

[0049] Finally, the Gaussian coefficient is substituted into the expression of the external magnetic potential The magnetic potential U is obtained. In the formula, R is the length of the earth's radius and r is the distance from the center of the earth. The inclusion of the (r / R)n factor indicates that the magnetic field source is outside the earth, that is, the external source field. According to the relationship between magnetic potential and equivalent current Where a is the radius of the earth, r is the distance from the center of the earth, and μ0 is the dielectric constant, and the equivalent current J can be calculated.

[0050] In the above embodiment, the specific data distribution is as follows: Figure 3-Figure 9 As shown, Figure 3The schematic diagram of the extraction of the geomagnetic diurnal variation field Sq during the magnetic storm using principal component analysis. In which, a separate operation is performed for each station data, the sample length n = 24, and the number of samples m = 8. As can be seen from the figure, its shape and amplitude are consistent with the variation characteristics of the average Sq, indicating that the principal component analysis method can effectively extract the geomagnetic diurnal variation field from the geomagnetic record.

[0051] Figure 4-6 The schematic diagram of the sample test of extracting the geomagnetic diurnal variation field Sq by principal component analysis of the original data in the embodiment. The test data is the hourly data of the horizontal component of the geomagnetic field of the Ming Tombs geomagnetic observatory in Beijing in 2008. In the principal component analysis, the sample length n=24, the sample number is selected with a sliding window of 60 days, and a total of 306 (=366-60) test calculations are performed. Figure 4 The correlation between the first three principal components and Sq varies with the number of samples. It can be seen that the correlation coefficient between the first principal component and Sq in most tests is greater than 0.8, indicating that the first principal component of NOC can reflect the change of Sq in most cases.

[0052] Figure 5 The data quality analysis diagram is obtained by correlation analysis of the principal component analysis sample data and the average Sq. For the 366 data samples of the Beijing Ming Tombs Geomagnetic Observatory in 2008, the sample number is used for sliding and increasing the sample number for comparison, that is, the sample number is m = 5, 6, 7, ... 366, and the sample is the data sample selected by sliding throughout the year with the sample number as the window and 1 day as the step length. For example, when m = 5, 362 calculations will be performed from the first 5 days of the year to the last 5 days, so the data of the whole year will be sampled and slid 5 to 366 times, and a total of 65340 calculations will be performed. As shown in the figure, the first-order principal component is positively correlated with the data quality, indicating that the better the data quality or the more similar it is to the standard Sq, the more the first-order principal component can reflect the Sq.

[0053] Figure 6 The relationship between the correlation between the first-order principal component and Sq calculated above and the number of data samples and sample quality, where the black points represent the sample quality and the gray points represent the number of samples. It can be seen that there is a significant increasing relationship between the correlation between the first-order principal component and Sq and the number of data samples and sample quality, indicating that in most cases the first-order principal component can reflect Sq and can stably extract the geomagnetic diurnal variation field.

[0054] Figure 7 Schematic diagram of the spatial distribution of global stations and CM model grid points in the embodiment. In the figure, the black points are global geomagnetic stations, the asterisk points are 30°×30° model grid points, and the diamond points represent 15°×30° model grid points.

[0055] Figure 8This is a comparative diagram of the ionospheric equivalent current system obtained by inverting measured data or / and two model data (January 1, 2008, UT=12). It can be seen that the current vortex structure obtained using only measured data is not obvious, and the contribution of individual stations is outstanding. After adding 15°×30° and 30°×30° grid model data, the vortex current structure of the two is significant, and the spatial distribution characteristics are basically the same. However, the vortex current intensity obtained by inverting the 15°×30° grid is generally larger, so it is believed that the 30°×30° model grid can provide stable model data.

[0056] Fig. 9 Schematic diagram for comparing the spatial distribution of ionospheric equivalent current obtained by inversion before and after the introduction of CM model data in the embodiment (January 1, 2008, UT=00, 01, 02, 03 from top to bottom).

[0057] It can be seen that the current vortex structure is more prominent after the introduction of the CM model data, and its spatial structure changes in longitude over time, reflecting the effect of the earth's rotation, indicating that the results of this inversion method are stable and reliable.

[0058] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. An hourly inversion method for an ionospheric equivalent current system using a CM model, characterized in that: The following steps are involved: S1 collects and records the original data of geomagnetic field vectors from global geomagnetic stations; S2 eliminates the background changes of long-term items and trend items in the original data. Specifically, the daily data series is first converted from the world time coordinates to the local time coordinates according to the conversion relationship; then, the long-term items are calculated by the average data from 23:00 to 01:00 local time of the station on the same day and deducted from the data, and the gradient values ​​of the data at 00:00 and 23:00 local time are further used to eliminate the trend items of the data at each local time; then the daily data series is converted from the local time coordinates to the world time coordinate system; S3 processes the single station data to obtain the geomagnetic diurnal variation field component, which is used as the geomagnetic diurnal variation field basic data for hourly inversion of ionospheric equivalent current; S4 uses the CM magnetic field model, inputs the operation parameters, and calculates the model magnetic field data of the ionospheric current magnetic field components at the global grid points; combines the measured data and the model data for grid interpolation to obtain the global geomagnetic diurnal field grid data of the ionospheric equivalent current hourly inversion; S5, based on the spherical harmonic analysis method, separates the internal and external source fields from the global geomagnetic diurnal field grid data, and obtains the external source field corresponding to the ionospheric equivalent current magnetic potential; S6 Inversely calculates the hourly ionospheric equivalent current based on the relationship between magnetic field and current given by the Biot-Saffar law; The conversion relationship in S2 is local time = universal time + geographical longitude / 15; The method for processing single data in S3 is principal component analysis method; S4 The CM magnetic field model operation parameters are: time, longitude and latitude coordinates and geomagnetic activity index Dst; The analysis method described in S5 specifically comprises the following steps: S51. Calculate the grid data of the geomagnetic daily variation field, and finally obtain network points evenly distributed in longitude and latitude on a global scale; S52. Solve the Laplace equation for the geomagnetic field components in the spherical coordinate system to obtain the harmonic superposition of a series of sine and cosine functions containing associated Legendre coefficients; The equivalent relationship is established by the Fourier coefficients of the observed data and the harmonic coefficients of the spherical harmonic expansion, that is, HA = B; S53. Establish the expression of Gaussian coefficient through the relationship between the Fourier harmonic component sine and cosine coefficients and the associated Legendre coefficients, and solve for the Gaussian coefficient; The data calculation method of S51 is: S511 expands the grid data of the geomagnetic daily variation field by Fourier series according to the longitude direction to obtain the Fourier coefficients corresponding to each latitude; S522. Interpolate all Fourier coefficients in the latitude direction to make them continuously smooth within a certain latitude range, thereby obtaining network points uniformly distributed in longitude and latitude on a global scale; In the S52, in the equivalent relationship HA=B, H is a matrix composed of basis functions, and A is , , , The associated Legendre coefficient matrix is ​​composed of the triangular coefficient matrix, and B is the matrix composed of triangular coefficients. It can be obtained that A = H\B, that is, the associated Legendre coefficient is obtained , , , .

2. The hourly inversion method of the ionospheric equivalent current system introducing the CM model according to claim 1 is characterized in that: The method for solving the Gaussian coefficient in S53 is: least square method.

3. The hourly inversion method of the ionospheric equivalent current system introducing the CM model according to claim 2 is characterized in that: The interpolation method in S522 is specifically: interpolation is performed using a cubic spline function in the latitude direction.

4. The hourly inversion method of the ionospheric equivalent current system introducing the CM model according to claim 2 is characterized in that: The S522 mid-latitude range is -90° to 90°.

Images