Method for constructing geoid based on least squares spectral combination and spherical wavelet

Abstract

The application discloses a geoid construction method based on least square spectrum combination and spherical wavelet, and specifically comprises the following steps: obtaining multi-source gravity data of a region to be constructed, including global gravity field model data, aviation gravity data and ground gravity data; then, according to signal and error spectrum characteristics of each data source, effectively contributing to low-frequency, medium-frequency and high-frequency components is respectively attributed; first, a spectrum weight function is constructed based on a modified spherical Butterworth function form to accurately define low-frequency, medium-frequency and high-frequency bands; then, the determined spectrum weight is used to weight model multi-source gravity observation data, and low-frequency, medium-frequency and high-frequency geoid component models are respectively established; first, a zero-order term and a conversion parameter are obtained, then the optimal weight of the low-frequency component, the medium-frequency component and the high-frequency component, the zero-order term and the conversion parameter are combined to obtain a final geoid result. The application has good stability in geoid modeling in a complex terrain region.

Description

Technical Field

[0001] This invention relates to an improvement in geoid construction technology, belonging to the field of geodesy, and particularly to a geoid construction method based on least squares spectral combination and spherical wavelet. Background Technology

[0002] The geoid is one of the geodetic benchmarks. Determining the geoid is an important project in basic surveying. The shape of the geoid reflects information such as the structure, density, and distribution of materials inside the Earth, and plays an important role in research and application in related Earth science fields such as oceanography, seismology, geophysics, geological exploration, and petroleum exploration.

[0003] Multi-source gravity data, as the core foundation of geoid modeling, inherently suffers from physical inconsistencies in the spectral domain. Satellite gravity primarily provides low-frequency (long-wavelength) information, airborne gravity mainly covers the mid-frequency band, and ground gravity provides high-frequency information. These three types of data differ significantly in terms of observation noise, coverage density, resolution, and spectral energy attenuation patterns. Directly stitching and fusing such data can easily lead to spectral leakage, spectral discontinuities, and noise amplification, making it difficult to guarantee the continuity and physical consistency of regional geoid models. Furthermore, complex terrain areas such as mountainous regions, canyons, coastal areas, and islands exhibit significant short- to medium-scale geoid energy variations. Relying solely on satellite gravity data, or the traditional approach of combining satellite and ground gravity data, cannot effectively recover the mid-frequency information of these areas, resulting in poor stability of geoid modeling in complex terrain regions.

[0004] Chinese patent application CN202311477549.5, filed on November 8, 2023, discloses a regional geoid refinement method and system based on the least squares collocation method. This method subtracts normal gravity, free-space anomalies, and gravity anomalies calculated by a global gravity field model from measured regional gravity data. The residual gravity anomaly values ​​at each point are obtained by subtracting the topographic gravity calculated using a multi-core parallel computation based on prism integrals from a high-resolution digital elevation model from the above results. The least squares collocation method and Hirvonen covariance are then used. The model calculates autocovariance and crosscovariance, and fits a smooth residual gravity anomaly grid. The residual gravity anomaly grid is then calculated using Stokes integration or fast Fourier transform to obtain a residual elevation anomaly grid. The regional refined geoid is obtained by adding the elevation anomaly and the elevation influenced by topographic gravity to the residual elevation anomaly grid. The above scheme can more accurately achieve grid interpolation and error estimation through variance and covariance models, and achieve optimal geoid estimation in the region. However, the above scheme does not solve the problem of poor stability in geoid modeling in complex terrain areas.

[0005] The information disclosed in this background section is intended only to enhance the understanding of the overall background of this patent application and should not be construed as an admission or in any way implying that the information constitutes prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to overcome the shortcomings of existing geoid modeling methods in complex terrain areas, which suffer from insufficient model stability due to incomplete recovery of mid-to-high frequency gravity signals. This invention provides a geoid construction method using least squares spectrum combination and spherical wavelet, which offers better model stability in complex terrain areas where mid-to-high frequency gravity signal recovery is incomplete.

[0007] To achieve the above objectives, the technical solution of the present invention is: a method for constructing a geoid based on least squares spectral combination and spherical wavelet, wherein the method for constructing a geoid based on least squares spectral combination and spherical wavelet includes the following steps:

[0008] The first step is to acquire multi-source gravity data for the region to be constructed, including global gravity field model data, airborne gravity data, and ground gravity data; then, based on the signal and error spectrum characteristics of each data source, their effective contributions are attributed to low-frequency, mid-frequency, and high-frequency components, respectively.

[0009] The second step is to first construct a spectral weighting function based on the modified spherical Butterworth function to accurately define the low-frequency, mid-frequency and high-frequency bands; then, using the determined spectral weights, perform weighted modeling on the multi-source gravity observation data to establish low-frequency, mid-frequency and high-frequency geoid component models respectively.

[0010] The third step is to first obtain the zero-order term and transformation parameters, and then combine the optimal weights of the low-frequency, mid-frequency, and high-frequency components with the zero-order term and transformation parameters to obtain the final geoid result.

[0011] After obtaining the final geoid result, the accuracy verification is also included.

[0012] The accuracy verification specifically involves using GNSS leveling data to check the accuracy of the constructed geoid model, with an accuracy of 1-2 cm in plain areas and 2-3 cm in mountainous or plateau areas.

[0013] The multi-source gravity data also includes terrain data and GNSS leveling data.

[0014] The low-frequency, mid-frequency, and high-frequency components are specifically:

[0015] The range of spherical harmonic orders corresponding to the low-frequency components is: 2≤n≤190;

[0016] The range of spherical harmonic orders corresponding to the intermediate frequency components is: 190 < n ≤ 545;

[0017] The range of spherical harmonic orders corresponding to high-frequency components is: n > 545;

[0018] The dominant data source for the low-frequency component is global gravity field model data, the dominant data source for the mid-frequency component is airborne gravity data, and the dominant data source for the high-frequency component is ground gravity data.

[0019] The low-frequency geoid component model is expressed in a spherical harmonic expansion, specifically as follows:

[0020] ;

[0021] in, The distance to the Earth's center. The gravitational constant is the constant of gravity. For reference to Earth's radius, and Let the order and degree of the spherical harmonic expansion be given. The maximum number of low-frequency components. and These are the fully normalized spherical harmonic coefficients. For a fully normalized association Legendre function, , These are spherical latitude and longitude, respectively.

[0022] The corresponding gravitational perturbation is obtained by radial differentiation of the perturbation potential:

[0023] ;

[0024] This section describes the long-wave gravity characteristics and serves as a benchmark for the enhancement of mid-to-high frequency components.

[0025] The mid-frequency geoid component modeling is specifically as follows: To fill the spectral gap between satellite and ground gravity measurements, airborne gravity data is used to extract mid-band disturbances. First, the low-frequency component and the topographic residual model are removed from the airborne gravity observations to obtain the mid-frequency component residual field.

[0026] ;

[0027] in, This represents the residual gravity perturbation in the mid-frequency component. For airborne gravity disturbance observation, For the middle arrive RTM correction;

[0028] To achieve optimal spectral localization, the mid-frequency component perturbation is represented using spherical wavelet expansion:

[0029] ;

[0030] in, To control the scale parameters of spatial resolution. For the position index of the wavelet center, For scale The number of wavelet basis functions corresponding to a location. The wavelet expansion coefficients are given, and the wavelet basis functions are defined as follows:

[0031] .

[0032] The high-frequency geoid component model is as follows: the high-frequency component reflects local density changes and is obtained by processing ground gravity observations. The high-frequency component residual gravity disturbance is expressed as follows:

[0033] ;

[0034] in, For ground gravity disturbance observation, For the middle to RTM correction;

[0035] High-frequency signals are represented using band-limited wavelet expansion:

[0036] ;

[0037] in, These are the wavelet expansion coefficients. For high-frequency wavelet basis functions:

[0038] .

[0039] In the second step, the determined spectral weights are specifically: to ensure that the sum of the spectral weights of the three frequency bands is always 1.

[0040] ;

[0041] in, , and These are the spectral weights for GGM, airborne and ground gravity disturbances, respectively.

[0042] Spectral weights are calculated using the following formula:

[0043] ;

[0044] in, and It is calculated based on an improved spherical Butterworth filter;

[0045] ;

[0046] in, .

[0047] In the third step, the zero-order term and transformation parameters are first obtained. Then, the optimal weights of the low-frequency, mid-frequency, and high-frequency components are combined with the zero-order term and transformation parameters to obtain the final geoid result. Specifically, the perturbation potential function can be expressed as a weighted superposition of multi-band spectral components.

[0048] ;

[0049] in, , and Represent the perturbation potential spectral components in the low, medium, and high frequency bands, respectively, and the weighting function... To achieve optimal spectral combination, a rationally designed spectral weighting function is used to ensure that the contributions of each frequency band are evenly distributed across different spatial scales.

[0050] The geoid height can be obtained from the disturbance potential:

[0051] ;

[0052] in, Calculated by GGM, Calculated from airborne gravity data, Calculated from land gravity data, The zero-degree term, The difference between the geoid and the quasi-geoid is expressed as:

[0053] ;

[0054] The wavelet basis function with spectral weights is defined as:

[0055] ;

[0056] in, and These are wavelet kernel functions for the frequency bands;

[0057] ;

[0058] in, The gravitational constant of the Earth's center of mass. for The gravitational constant of the reference ellipsoid, The average Earth radius, The zero-degree geoid term is used to explain the difference between the reference Earth model and the actual gravity field. To reference the normal gravity on the surface of the ellipsoid, The position of the geoid. For reference ellipsoid constant, and For the fully normalized spherical harmonic coefficients of the perturbation position, To fully normalize the associated Legendre function, and These are the coplanar latitude and longitude of the sphere, respectively. and These are the wavelet expansion coefficients for the mid-frequency and high-frequency bands. and The scale and location of wavelet decomposition. , and The highest order of spherical harmonics for each spectral band. , and The maximum scale parameter for waveform decomposition. For scale The number of wavelet basis functions at that point For calculation points and scales The first The spherical distance between the centers of each wavelet For scale The corresponding central spherical harmonicity, For bandwidth parameters, For Bouguer anomalies, This represents the average normal gravity.

[0059] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0060] 1. In this invention, a geoid construction method based on least squares spectral combination and spherical wavelet is proposed. First, multi-source gravity data related to the region to be constructed is acquired. Then, the spectrum of the multi-source gravity data is divided into low-frequency, mid-frequency, and high-frequency components. First, a spectral weighting function is constructed based on the modified spherical Butterworth function to accurately define the low-frequency, mid-frequency, and high-frequency bands. Then, using the determined spectral weights, the multi-source gravity observation data is weighted and modeled to establish low-frequency, mid-frequency, and high-frequency geoid component models respectively. First, the zero-order term and transformation parameters are obtained. Then, the optimal weights of the low-frequency, mid-frequency, and high-frequency components are combined with the zero-order term and transformation parameters to obtain the final geoid result. The advantages of this design are as follows:

[0061] First, by using a collaborative modeling framework of low-frequency, mid-frequency, and high-frequency spectrum domains, we can adapt to the physical characteristics of multi-source gravity data, avoid the problems of spectrum leakage, spectral discontinuity, and noise amplification caused by direct splicing, and ensure the frequency domain continuity and physical consistency of the geoid model.

[0062] Secondly, by using wavelet multi-scale decomposition to enhance local features in the spatial domain, the recovery capability of medium and short scales, especially complex terrain areas, is significantly improved, resulting in better synergy between spatial local scale and frequency domain.

[0063] Thirdly, by optimizing the weights of the least squares spectrum combination, the weights of each frequency band are estimated to minimize the residuals in a statistical sense, rather than by empirical weighting or rule splicing, thus minimizing the model residuals in a statistical sense and avoiding the accumulation of errors caused by empirical weighting.

[0064] Fourthly, this design is suitable for high-precision modeling of key areas with abundant data and complex terrain, and can also meet the basic accuracy requirements of areas with sparse data. Through spectral continuity, local scale adaptability and error controllability, it can improve the accuracy of geoid models and enhance the robustness of complex terrain areas.

[0065] Therefore, the model of this invention has good stability in complex terrain areas and when the recovery of medium and high frequency gravity signals is incomplete.

[0066] 2. In this invention, a geoid construction method based on least squares spectral combination and spherical wavelets, firstly, the frequency band weight function derived from least squares spectral combination is used to ensure that satellite, airborne, and ground data dominate within their respective dominant spherical harmonic orders, and that the weights naturally decay to near zero in non-dominant frequency bands, achieving segmented contribution at the frequency domain level. Secondly, spherical wavelet expansion possesses dual locality characteristics in both space and frequency, strictly confining wavelet functions of different scales within their corresponding frequency bands, effectively blocking spectral leakage and spatial aliasing. Through segmented modeling in the frequency domain, least squares spectral combination, and multi-scale wavelet decomposition, effective isolation of data from different frequency bands is achieved, mathematically avoiding spectral leakage and cross-interference. Therefore, this invention avoids spectral interference, resulting in more accurate results.

[0067] 3. This invention discloses a geoid construction method based on least-squares spectral combination and spherical wavelets. The long-wavelength portion of the geoid can be recovered using the GRACE / GOCE satellite gravity field model, and its accuracy is generally 10-30 cm when verified with GNSS / leveling data. By fusing satellite gravity (long-wavelength information) and ground gravity (high-frequency information), and verifying with GNSS / leveling data, the geoid modeling accuracy can reach 2-6 cm. Furthermore, introducing airborne gravity data, especially in areas with sparse ground observations or complex terrain, can effectively improve the recovery of mid-frequency gravity field signals. After verification with GNSS / leveling data, the geoid constructed by fusing multiple types of gravity data can achieve an accuracy of 1-3 cm. Therefore, this invention improves the accuracy of geoid modeling. Attached Figure Description

[0068] Figure 1 This is a flowchart of the present invention.

[0069] Figure 2 This is a flowchart of the wavelet multi-scale decomposition and least squares spectrum combination of the present invention. Detailed Implementation

[0070] The present invention will be further described in detail below with reference to the accompanying drawings and specific embodiments.

[0071] See Figures 1 to 2 A method for constructing a geoid based on least squares spectral combination and spherical wavelet, comprising the following steps:

[0072] The first step is to acquire multi-source gravity data for the region to be constructed, including global gravity field model data, airborne gravity data, and ground gravity data; then, based on the signal and error spectrum characteristics of each data source, their effective contributions are attributed to low-frequency, mid-frequency, and high-frequency components, respectively.

[0073] The second step is to first construct a spectral weighting function based on the modified spherical Butterworth function to accurately define the low-frequency, mid-frequency and high-frequency bands; then, using the determined spectral weights, perform weighted modeling on the multi-source gravity observation data to establish low-frequency, mid-frequency and high-frequency geoid component models respectively.

[0074] The third step is to first obtain the zero-order term and transformation parameters, and then combine the optimal weights of the low-frequency, mid-frequency, and high-frequency components with the zero-order term and transformation parameters to obtain the final geoid result.

[0075] After obtaining the final geoid result, the accuracy verification is also included.

[0076] The accuracy verification specifically involves using GNSS leveling data to check the accuracy of the constructed geoid model, with an accuracy of 1-2 cm in plain areas and 2-3 cm in mountainous or plateau areas.

[0077] The multi-source gravity data also includes terrain data and GNSS leveling data.

[0078] The low-frequency, mid-frequency, and high-frequency components are specifically:

[0079] The range of spherical harmonic orders corresponding to the low-frequency components is: 2≤n≤190;

[0080] The range of spherical harmonic orders corresponding to the intermediate frequency components is: 190 < n ≤ 545;

[0081] The range of spherical harmonic orders corresponding to high-frequency components is: n > 545;

[0082] The dominant data source for the low-frequency component is global gravity field model data, the dominant data source for the mid-frequency component is airborne gravity data, and the dominant data source for the high-frequency component is ground gravity data.

[0083] The low-frequency geoid component model is expressed in a spherical harmonic expansion, specifically as follows:

[0084] ;

[0085] in, The distance to the Earth's center. The gravitational constant is the constant of gravity. For reference to Earth's radius, and Let the order and degree of the spherical harmonic expansion be given. The maximum number of low-frequency components. and These are the fully normalized spherical harmonic coefficients. For a fully normalized association Legendre function, , These are spherical latitude and longitude, respectively.

[0086] The corresponding gravitational perturbation is obtained by radial differentiation of the perturbation potential:

[0087] ;

[0088] This section describes the long-wave gravity characteristics and serves as a benchmark for the enhancement of mid-to-high frequency components.

[0089] The mid-frequency geoid component modeling is specifically as follows: To fill the spectral gap between satellite and ground gravity measurements, airborne gravity data is used to extract mid-band disturbances. First, the low-frequency component and the topographic residual model are removed from the airborne gravity observations to obtain the mid-frequency component residual field.

[0090] ;

[0091] in, This represents the residual gravity perturbation in the mid-frequency component. For airborne gravity disturbance observation, For the middle arrive RTM correction;

[0092] To achieve optimal spectral localization, the mid-frequency component perturbation is represented using spherical wavelet expansion:

[0093] ;

[0094] in, To control the scale parameters of spatial resolution. For the position index of the wavelet center, For scale The number of wavelet basis functions corresponding to a location. The wavelet expansion coefficients are given, and the wavelet basis functions are defined as follows:

[0095] .

[0096] The high-frequency geoid component model is as follows: the high-frequency component reflects local density changes and is obtained by processing ground gravity observations. The high-frequency component residual gravity disturbance is expressed as follows:

[0097] ;

[0098] in, For ground gravity disturbance observation, For the middle to RTM correction;

[0099] High-frequency signals are represented using band-limited wavelet expansion:

[0100] ;

[0101] in, These are the wavelet expansion coefficients. For high-frequency wavelet basis functions:

[0102] .

[0103] In the second step, the determined spectral weights are specifically: to ensure that the sum of the spectral weights of the three frequency bands is always 1.

[0104] ;

[0105] in, , and These are the spectral weights for GGM, airborne and ground gravity disturbances, respectively.

[0106] Spectral weights are calculated using the following formula:

[0107] ;

[0108] in, and It is calculated based on an improved spherical Butterworth filter;

[0109] ;

[0110] in, .

[0111] In the third step, the zero-order term and transformation parameters are first obtained. Then, the optimal weights of the low-frequency, mid-frequency, and high-frequency components are combined with the zero-order term and transformation parameters to obtain the final geoid result. Specifically, the perturbation potential function can be expressed as a weighted superposition of multi-band spectral components.

[0112] ;

[0113] in, , and Represent the perturbation potential spectral components in the low, medium, and high frequency bands, respectively, and the weighting function... To achieve optimal spectral combination, a rationally designed spectral weighting function is used to ensure that the contributions of each frequency band are evenly distributed across different spatial scales.

[0114] The geoid height can be obtained from the disturbance potential:

[0115] ;

[0116] in, Calculated by GGM, Calculated from airborne gravity data, Calculated from land gravity data, The zero-degree term, The difference between the geoid and the quasi-geoid is expressed as:

[0117] ;

[0118] The wavelet basis function with spectral weights is defined as:

[0119] ;

[0120] in, and These are wavelet kernel functions for the frequency bands;

[0121] ;

[0122] in, The gravitational constant of the Earth's center of mass. for The gravitational constant of the reference ellipsoid, The average Earth radius, The zero-degree geoid term is used to explain the difference between the reference Earth model and the actual gravity field. To reference the normal gravity on the surface of the ellipsoid, The position of the geoid. For reference ellipsoid constant, and For the fully normalized spherical harmonic coefficients of the perturbation position, To fully normalize the associated Legendre function, and These are the coplanar latitude and longitude of the sphere, respectively. and These are the wavelet expansion coefficients for the mid-frequency and high-frequency bands. and The scale and location of wavelet decomposition. , and The highest order of spherical harmonics for each spectral band. , and The maximum scale parameter for waveform decomposition. For scale The number of wavelet basis functions at that point For calculation points and scales The first The spherical distance between the centers of each wavelet For scale The corresponding central spherical harmonicity, For bandwidth parameters, For Bouguer anomalies, This represents the average normal gravity.

[0123] The supplementary technical features of this invention are as follows:

[0124] Satellites, limited by their high orbital altitude, cannot observe mid-to-high frequencies and therefore cannot replace airborne or ground-based data. While airborne gravity data performs well in the mid-frequency range, high-frequency signals are significantly attenuated at flight altitudes, failing to achieve the resolution of ground-based gravity data. Furthermore, its regional coverage cannot provide a global longwave reference, thus it cannot replace satellite data either. Although ground-based data possesses high-frequency advantages, its point-like, discrete, and discontinuous coverage leads to aliasing and interpolation errors when directly used for mid-frequency modeling. It cannot provide the continuous mid-frequency support required by airborne data, nor can it provide a stable global longwave framework. Therefore, these three types of data are irreplaceable in their respective frequency bands and are essential for ensuring spectral continuity, consistency of multi-source information, and centimeter-level geoid modeling accuracy.

[0125] Example 1:

[0126] A method for constructing a geoid based on least squares spectral combination and spherical wavelet, comprising the following steps:

[0127] Step 1: Data Acquisition and Frequency Band Division: First, acquire multi-source gravity data for the region to be constructed, including global gravity field model data, airborne gravity data, and ground gravity data; then, based on the signal and error spectrum characteristics of each data source, assign their effective contributions to the low-frequency, mid-frequency, and high-frequency components respectively.

[0128] The second step, spectral weight determination and frequency band modeling: First, a spectral weight function is constructed based on the modified spherical Butterworth function to accurately define the low-frequency, mid-frequency and high-frequency bands; then, using the determined spectral weights, multi-source gravity observation data are weighted and modeled to establish low-frequency, mid-frequency and high-frequency geoid component models respectively.

[0129] The third step is geoid fusion: First, obtain the zero-order term and transformation parameters, and then combine the optimal weights of the low-frequency, mid-frequency, and high-frequency components with the zero-order term and transformation parameters to obtain the final geoid result.

[0130] Example 2:

[0131] Example 2 is basically the same as Example 1, except that:

[0132] A method for constructing a geoid based on least squares spectral combination and spherical wavelet, comprising the following steps:

[0133] Step 1: Data Acquisition and Frequency Band Division: First, acquire multi-source gravity data for the region to be constructed, including global gravity field model data, airborne gravity data, and ground gravity data; then, based on the signal and error spectrum characteristics of each data source, assign their effective contributions to the low-frequency, mid-frequency, and high-frequency components respectively.

[0134] The second step, spectral weight determination and frequency band modeling: First, a spectral weight function is constructed based on the modified spherical Butterworth function to accurately define the low-frequency, mid-frequency and high-frequency bands; then, using the determined spectral weights, multi-source gravity observation data are weighted and modeled to establish low-frequency, mid-frequency and high-frequency geoid component models respectively.

[0135] The third step is geoid fusion: First, obtain the zero-order term and transformation parameters, and then combine the optimal weights of the low-frequency, mid-frequency, and high-frequency components with the zero-order term and transformation parameters to obtain the final geoid result.

[0136] The fourth step is to use GNSS leveling data to check the accuracy of the constructed geoid model, ensuring that the accuracy reaches 1-2 cm in plain areas and 2-3 cm in mountainous or plateau areas.

[0137] Example 3:

[0138] Example 3 is basically the same as Example 1, except that:

[0139] The dominant data source for the low-frequency component is global gravity field model data, the dominant data source for the mid-frequency component is airborne gravity data, and the dominant data source for the high-frequency component is ground gravity data.

[0140] Satellite gravity, airborne gravity, and ground gravity differ fundamentally in observation altitude, spatial resolution, and noise propagation characteristics, thus naturally corresponding to different frequency bands of the gravity field. Satellite gravity, situated at altitudes of hundreds of kilometers, experiences an exponential decay in high-frequency components with altitude, only able to reliably recover long-wavelength information with spherical harmonic orders approximately n < 190. Airborne gravity, at altitudes of hundreds to thousands of meters, offers continuous regional coverage and moderate spatial resolution, effectively recovering mid-frequency gravity variations with n ≈ 190–545, effectively compensating for the spectral gap between satellites and the ground. Ground gravity, close to the Earth's surface, directly reflects local density and topographic changes, reliably providing fine short-wavelength structures with n > 545. Therefore, the optimal matching relationship between these three observation methods and the low-frequency, mid-frequency, and high-frequency segmentation of the gravity field is determined by their physical observability and is theoretically inevitable. The high-frequency band is limited to orders > 545.

[0141] Example 4:

[0142] Example 4 is basically the same as Example 1, except that:

[0143] The low-frequency component model is expressed in spherical harmonic expansion, specifically as follows:

[0144] ;

[0145] in, The distance to the Earth's center. The gravitational constant is the constant of gravity. For reference to Earth's radius, and Let the order and degree of the spherical harmonic expansion be given. The maximum number of low-frequency components. and These are the fully normalized spherical harmonic coefficients. For a fully normalized association Legendre function, , These are spherical latitude and longitude, respectively.

[0146] The corresponding gravitational perturbation is obtained by radial differentiation of the perturbation potential:

[0147] ;

[0148] This section describes the long-wave gravity characteristics and serves as a benchmark for the enhancement of mid-to-high frequency components.

[0149] The modeling of the mid-frequency component specifically involves: to fill the spectral gap between satellite and ground gravity measurements, airborne gravity data is used to extract mid-band perturbations. First, the low-frequency component and the terrain residual model are removed from the airborne gravity observations to obtain the mid-frequency component residual field.

[0150] ;

[0151] in, This represents the residual gravity perturbation in the mid-frequency component. For airborne gravity disturbance observation, For the middle arrive RTM correction;

[0152] To achieve optimal spectral localization, the mid-frequency component perturbation is represented using spherical wavelet expansion:

[0153] ;

[0154] in, To control the scale parameters of spatial resolution. For the position index of the wavelet center, For scale The number of wavelet basis functions corresponding to a location. The wavelet expansion coefficients are given, and the wavelet basis functions are defined as follows:

[0155] .

[0156] The high-frequency components are modeled, reflecting local density changes. These components are obtained by processing ground gravity observations, and the residual gravity perturbation of the high-frequency components is expressed as follows:

[0157] ;

[0158] in, For ground gravity disturbance observation, For the middle to RTM correction;

[0159] High-frequency signals are represented using band-limited wavelet expansion:

[0160] ;

[0161] in, These are the wavelet expansion coefficients. For high-frequency wavelet basis functions:

[0162] .

[0163] The low-frequency satellite model suppresses high-order noise using an MSB filter, ensuring that only robust long-wave components are retained. Airborne data is corrected using a bandpass wavelet kernel combined with residual terrain effect (RTM). RTM parameters (such as the integration radius) have a slight impact on accuracy; the optimal parameters were determined experimentally, with an integration radius of 150 km selected to reduce strip noise and attitude errors, highlighting the effective mid-frequency signal. Ground data utilizes localized high-frequency wavelets and high-frequency residuals to eliminate trend terms and mid-frequency residues, outputting only a reliable short-wave structure. Combined with differentiated denoising methods across the three frequency bands, the fusion process ultimately achieves noise-controlled, energy-clear, multi-source gravity information integration across frequency bands, thus avoiding mutual interference between different frequency bands.

[0164] Example 5:

[0165] Example 5 is basically the same as Example 1, except that:

[0166] In the second step, the frequency bands are normalized at each spherical harmonic order to generate the optimal weights for each frequency band at different spherical harmonic orders. Specifically, the sum of the spectral weights of the three frequency bands is kept constant at 1.

[0167] ;

[0168] in, , and These are the spectral weights for GGM, airborne and ground gravity disturbances, respectively.

[0169] Spectral weights are calculated using the following formula:

[0170] ;

[0171] in, and It is calculated based on an improved spherical Butterworth filter;

[0172] ;

[0173] in, .

[0174] and The values ​​are derived from the dominant frequency bands and noise characteristics of different data sources in the spherical harmonic domain, among which... Corresponding to the transition from low frequency to mid frequency, Corresponding to the transition from intermediate frequency to high frequency, they are generated by a modified spherical Butterworth filter to form a smooth, monotonic frequency domain response, ensuring a continuous and ringless transition region.

[0175] In the least squares spectral combination method, the initial spectral response of each data source (satellite, airborne, ground) is first generated by the MSB-type filter function according to its noise characteristics and effective frequency band. Then, the final three-band spectral weight is constructed by normalization at each spherical harmonic order. That is, the original filter response of each frequency band is used as the numerator and the sum of the three is used as the denominator, thus forming a result that satisfies Pnlow + Pnmed + Pnhigh = 1. This normalization strategy ensures that regardless of the spherical harmonic order, all three types of data contribute within their dominant frequency bands, collectively forming a continuous and stable spectral energy allocation. To avoid situations where a certain frequency band has excessively low or even near-zero weights after normalization, this method ensures the smoothness and stability of the weights through three mechanisms: First, the MSB filter has smooth transition characteristics, making the original spectral weights monotonically change with the spherical harmonic order without drastic jumps; second, the method ensures that each frequency band contains effective energy through the stepwise construction of low-frequency, mid-frequency, and high-frequency residuals, preventing the "signal itself from approaching zero" situation; third, the spectral weights are obtained based on a least-squares noise model, and the asymptotic characteristics of noise in different frequency bands ensure that the weights naturally exhibit a smooth change. Considering these factors, the final weights of the three frequency bands remain continuous and without abrupt changes throughout the entire frequency domain, avoiding the problem of extreme near-zero weights.

[0176] Example 6:

[0177] Example 6 is basically the same as Example 1, except that:

[0178] The perturbation potential function can be expressed as a weighted superposition of multi-band spectral components:

[0179] ;

[0180] in, , and Represent the perturbation potential spectral components in the low, medium, and high frequency bands, respectively, and the weighting function... To achieve optimal spectral combination, a rationally designed spectral weighting function is used to ensure that the contributions of each frequency band are evenly distributed across different spatial scales.

[0181] The geoid height can be obtained from the disturbance potential:

[0182] ;

[0183] in, Calculated by GGM, Calculated from airborne gravity data, Calculated from land gravity data, The zero-degree term, The difference between the geoid and the quasi-geoid is expressed as:

[0184] ;

[0185] The wavelet basis function with spectral weights is defined as:

[0186] ;

[0187] in, and These are wavelet kernel functions for the respective frequency bands.

[0188] The fusion function is as follows:

[0189] ;

[0190] in, The gravitational constant of the Earth's center of mass. for The gravitational constant of the reference ellipsoid, The average Earth radius, The zero-degree geoid term is used to explain the difference between the reference Earth model and the actual gravity field. To reference the normal gravity on the surface of the ellipsoid, The position of the geoid. For reference ellipsoid constant, and For the fully normalized spherical harmonic coefficients of the perturbation position, To fully normalize the associated Legendre function, and These are the coplanar latitude and longitude of the sphere, respectively. and These are the wavelet expansion coefficients for the mid-frequency and high-frequency bands. and The scale and location of wavelet decomposition. , and The highest order of spherical harmonics for each spectral band. , and The maximum scale parameter for waveform decomposition. For scale The number of wavelet basis functions at that point For calculation points and scales The first The spherical distance between the centers of each wavelet For scale The corresponding central spherical harmonicity, For bandwidth parameters, For Bouguer anomalies, To average normal gravity, a small regularization term and threshold constraint can be added during the normalization process to ensure that each frequency band maintains a non-zero contribution and transitions smoothly in its dominant range, thereby forming a stable, continuous and physically reasonable weight allocation.

[0191] The wavelet coefficients are solved using the least squares method, and the functional model of the residual gravity observation value l can be expressed as:

[0192] ;

[0193] in, For the residual vector, To observe the residual gravity perturbation vector;

[0194] ;

[0195] A is the design matrix, which consists of the values ​​of the wavelet basis functions at the observation points and is used to map the wavelet coefficients to gravity perturbations.

[0196] For the unknown wavelet coefficient vector: ;

[0197] The least squares objective function is: ;

[0198] To enhance the stability of the solution and avoid overfitting, especially when the design matrix is ​​ill-conditioned or the data distribution is uneven, a regularization term can be introduced:

[0199] ;

[0200] in, This is a regularization parameter, which can be determined through cross-validation or the L-curve method;

[0201] The corresponding normal equation is: ;

[0202] Solving this system of equations yields the estimated wavelet coefficients:

[0203] ;

[0204] Calculate the coefficients and Then, by substituting it into the spectral weight calculation formula, the mid-frequency and high-frequency geoid elevation components can be calculated.

[0205] The above description is only a preferred embodiment of the present invention. The scope of protection of the present invention is not limited to the above embodiments. Any equivalent modifications or changes made by those skilled in the art based on the content disclosed in the present invention should be included within the scope of protection set forth in the claims.

Claims

1. A method for constructing a geoid based on least squares spectral combination and spherical wavelet, characterized in that: The geoid construction method based on least squares spectral combination and spherical wavelet includes the following steps: The first step is to acquire multi-source gravity data for the region to be constructed, including global gravity field model data, airborne gravity data, and ground gravity data. Then, based on the signal and error spectrum characteristics of each data source, their effective contributions are attributed to low-frequency, mid-frequency, and high-frequency components, respectively. The second step is to first construct a spectral weighting function based on the improved spherical Butterworth function form to accurately define the low-frequency, mid-frequency and high-frequency bands. Then, using the determined spectral weights, the multi-source gravity observation data is weighted and modeled to establish geoid component models for low-frequency, mid-frequency and high-frequency components respectively. The third step is to first obtain the zero-order term and transformation parameters, and then combine the optimal weights of the low-frequency, mid-frequency, and high-frequency components with the zero-order term and transformation parameters to obtain the final geoid result. Specifically, the perturbation potential function is represented as a weighted superposition of multi-band spectral components: ； in, , and These represent the perturbation bits in the low, mid, and high frequency bands, respectively, and the weighting functions are... To achieve optimal spectral combination, a rationally designed spectral weighting function is used to ensure that the contributions of each frequency band are evenly distributed across different spatial scales. The geoid height can be obtained from the disturbance potential: ； in, This is a low-frequency elevation anomaly, calculated using a global gravity field model. This is a mid-frequency elevation anomaly, calculated from airborne gravity data. This is a high-frequency elevation anomaly, calculated from land gravity data. For zero-order terms, The difference between the geoid and the quasi-geoid is expressed as: ； The wavelet basis function with spectral weights is defined as: ； in, and Wavelet kernel functions for the mid- and high-frequency bands are respectively: ； in, The distance to the Earth's center. The gravitational constant of the Earth's center of mass. for The gravitational constant of the reference ellipsoid, The average Earth radius, To reference the normal gravity on the surface of the ellipsoid, The potential constant of the geoid. For the potential constant of the reference ellipsoid, and For the fully normalized spherical harmonic coefficients of the perturbation position, To fully normalize the associated Legendre function, and These are the coplanar latitude and longitude of the sphere, respectively. and These are the wavelet expansion coefficients for the mid-frequency and high-frequency bands. and The scale and location of wavelet decomposition. and Let the order and degree of the spherical harmonic expansion be given. , and The highest order of the spherical harmonics in each spectral band. , and The maximum scale parameter for waveform decomposition. For scale The number of wavelet basis functions at that point For calculation points and scales The first The spherical distance between the centers of each wavelet For scale The corresponding central spherical harmonicity, For bandwidth parameters, For Bouguer anomalies, This represents the average normal gravity.

2. The geoid construction method based on least squares spectral combination and spherical wavelet as described in claim 1, characterized in that: After obtaining the final geoid result, the accuracy verification is also included.

3. The geoid construction method based on least squares spectral combination and spherical wavelet according to claim 2, characterized in that: The accuracy verification specifically involves using GNSS leveling data to check the accuracy of the constructed geoid model, with an accuracy of 1-2 cm in plain areas and 2-3 cm in mountainous or plateau areas.

4. The geoid construction method based on least squares spectral combination and spherical wavelet as described in claim 1, characterized in that: The multi-source gravity data also includes terrain data and GNSS leveling data.

5. The geoid construction method based on least squares spectral combination and spherical wavelet according to claim 1, characterized in that: The low-frequency, mid-frequency, and high-frequency components are specifically: The range of spherical harmonic orders corresponding to the low-frequency components is: 2≤n≤190; The range of spherical harmonic orders corresponding to the intermediate frequency components is: 190 < n ≤ 545; The range of spherical harmonic orders corresponding to high-frequency components is: n > 545; The dominant data source for the low-frequency component is global gravity field model data, the dominant data source for the mid-frequency component is airborne gravity data, and the dominant data source for the high-frequency component is ground gravity data.

6. The geoid construction method based on least squares spectral combination and spherical wavelet according to claim 5, characterized in that: The low-frequency geoid component model is expressed in a spherical harmonic expansion, specifically as follows: ； in, The maximum number of low-frequency components. and These are the fully normalized spherical harmonic coefficients. For a fully normalized association Legendre function, , These are the coplanar latitude and longitude of the sphere, respectively. The corresponding gravitational perturbation is obtained by radial differentiation of the perturbation potential: ； This section describes the long-wave gravity characteristics and serves as a benchmark for the enhancement of mid-to-high frequency components.

7. The geoid construction method based on least squares spectral combination and spherical wavelet according to claim 6, characterized in that: The modeling of the mid-frequency geoid component is specifically as follows: To fill the spectral gap between satellite and ground gravity measurements, airborne gravity data is used to extract mid-band disturbances. First, the low-frequency component and the topographic residual model are removed from the airborne gravity observations to obtain the mid-frequency component residual field. ； in, This represents the residual gravity perturbation in the mid-frequency component. For airborne gravity disturbance observation, For the middle arrive Residual terrain model correction; To achieve optimal spectral localization, the mid-frequency component perturbation is represented using spherical wavelet expansion: ； in, To control the scale parameters of spatial resolution. For the position index of the wavelet center, For scale The number of wavelet basis functions corresponding to a location. These are the wavelet expansion coefficients. This represents the maximum scale for mid-frequency waveform decomposition. For the maximum scale of low-frequency waveform decomposition, the wavelet basis function is defined as: ； in, This represents the maximum spherical harmonic order in the mid-frequency range. For the intermediate frequency spectrum weighting function, For mid-frequency wavelet kernel functions, This is used to calculate the spherical distance between a point and the k-th wavelet center at scale j.

8. The geoid construction method based on least squares spectral combination and spherical wavelet according to claim 7, characterized in that: The high-frequency geoid component model is as follows: the high-frequency component reflects local density changes and is obtained by processing ground gravity observations. The high-frequency component residual gravity disturbance is expressed as follows: ； in, For ground gravity disturbance observation, For the middle to Residual terrain model correction; High-frequency signals are represented using band-limited wavelet expansion: ； in, These are the wavelet expansion coefficients. This represents the maximum scale for high-frequency waveform decomposition. For high-frequency wavelet basis functions: ； in, This represents the maximum spherical harmonic order in the high-frequency band. For high-frequency spectrum weighting functions, This is a high-frequency wavelet kernel function.

9. The geoid construction method based on least squares spectral combination and spherical wavelet according to claim 8, characterized in that: In the second step, the determined spectral weights are specifically: to ensure that the sum of the spectral weights of the three frequency bands is always 1. ； in, , and These are the spectral weighting functions for the low, medium, and high frequency bands, respectively. The spectral weighting function is calculated using the following formula: ； in, and It is calculated based on the improved spherical Butterworth function; ； in, , This represents the cutoff order of the frequency band characteristics.

Images