Gravitational field model establishment method of new-generation low-low satellite tracking satellite scheme based on Chinese region enhancement

Abstract

The invention discloses a time-varying gravity field model establishment method of a new generation of low-low satellite tracking satellite scheme based on China region enhancement, and belongs to the field of satellite gravitational mechanics. The method refers to a selection index of a global optimal double-track four-satellite formation scheme, determines a plurality of double-track four-satellite formation schemes in a China region, selects an optimal scheme from the double-track four-satellite formation schemes, and establishes a time-varying gravity field model of the new generation of low-low satellite tracking satellite scheme. Different data in the optimal scheme are adopted to independently solve low-order and high-order gravitational potential spherical harmonic coefficients, the negative influence of the time domain aliasing effect on the precision of the time-varying gravity field model is weakened, the error of the high-order gravitational potential spherical harmonic coefficients can be reduced to a certain degree, and finally the Chinese region gravity field model with higher precision is obtained. A final result shows that compared with an independent polar orbit formation, the method not only greatly weakens south-north stripe errors caused by a south-north tracking mode in a Chinese region, but also avoids the problems of gravity signal leakage and the like caused by a post-processing process, and finally obtains a Chinese region gravity field model with smaller errors and larger signal-to-noise ratio.

Description

Technical Field

[0001] The present invention belongs to the field of satellite gravity, and more specifically, relates to a method for establishing a gravity field model of a new generation of low-low satellite tracking satellite solution based on China regional enhancement. Background Art

[0002] Satellite gravity technology is an important means to achieve high-precision global gravity field modeling. The gravity satellite missions that have been launched include CNAMP, GRACE, GRACE Follow-On, GOCE, etc. They all use the independent polar formation mode (PSF, Polar Satellite Formation), that is, the angle between the satellite orbit inclination and the equatorial plane is close to 90°, which makes the gravity satellite only sensitive to the north-south gravity signal, but not to the east-west gravity signal, which eventually causes the north-south strip noise in the Earth Gravity Field Model (EGM: Earth Gravity). In order to obtain a more accurate Earth gravity field model, the international community has carried out research on a new generation of satellite gravity detection technology, such as the GRACE-FO satellite, which is based on the GRACE satellite. It is additionally equipped with a laser ranging system to complete the measurement of inter-satellite distance and its changes, which can increase the accuracy of inter-satellite observations by about 100 times; at the same time, a drag-free compensation system can be used to compensate for the non-conservative forces on the satellite in real time. However, under this trend, the accuracy of the solution of the Earth's gravity field model has not been greatly improved. This is because there are still important error sources that affect the accuracy of the gravity field model solution in the new generation of gravity satellite missions, namely, time-domain mixing errors. In addition to the optimization of onboard instruments, many scholars have also conducted corresponding research on the new satellite formation configuration. The above-mentioned gravity satellites have major defects in the satellite formation configuration, that is, the gravity satellite can only detect gravity signals in the north-south direction, but lacks detection of east-west signals; at the same time, since there is only one group or a single satellite, its signal sampling frequency is low. According to the Nyquist sampling theorem, this will make it difficult for the gravity satellite's observation data to correctly model short-period time-varying gravity signals, thereby affecting the accuracy of the gravity field model solution.

[0003] At present, the international community is considering launching a dual-orbit four-star formation gravity satellite program to obtain a more accurate global gravity field model. However, for high-precision regional gravity field models, the development of a satellite program with regional enhancement capabilities has not been fully considered. In order to further improve my country's ability to characterize the Earth's gravity field, it is urgent to develop a new generation of gravity satellite programs for China's regional enhancement. Summary of the invention

[0004] In view of the above defects or improvement needs of the prior art, the present invention provides a method for establishing a gravity field model of a new generation of low-low satellite tracking satellite solution based on China's regional enhancement, which can improve the accuracy of China's regional gravity field model. To achieve the above purpose, according to a first aspect of the present invention, a method for establishing a time-varying gravity field model of a new generation of low-low satellite tracking satellite solution based on China's regional enhancement is provided, comprising:

[0005] S1, determining a tracking satellite plan; wherein the tracking satellite plan includes 6 sets of dual-orbit four-satellite formation plans; in each formation plan, the combinations of polar orbit inclination and inclined orbit inclination are [89° / 55°], [89° / 60°], [89° / 65°], [89° / 70°], [89° / 75°], and [89° / 80°] respectively;

[0006] S2, respectively process the satellite trajectories of each formation scheme to obtain the corresponding time-varying gravity field model;

[0007] S3, taking the formation plan with the best observation uniformity of the gravity satellite and the best accuracy of the time-varying gravity field model as the optimal formation plan; using a preset method to determine the final time-varying gravity field model based on the simulated observation data and reference observation data of the optimal formation plan.

[0008] The preset method includes: according to the high-order solution period n, the high-order maximum solution order n max , low-order maximum solution order n int , using n days of simulated observation data and reference observation data to solve the order n int +1 to n max Spherical harmonic coefficients; according to n int And the low-order solution period T int , respectively using T int The simulated observation data and reference observation data of the day are solved to obtain n groups of 2nd to nth order int The spherical harmonic coefficients of order i are obtained by weighted averaging them to obtain a set of 2nd to i nit The average spherical harmonic coefficient of order i nit +1 to n max The spherical harmonic coefficients of T are used together as the final spherical harmonic coefficients of the time-varying gravity field model; int <n。

[0009] According to a second aspect of the present invention, there is provided an electronic device, comprising: a computer-readable storage medium and a processor;

[0010] The computer-readable storage medium is used to store executable instructions;

[0011] The processor is used to read the executable instructions stored in the computer-readable storage medium and execute the method as described in the first aspect.

[0012] According to a third aspect of the present invention, a computer-readable storage medium is provided, wherein the computer-readable storage medium stores computer instructions, and the computer instructions are used to enable a processor to execute the method described in the first aspect.

[0013] In general, the above technical solutions conceived by the present invention can achieve the following beneficial effects compared with the prior art:

[0014] Considering that when detecting the earth's gravity field based on the existing independent polar orbit formation, there are inevitably systematic errors such as north-south strip errors; the temporal and spatial resolution and accuracy of the current gravity satellites can no longer meet the needs of refined gravity signal detection, and there is a lack of high-precision local regional gravity field models. Based on this, the present invention proposes a gravity field model establishment method based on a new generation of low-low satellite tracking satellite scheme enhanced in the Chinese region. Referring to the selection index of the global dual-orbit four-star formation scheme, multiple dual-orbit four-star formation schemes in the Chinese region are determined, and the optimal scheme is selected from them. Different data in the optimal scheme are used to independently solve the low- and high-order gravity potential spherical harmonic coefficients, weakening the negative impact of the time domain aliasing effect on the accuracy of the time-varying gravity field model, and can reduce the error of the high-order gravity potential spherical harmonic coefficients to a certain extent, and finally obtain a higher-precision Chinese regional gravity field model. The final results show that compared with the independent polar orbit formation, the present invention not only greatly weakens the north-south strip error caused by the north-south tracking mode in the Chinese region, but also avoids the problems of gravity signal leakage caused by the post-processing process, and finally obtains a Chinese regional gravity field model with smaller errors and higher signal-to-noise ratio.

[0015] The method provided by the present invention separately solves the low-order daily solution time-varying gravity signal. In the general gravity field solution process, the same set of data is usually used to simultaneously solve the high- and low-order gravity potential spherical harmonic coefficients, which will cause aliasing effects on signals of different time and space scales, that is, high-order long-period time-varying gravity signals will be seriously contaminated by low-order short-period time-varying signals. The present invention uses different data to independently solve the low- and high-order gravity potential spherical harmonic coefficients, weakens the negative influence of the time domain aliasing effect on the accuracy of the time-varying gravity field model, and can reduce the error of the high-order gravity potential spherical harmonic coefficients to a certain extent. BRIEF DESCRIPTION OF THE DRAWINGS

[0016] Figure 1 A schematic diagram of the flow of gravity field inversion provided for an advanced embodiment of the present invention.

[0017] Figure 2The following are sub-satellite point trajectory diagrams of six dual-orbit four-satellite formation schemes in the Chinese region provided for the further embodiments of the present invention (wherein the red line and the blue line represent the sub-satellite point trajectories of the inclined orbit satellite group and the polar orbit satellite group, respectively).

[0018] Figure 3 A schematic diagram of the changes in terrestrial water reserves in China using the gravity field model (without spatial filtering) of the six dual-track four-satellite formation schemes provided in the embodiments of the present invention and the “HIS” model.

[0019] Figure 4 A schematic diagram of the interannual trend signal fitting of the "HIS" model and the changes in land water reserves in the Chinese region corresponding to six dual-track four-star formation schemes (without spatial smoothing filtering) are provided for the further implementation of the present invention.

[0020] Figure 5 The present invention provides six dual-track four-star formation schemes corresponding to the changes in land water reserves in China (without spatial smoothing filtering) and a schematic diagram of the interannual amplitude signal fitting of the "HIS" model.

[0021] Figure 6 A schematic diagram of the fitting residuals of the "HIS" model and the changes in terrestrial water reserves in the Chinese region corresponding to six dual-track four-star formation schemes are provided for the further implementation of the present invention (without spatial smoothing filtering).

[0022] Figure 7 A schematic diagram of the RMS (Root Mean Square) difference between the changes in land water reserves in the Chinese region corresponding to six dual-track four-star formation schemes (without spatial smoothing filtering) and the "HIS" model is provided for the further implementation of the present invention.

[0023] Figure 8 A flow chart of determining the final time-varying gravity field model based on simulated observation data and reference observation data of the optimal formation scheme provided in an embodiment of the present invention. DETAILED DESCRIPTION

[0024] In order to make the purpose, technical solutions and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not intended to limit the present invention. In addition, the technical features involved in the various embodiments of the present invention described below can be combined with each other as long as they do not conflict with each other.

[0025] The embodiment of the present invention provides a method for establishing a gravity field model of a new generation low-low satellite tracking satellite solution based on China regional enhancement, comprising:

[0026] S1, determine a tracking satellite plan; wherein the tracking satellite plan includes 6 sets of dual-orbit four-star formation plans; in each formation plan, the combinations of polar orbit inclination and inclined orbit inclination are [89° / 55°], [89° / 60°], [89° / 65°], [89° / 70°], [89° / 75°], [89° / 80°] respectively.

[0027] Specifically, the dual-orbit four-satellite formation is based on the independent polar orbit formation, and adds a group of inclined orbit satellites, which not only increases the detection of gravity satellites to the east-west gravity signal, but also increases the sampling frequency of satellites in the repeated observation area. In addition, if the global scope is used as the application scale, [89° / 70°] is the optimal dual-orbit four-satellite formation solution.

[0028] Based on this, the orbital inclination of the polar-orbiting satellite group in the dual-orbit four-satellite formation is set to 89° to ensure that observation data are available worldwide. Considering that the northernmost latitude in the Chinese region is about 54°, and the observation range of the inclined-orbit satellite depends on its orbital inclination, and considering that the orbital inclination range of the inclined-orbit satellite group must be [55°, 80°], otherwise the inversion performance of the dual-orbit four-satellite formation will be reduced. In view of the above two points, the orbital inclination selection range of the inclined-orbit satellite group is determined to be [55°, 80°] to ensure that the inclined-orbit satellite has observation data in the Chinese region.

[0029] The orbital inclination interval of the inclined-orbit satellite group is 5°, specifically: [55°, 60°, 65°, 70°, 75°, 80°], which forms a dual-orbit four-star formation with the polar-orbit satellites, thus obtaining 6 dual-orbit four-star formation schemes: [89° / 55°], [89° / 60°], [89° / 65°], [89° / 70°], [89° / 75°], [89° / 80°].

[0030] According to the orbit height h and orbit inclination i and the satellite orbit repetition period (Considering the limitation of the temporal and spatial resolution of current gravity satellites, the present invention sets it to 30 days) The relationship confirms that the orbital altitude of the satellite is h = aR e , to obtain the satellite's orbital semi-major axis a, the calculation formula for a is as follows:

[0031]

[0032] Among them, α is the number of repeated days of satellite flight, β is the number of repeated circles of satellite flight (its geometric meaning is the number of intersections of the satellite flight trajectory and the equator), n is the average angular velocity of the orbit, ω e is the average angular velocity of the Earth, R e is the average radius of the earth, a is the semi-major axis of the satellite orbit, i is the orbital inclination of the satellite, J 2 =-C2,0 (Here C 2,0 is regularized, C 2,0 are the spherical harmonic coefficients that characterize the Earth’s oblateness).

[0033] Finally, taking 30 days as the inversion period of the gravity field model (that is, the time scale of the gravity field model is 30 days), it is determined that the orbital altitudes of the inclined orbit satellites corresponding to [55°, 60°, 65°, 70°, 75°, 80°] are [511.0km, 495.0km, 501.5km, 497.5km, 494.5v, 499.5km] respectively, and the orbital altitude corresponding to the polar orbit satellite group (89°) is 496.5v.

[0034] S2, processing all satellite trajectories of each formation scheme respectively to obtain the gravity field model corresponding to each formation scheme.

[0035] Specifically, all satellite trajectories of each formation scheme are processed based on the orbit integrator to obtain the gravity field model corresponding to each formation scheme.

[0036] Preferably, the processing includes: dividing the flight trajectory of the satellite within a preset time into N arc segments evenly according to a preset time length, inputting the initial state of each arc segment into the first and second orbit integrators respectively, and obtaining simulated observation data and reference observation data of the orbit of the satellite within the preset time; and calculating the second to nth order of the gravity field model by the least squares adjustment method according to the simulated observation data and the reference observation data. max Order spherical harmonic coefficients; wherein the background force model of the first orbit integrator includes a static gravity field model, an ocean tide model and a time-varying gravity field model, and the background force model of the second orbit integrator includes a static gravity field model, an ocean tide model and a non-tidal atmosphere and ocean model.

[0037] Taking into account the integration error, the preset duration is usually 6 hours.

[0038] The preset time can be set according to actual needs.

[0039] The embodiment of the present invention is described by taking the preset time of 60 months and the preset duration of 6 hours as an example.

[0040] Specifically, step S2 includes:

[0041] S21, evenly divide the flight trajectory of the satellite within the preset time into N arc segments according to the preset duration, and calculate the initial state vector of each orbital arc segment of the two groups of satellites in each scheme based on the orbital parameters of each scheme.

[0042] Since the integration error will increase with the growth of time, the orbit will gradually deviate from the set orbital parameters. Therefore, the present invention will predetermine the initial state of each orbital arc segment, and the above initial state is obtained by the third Gauss-Jackson integrator (i.e., the third orbital integrator), wherein the background force model of the integrator is set to the gravity of the center of the earth, so as to ensure that the integrated orbit is always the pre-set orbit at the initial position of each arc segment (here in actual practice, this can be achieved by adopting a resistance compensation system).

[0043] The calculation method of the initial state vector of each orbital arc segment is as follows: the initial mean anomaly angle of the satellite (i.e., the starting position of the first arc segment) is set in the third Gauss-Jackson integrator (i.e., the third orbital integrator). Taking the preset duration of 6 hours and the preset time of 60 months as an example, the flight duration of the two groups of satellites in one day is divided into 4 segments. Accordingly, the flight trajectories of the two groups of satellites within 60 months are divided into 4×60×30=7200 orbital arc segments. Determine the six Kepler orbital elements of the initial state of the two groups of satellites for a total of 7200 orbital arcs in one month (including the semi-major axis, inclination, flattening, right ascension of the ascending node, perigee angular distance and mean anomaly of the satellite orbit, where the mean anomaly is used to determine the position of the satellite on the orbit), and use 6j and 5s as the arc length unit and integration step of the orbital integration respectively; after setting the satellite orbital parameters (i.e. the six Kepler orbital elements of the initial state of each arc segment), the integrator will use 6h and 5s as the arc length unit and integration step of the orbital integration respectively based on the pre-set background force model, and obtain the six Kepler orbital elements of the initial state of each arc segment except the starting arc segment, and convert it into a vector form, i.e. the initial state vector. Among them, 6h is used as the arc length unit of the orbital integration, so as to control the growth of the orbital integration cumulative error.

[0044] S22, gravity field inversion is performed on each of the six schemes to obtain a monthly time-varying gravity field model (displayed in the form of spherical harmonic coefficients), and the initial state vector of the orbital arc is used as the initial value of the orbital integral to invert the time-varying gravity field model for a total of 60 months from 2003 to 2007, and the time-varying gravity field model is converted into terrestrial water storage changes (Terrestrial Water Storage Changing) according to formula (2), which is specifically expressed in the form of equivalent water height (EWH: cm). Equivalent water height means using equivalent water thickness changes to represent gravity signals of corresponding sizes. In the present invention, the equivalent water height is a water column of 1°×1°, and one grid point corresponds to one water column, i.e., the equivalent water height.

[0045]

[0046] Among them, ρ ave represents the average density of the Earth, ρ w Indicates water density, N maxand m represent the cutoff order and degree of the spherical harmonic coefficients, respectively, k l represents the Love number, is the fully normalized Legendre function of order l and order m, C lm and C lm represents the solved l-order m-order gravitational potential spherical harmonic coefficients, and represents the fully normalized l-th order m-th gravitational potential spherical harmonic coefficient, θ and λ represent the co-latitude and longitude of the grid point, and μ represents the gravitational potential spherical harmonic coefficient corresponding to the grid point at (θ, λ).

[0047] The above gravity field inversion includes forward process and inversion process, such as Figure 1 shown.

[0048] In the forward modeling process, based on the initial state vector of the orbital arc segment, the first Gauss-Jackson integrator (i.e., the first orbital integrator) is used to generate simulated values ​​of orbital observation data (including orbital coordinate data of the two satellites at each sampling point and inter-satellite distance variation data), which are called "real observation data." It can be understood that after setting the initial state vector of the orbital arc segment, the first orbital integrator will use 6h and 5s as the arc length unit and integration step of the orbital integration based on the pre-set background force model to obtain the six Kepler orbital numbers of each arc segment sampling point, i.e., "real observation data." Among them, the background force model of the first orbital integrator is a set of "real models."

[0049] The "real model" in the forward simulation process mainly includes the static gravity field model, the ocean tide model and the time-varying gravity field model, which are GGM05s, EOT11a and ESM AOHIS respectively. The time-varying gravity field model ESM AOHIS includes components such as atmosphere ("A"), ocean ("O"), hydrology ("H"), ice ("I") and solid earth ("S"), and its time resolution is 6h. The simulation sampling interval and integral arc length are set to 5s and 6h respectively.

[0050] Furthermore, in order to simulate more realistic observation data, it is necessary to add noise from sensors such as laser interferometers, accelerometers, and GNSS receivers to the error-free simulated observation data. During the simulation process, the observation data of the laser interferometer is obtained by differentially analyzing the satellite orbit data (including position and velocity), which is the simulated satellite observation data; the observation data of the accelerometer is obtained by using the observation data of the accelerometer during the actual flight of the satellite, which is used to simulate the non-conservative force on the satellite during flight; the observation data of the GNSS receiver corresponds to the satellite orbit data (including position and velocity). These noise-contaminated observation data are called "simulated observation data."

[0051] In the inversion process, based on the initial state vector of the orbital arc segment, the second Gauss-Jackson integrator (i.e., the second orbital integrator) is used to generate simulated values ​​of the orbital observation data (including the orbital coordinate data of the two satellites at each sampling point and the inter-satellite distance variation rate data), which are called "reference observation data". It can be understood that after setting the initial state vector of the orbital arc segment, the second orbital integrator will use 6h and 5s as the arc length unit and integration step of the orbital integration based on the pre-set background force model to obtain the six Kepler orbit numbers of each arc segment sampling point, i.e., "reference observation data". Among them, the background force model of the second orbital integrator is a set of "reference models".

[0052] The "reference model" in the above inversion process mainly includes a static gravity field model, an ocean tide model and a non-tidal atmosphere-ocean model. Among them, the static gravity field model and the ocean tide model respectively adopt the GGM05s model and the FES2004 model. For the non-tidal atmosphere-ocean model (AOD), the present invention adopts the DEAL model and the AOerr model, which correspond to the signals and model errors of the non-tidal mass changes of the time-varying gravity signal components of the atmosphere (A) and the ocean (O), respectively.

[0053] Based on the above “simulated observation data” and “reference observation data”, the time-varying gravity field model is estimated using the least squares adjustment method:

[0054] The least squares adjustment observation equation is established, and the observation vector l, design matrix A and unknown vector x are established in turn, as shown in formula (2), to obtain the spherical harmonic coefficients of the time-varying gravity field model. That is, the total amount of the three time-varying signal components of hydrology ("H"), ice ("I") and solid soil ("S") in the time-varying gravity field model ESMAOHIS is estimated, referred to as "HIS".

[0055] Among them, the observation vector l represents the residual between the “real observation data” and the “reference observation data”;

[0056] The design matrix A is composed of the variational equation (i.e., the partial derivative equation of the satellite instantaneous state vector with respect to the spherical harmonic coefficients to be solved) The derivative with respect to time is a function of The solution of the first-order ordinary differential equation of is obtained by numerical integration;

[0057] The unknown vector x includes the local parameters and global parameters of each arc segment (i.e., the spherical harmonic coefficients to be solved representing the earth's gravity field model), which are calculated by formula (3), where P represents the weight matrix (in the present invention, P is the unit matrix). The local parameters include the scale factor and deviation correction of the accelerometer and the position and velocity correction of the initial state of each arc segment, and the x, y, and z directions all need to be corrected. In the present invention, the scale factor is set to 1, and the deviation includes the constant term, the first term, and the second term correction, so there are actually 15 local parameters that need to be corrected.

[0058] Finally, a set of spherical harmonic coefficients x representing the earth's time-varying gravity signal is calculated according to the principle of least square adjustment. The spherical harmonic coefficients calculated by the present invention are of order 2-60.

[0059] Ax≈l (3)

[0060] x=(A T PA) -1 A T Pl (4)

[0061] The above calculation process is specifically shown as follows.

[0062] Since the observation vector l and its corresponding Ax are not completely equal, a correction V needs to be added, as shown in formula (5):

[0063] V=Ax-l (5)

[0064] The criterion for least squares adjustment is:

[0065] V T PV=min (6)

[0066] According to the least squares principle, x must satisfy the above criteria, so the partial derivative of the above formula with respect to x is obtained:

[0067]

[0068] Substituting formula (5) into formula (7), we can get:

[0069] A T PAx-A T Pl=0 (8)

[0070] Finally, x is calculated as:

[0071] x=(A T PA) -1 A T Pl (9)

[0072] S3, taking the formation scheme with the best observation uniformity of the gravity satellite and the best accuracy of the gravity field model as the optimal formation scheme; using a preset method to determine the final gravity field model according to the simulated observation data and reference observation data of the optimal formation scheme;

[0073] The preset method includes: according to the high-order solution period n, the high-order maximum solution order n max , low-order maximum solution order n int , using n days of simulated observation data and reference observation data to solve the order n int +1 to n max Spherical harmonic coefficients; according to n int And the low-order solution period T int , respectively using T int The simulated observation data and reference observation data of the day are solved to obtain n groups of 2nd to nth order int The spherical harmonic coefficients of order 2 are obtained by weighted averaging. int The average spherical harmonic coefficients of order n int +1 to n max The spherical harmonic coefficients of T are used together as the final spherical harmonic coefficients of the time-varying gravity field model; int <n。

[0074] The observation uniformity of the gravity satellite is mainly reflected in the sub-satellite point trajectory of the dual-track four-satellite formation scheme in the Chinese region.

[0075] The sub-satellite point trajectory diagrams of six dual-track four-star formation schemes in the Chinese region are drawn, and the distribution uniformity of each scheme in the Chinese region, that is, the observation uniformity, is compared.

[0076] The accuracy of the gravity field model is determined by:

[0077] The smaller the difference between the gravity field model and the "HIS" model, the higher the model accuracy of the gravity field.

[0078] Specifically, the accuracy level of the time-varying gravity field model can be reflected by comparing the differences between the terrestrial water storage changes of the dual-track four-star formation scheme and the "HIS" model during the five years from 2003 to 2007.

[0079] The terrestrial water storage changes of each scheme in the five years from 2003 to 2007 were compared with the "HIS" model (it can be understood that the HIS model is a set of spherical harmonic coefficients, which can be converted into equivalent water heights for comparison with terrestrial water storage changes) to evaluate the accuracy of the time-varying gravity field model. The specific contents are: 1) Calculate the fitted interannual amplitude term and fitted trend term of terrestrial water storage changes, and compare them with the fitted interannual amplitude term and fitted trend term of the "HIS" model; 2) Calculate the RMS difference between the terrestrial water storage changes and the "HIS" model.

[0080] The closer the fitted interannual amplitude term and trend term are to the “HIS” model, and the smaller the RMS difference from the “HIS” model, the higher the accuracy and more reliable the gravity field model.

[0081] The fitted interannual amplitude term and fitted trend term of terrestrial water storage change are calculated based on formula (8):

[0082] TWS=a+bt+c 1 cos(2πt)+c 2 sin(2πt)+e 1 cos(4πt)+e 2 sin(4πt) (8)

[0083] Among them, b is the linear trend, that is, the fitting trend term, (c 1 cos(2πt)+c 2 sin(2πt)+e 1 cos(4πt)+e 2 sin(4πt)) is the fitting annual and semi-annual amplitude terms, c 1 and c 2 represents the phase information of the annual time-varying gravity signal, e 1 and e 2 It represents the phase information of the semi-annual time-varying gravity signal. The calculation process of the RMS difference between the terrestrial water storage change at each grid point and the "HIS" model is as follows: first remove the equivalent water medium amplitude (including annual amplitude and semi-annual amplitude) and trend term of each grid point to obtain its monthly residual value, then calculate the RMS value of 60 months, and finally obtain the residual RMS value of each grid point.

[0084] After obtaining the optimal dual-track four-star formation plan for the Chinese region, it is necessary to further adopt the preset method to weaken the negative impact of the time-domain aliasing effect on the time-varying gravity field model.

[0085] The above preset method is specifically as follows: the total arc length of the satellite's flight for one day is taken as the minimum solution object, that is, one day is taken as the unit of data used for solution, and the low-order daily solution gravity potential spherical harmonic coefficients and the correction amount of the satellite's initial state vector (including position and velocity), the scale factor of the accelerometer and the correction amount of the deviation are taken as local parameters in the gravity field solution process for solution.

[0086] The solution period n (n is also the time scale of the time-varying gravity field model, which is set to 30 days in this example) and the solution order n of the high-order long-period solution are max 、The calculation period of low-order daily solution T int Solve for order n int Assignment; According to the low-order daily solution of gravity potential spherical harmonic coefficients, the order n int , using n days of "real observation data" and "reference observation data" to solve the high-order gravity field model, that is, to solve the n-order int +1 to n max The spherical harmonic coefficients of the gravity field model are then solved using the "real observation data" and "reference observation data" of one day in n days to obtain n groups of daily solutions from 2nd to nth order. int The gravity potential coefficient (for example, when using a single-day solution, T int =1). The weighted average of n sets of daily solution low-order gravity potential spherical harmonic coefficients is used to obtain a set of 2nd to nth order spherical harmonic coefficients. int The gravitational potential spherical harmonic coefficients are finally compared with the n-day average high-order n int +1 level to n max The spherical harmonic coefficients of the gravity potential are composed of n-day average 2nd to nth order max The spherical harmonic coefficients of the gravitational potential.

[0087] For example, after comparing the time-varying gravity field model inversion capabilities of different schemes, the optimal dual-track four-satellite formation scheme is determined to be [89° / 55°]. Based on this, the present invention further adopts a preset method to weaken the time domain aliasing error. int 、n max And the solution period T of the low-order daily solution int (e.g. one day), where n corresponds to the time scale of the time-varying gravity field model. In the present invention, the time scale is the monthly scale, so the parameter n is assigned a value of 30; at the same time, based on the excellent performance of the dual-track four-star formation in noise reduction, n max In the present invention, it is set to 60 (no spatial filtering is performed). int and T int The optimal value of n is unknown. int and T int Set up different experimental groups and use the n of each experimental group int and T intPerform anti-aliasing processing and compare the solution accuracy of the time-varying gravity field model of each group of experiments to determine n int and T int The optimal value of n int and T int The optimal value is used to perform anti-aliasing processing in the process of solving the time-varying gravity field, thereby further improving the solution accuracy of the time-varying gravity field model.

[0088] Assume that n int and T int The optimal values ​​are 15 and 1 respectively. The 16th to 60th order spherical harmonic coefficients are calculated using the 30-day "real observation data" and "reference observation data" in the optimal dual-orbit four-star formation scheme, and then 30 groups of 2nd to 15th order spherical harmonic coefficients are calculated using the "real observation data" and "reference observation data" of each day in the 30 days. A group of 2nd to 15th order average spherical harmonic coefficients are weighted averaged, and are used together with the above-mentioned 16th to 60th order spherical harmonic coefficients as the final spherical harmonic coefficients of the time-varying gravity field model.

[0089] It can be understood that the embodiment of the present invention integrates the satellite's 60-month flight trajectory to obtain the corresponding "real observation data" and "reference observation data". Since the time scale of the time-varying gravity field model is 30 days, when solving the spherical harmonic coefficients of the time-varying gravity field model, the "real observation data" and "reference observation data" of the 1st to 60th months (each month has about 30 days, and the actual number of days is set according to the number of days in the month) are used to solve the spherical harmonic coefficients of the time-varying gravity field model, and finally the spherical harmonic coefficients of the time-varying gravity field model of 1-60 months are obtained.

[0090] Figure 2 The sub-satellite point trajectory diagram of the six dual-orbit four-satellite formation schemes in the Chinese region provided for the embodiments of the present invention (the red line and the blue line represent the sub-satellite point trajectory of the inclined orbit satellite group and the polar orbit satellite group, respectively). It can be seen that compared with the independent polar orbit formation, the six schemes all increase the observation in the east-west direction; among the six schemes, the sub-satellite point trajectory of the scheme [89.0° / 55.0°] is the most dense and uniform, and the observation uniformity is the best. Figure 3 The gravity field model (without spatial filtering) of the six dual-orbit four-satellite formation schemes provided in the embodiment of the present invention and the change of terrestrial water storage in China under the "HIS" model. It can be seen that the gravity field model corresponding to the independent polar orbit formation has the largest strip noise, and the gravity field signal is almost invisible; in the six dual-orbit four-satellite formation schemes, as the orbital inclination of the inclined orbit satellite group gradually increases, the time-varying gravity field model gradually shows strip noise. Figure 4The present invention provides six dual-track four-star formation schemes corresponding to the changes in land water storage in China (without spatial smoothing filtering) and the interannual trend signal fitting of the "HIS" model. It can be seen that among the six schemes, the schemes [[89° / 55°], [89° / 60°] are closest to the "HIS" model; Figure 5 The present invention provides six dual-track four-star formation schemes corresponding to the changes in land water storage in China (without spatial smoothing filtering) and the interannual amplitude signal fitting of the "HIS" model. It can be seen that among the six schemes, the schemes [[89° / 55°], [89° / 60°] are closest to the "HIS" model; Figure 6 The present invention provides six dual-track four-star formation schemes corresponding to the changes in land water storage in China (without spatial smoothing filtering) and the fitting residuals of the "HIS" model. It can be seen that the fitting residual of the scheme [89° / 55°] is the smallest; Figure 7 The present invention provides the RMS difference between the land water storage changes in China (without spatial smoothing filtering) and the "HIS" model for the six dual-track four-star formation schemes. It can be seen that the difference between the scheme [89° / 55°] and the "HIS" model is the smallest.

[0091] An embodiment of the present invention provides an electronic device, including: a computer-readable storage medium and a processor;

[0092] The computer-readable storage medium is used to store executable instructions;

[0093] The processor is used to read the executable instructions stored in the computer-readable storage medium and execute the method described in any of the above embodiments.

[0094] An embodiment of the present invention provides a computer-readable storage medium, wherein the computer-readable storage medium stores computer instructions, and the computer instructions are used to enable a processor to execute any of the above methods.

[0095] It will be easily understood by those skilled in the art that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for establishing a time-varying gravity field model based on a new generation of low-low satellite tracking satellite scheme based on China regional enhancement, characterized in that: include: S1, determining a tracking satellite plan; wherein the tracking satellite plan includes 6 sets of dual-orbit four-satellite formation plans; in each formation plan, the combinations of polar orbit inclination and inclined orbit inclination are [89° / 55°], [89° / 60°], [89° / 65°], [89° / 70°], [89° / 75°], and [89° / 80°] respectively; S2, respectively process the satellite trajectories of each formation scheme to obtain the corresponding time-varying gravity field model; S3, taking the formation scheme with the best observation uniformity of the gravity satellite and the best accuracy of the time-varying gravity field model as the optimal formation scheme; using a preset method to determine the final time-varying gravity field model according to the simulated observation data and reference observation data of the optimal formation scheme; The preset method includes: according to the high-order solution period n, the high-order maximum solution order n max , low-order maximum solution order n int , using n days of simulated observation data and reference observation data to solve the order n int +1 to n max Spherical harmonic coefficients; according to n int And the low-order solution period T int , respectively using T int The simulated observation data and reference observation data of the day are solved to obtain n groups of 2nd to nth order int The spherical harmonic coefficients of order 2 are obtained by weighted averaging. int The average spherical harmonic coefficients of order n int +1 to n max The spherical harmonic coefficients of T are used together as the final spherical harmonic coefficients of the time-varying gravity field model; int <n。 2. The method according to claim 1, characterized in that T int and n int The way to obtain is: Setting T int and n int Multiple value combinations of T init and n int The value combinations of T and n max , solve the 2nd to nth order time-varying gravity field model max The spherical harmonic coefficients of order 1 are used to select the combination with the best accuracy of the time-varying gravity field model as T int and n int The final value of .

3. The method according to claim 1, characterized in that In step S2, the processing includes: evenly dividing the flight trajectory of the satellite within a preset time into N arc segments according to a preset time length, inputting the initial state of each arc segment into the first and second orbit integrators respectively, and obtaining simulated observation data and reference observation data of the satellite's orbit within the preset time; according to the simulated observation data and the reference observation data, using the least squares adjustment method to solve the spherical harmonic coefficients of the time-varying gravity field model; wherein the background force model of the first orbit integrator includes a static gravity field model, an ocean tide model and a time-varying gravity field model, and the background force model of the second orbit integrator includes a static gravity field model, an ocean tide model and a non-tidal atmosphere and ocean model.

4. The method according to claim 3, characterized in that The initial state of each arc segment is obtained as follows: Set the initial mean anomaly of the satellite, which is the initial state of the starting arc segment, and input it into the third orbit integrator to obtain the initial states of other arc segments; Wherein, the background force model of the third orbit integrator is the earth center gravity model.

5. The method according to claim 1 or 3, characterized in that: The initial state of each arc segment is characterized by the six Kepler orbit numbers.

6. The method according to claim 3, characterized in that The static gravity field model, ocean tide model and time-varying gravity field model in the first orbit integrator adopt the GOCO05s model, the EOT11a model and the ESM AOHIS model respectively; The static gravity field model and ocean tide model in the second orbital integrator adopt the GGM05s model and the FES2004 model respectively, the non-tidal atmosphere and ocean model adopts the DEAL model, and the error of the model adopts the AOerr model.

7. The method according to claim 1, characterized in that The simulated observation data also includes noise data of the inter-satellite ranging system, accelerometer and GNSS receiver.

8. The method according to claim 1, characterized in that The accuracy of the time-varying gravity field model is determined by: The smaller the difference between the time-varying gravity field model and the "HIS" model, the higher the accuracy of the gravity field model; among them, the "HIS" model is the total amount of the three time-varying signal components of "H", "I" and "S" in the ESM AOHIS model.

9. An electronic device, characterized in that: include: A computer readable storage medium and a processor; The computer-readable storage medium is used to store executable instructions; The processor is used to read the executable instructions stored in the computer-readable storage medium and execute the method according to any one of claims 1 to 8.

10. A computer-readable storage medium, characterized in that: The computer-readable storage medium stores computer instructions, and the computer instructions are used to enable a processor to execute the method according to any one of claims 1 to 8.

Images