A gravimeter drift correction geodetic surveying method, device, medium and product based on drift rate function

Abstract

The application provides a drift rate function-based drift correction gravimeter drift geodetic survey method, device, medium and product, relates to the geophysical exploration technical field, and the method comprises the steps of setting a plurality of measuring points on a preset measuring line, using a gravimeter to carry out outgoing and return gravity measurement, and obtaining outgoing gravity values, return gravity values and corresponding measurement time of each measuring point; determining target instantaneous drift rates of each measuring point according to the time difference value of the outgoing and return gravity values, establishing a drift rate function with the measurement time interval as the independent variable and the target instantaneous drift rate as the dependent variable; performing definite integral operation on the drift rate function to obtain the target cumulative drift amount of each measuring point; combining the measurement time proportion of the outgoing and return to determine a time weight coefficient, multiplying the target cumulative drift amount by the time weight coefficient to obtain a time drift correction amount; and using the time drift correction amount to correct the drift of the outgoing and return gravity values, and obtaining target corrected gravity values.

Description

Technical Field

[0001] This application relates to the field of geophysical exploration technology, and in particular to a gravimeter drift correction geodetic mapping method, equipment, medium and product based on a drift rate function. Background Technology

[0002] With the rapid development of geophysical exploration technology and the continuous improvement of geodetic accuracy requirements, gravity measurement has become an important technical means for geological structure research, mineral resource exploration, and geodynamic monitoring.

[0003] In related technologies, a correction method based on the linear time drift assumption is typically used. In practice, surveyors first use a gravimeter to perform an initial gravity measurement at the starting point of the survey line, recording the initial gravity value and the measurement time. Then, they sequentially perform single-pass gravity measurements at each measuring point along the survey line, obtaining the observed gravity value at each point. Finally, they return to the starting point of the survey line to repeat the measurement. Based on the difference in gravity values ​​between the two measurements at the starting point and the time interval, they calculate the average drift rate. Then, using linear interpolation, they calculate the drift correction value for each measuring point according to the time difference between the measurement time and the starting time. While this method can achieve basic drift correction, the entire process heavily relies on the linear drift assumption.

[0004] However, the drift correction method described above, which estimates the uniform drift rate of the entire measurement line based solely on the two measurements at the starting point, is based on the ideal assumption that the instrument drift strictly follows a linear law. In reality, the drift of a gravimeter often exhibits nonlinear characteristics due to the influence of temperature changes and material fatigue on the quartz spring system. This drift correction method ignores the differences in the actual drift conditions at each measurement point, resulting in an inherent deviation between the linear model and the actual drift curve. More importantly, this deviation accumulates systematically over time, making the correction error larger for measurement points further away from the starting point, thus leading to lower correction accuracy for gravimeter drift measurements in related technologies. Summary of the Invention

[0005] This application provides a gravimeter drift correction geodetic mapping method, equipment, medium, and product based on a drift rate function, which is used to improve the correction accuracy of gravimeter drift measurements.

[0006] In a first aspect, this application provides a gravimeter drift correction geodetic mapping method based on a drift rate function, applied to the aforementioned electronic device. The method includes: setting multiple measuring points on a preset survey line; performing outward gravity measurements on each measuring point using a gravimeter in the outward sequence to obtain the outward gravity value and corresponding outward measurement time for each measuring point; and performing return gravity measurements on each measuring point using a gravimeter in the reverse order of the outward sequence to obtain the return gravity value and corresponding return measurement time for each measuring point; determining the target instantaneous drift rate of each measuring point based on the return gravity value, outward gravity value, outward measurement time, and return measurement time, using the time interval between the outward measurement time and the outward measurement start time as the independent variable, and the target instantaneous drift rate as the independent variable. With the drift rate as the dependent variable, a drift rate function is established. Starting from the measurement start time, the drift rate function is integrated along the time progression direction to obtain the target cumulative drift amount at each measuring point. The target cumulative drift amount represents the integral value of the drift amount at the measuring point from the measurement start time to the current measurement time along the measurement time axis. Based on the proportion of the outbound measurement time in the total outbound measurement time and the proportion of the return measurement time in the total return measurement time, the time weight coefficient of each measuring point is determined. The time drift correction amount of each measuring point is determined based on the target cumulative drift amount and the time weight coefficient. The time drift correction amount is used to correct the drift of the outbound gravity value and the return gravity value to obtain the target corrected gravity value.

[0007] By adopting the above technical solution, the outbound and return gravity measurements are performed in reverse order along the same measuring line, ensuring a symmetrical temporal distribution of the outbound and return measurement times at each measuring point. This symmetry provides a reliable time reference for subsequently determining the instantaneous drift rate of the target. A drift rate function is established based on the difference between the return and outbound gravity values ​​and the corresponding time relationships of the measurement times. This function accurately reflects the drift variation of the gravimeter throughout the measurement process. The cumulative drift is obtained through definite integral calculation starting from the measurement initiation time, enabling continuous cumulative calculation of the drift and avoiding discontinuities that may arise from discrete point correction. The time weighting coefficient, combined with the ratio of outbound and return measurement times, allows for a reasonable allocation of drift correction effects. Finally, the corrected target gravity value is obtained through comprehensive correction of the outbound and return gravity values, achieving precise compensation for drift errors. This solves the technical problem of low correction accuracy for gravimeter drift measurement in related technologies, achieving a significant improvement in the correction accuracy of gravimeter drift measurement.

[0008] Optionally, the target instantaneous drift rate of each measuring point is determined based on the return gravity value, the outbound gravity value, the outbound measurement time, and the return measurement time. Specifically, this includes: subtracting the outbound gravity value from the return gravity value to obtain the gravity difference at each measuring point; calculating the first difference between the return measurement time and the outbound measurement time to obtain the round-trip time interval at each measuring point; dividing the gravity difference by the round-trip time interval to obtain the average drift rate at each measuring point; obtaining the movement time of the gravimeter between every two adjacent measuring points, and obtaining the dwell observation time of the gravimeter at each measuring point; and determining the target instantaneous drift rate based on the movement time and dwell observation time. The actual observation period of each measuring point is determined; the average drift rate of each measuring point is correlated with the midpoint of the corresponding actual observation period to construct a discrete drift rate dataset; interpolation calculation is performed on the discrete drift rate dataset to obtain a continuous drift rate time function; the function value of the midpoint of the actual observation period of each measuring point is extracted from the drift rate time function, and the function value is used as the initial instantaneous drift rate of each measuring point; historical drift characteristic data of the gravimeter is obtained, and drift correction coefficients are determined based on the historical drift characteristic data; the initial instantaneous drift rate is corrected using the drift correction coefficients to obtain the target instantaneous drift rate.

[0009] By employing the above technical solution, the difference between the return gravity value and the outbound gravity value, divided by the round-trip time interval, yields the average drift rate, providing a basic drift estimate for each measuring point. The combination of travel time and dwell time determines the actual observation period, allowing the average drift rate to be accurately correlated to the midpoint of the corresponding observation period, forming a discrete drift rate dataset. Interpolation converts the discrete data into a continuous drift rate time function, eliminating abrupt changes in drift rate between measuring points and ensuring the continuity of drift changes. The initial instantaneous drift rate extracted from the continuous function reflects the real-time drift characteristics of each measuring point, while the drift correction coefficient determined by historical drift characteristic data corrects the initial instantaneous drift rate. This fully considers the influence of the gravimeter's historical operating state on the current drift, making the target instantaneous drift rate more consistent with the actual drift characteristics of the gravimeter and improving the accuracy of drift rate estimation.

[0010] Optionally, the drift correction coefficient is determined based on historical drift characteristic data, specifically including: extracting the drift rate variation curves of the gravimeter under different working durations from the historical drift characteristic data; determining the drift acceleration characteristic value of the gravimeter based on the drift rate variation curves; obtaining the cumulative working duration of the current measurement task, and determining the drift rate value corresponding to the cumulative working duration from the drift rate variation curves; comparing the drift rate value with a preset standard drift rate to obtain the drift rate ratio; and determining the drift correction coefficient K based on the drift acceleration characteristic value and the drift rate ratio using the following formula:

[0011] K = 1 + (R - 1) × (1 + α × T / T0)

[0012] Where K is the drift correction coefficient, α is the drift acceleration characteristic value, R is the drift rate ratio, T is the cumulative working time, and T0 is the preset standard working time.

[0013] By employing the above technical solution, the drift rate variation curve extracted from historical drift characteristic data can reflect the drift evolution law of the gravimeter under different working durations, and the drift acceleration characteristic value can quantify the changing trend of the drift rate. The ratio of the drift rate value corresponding to the cumulative working duration to the preset standard drift rate can reflect the degree of deviation of the current drift state from the standard state. In the drift correction coefficient formula, the drift rate ratio R can reflect the basic deviation level, and the product of the drift acceleration characteristic value α and the time ratio T / T0 can reflect the time-varying characteristics of the drift. The two interact through the form (R-1)×(1+α×T / T0), enabling the drift correction coefficient K to dynamically adapt to the drift characteristic changes of the gravimeter in different working stages. This correction method, which comprehensively considers both basic deviation and time-varying characteristics, can ensure that the drift correction coefficient accurately reflects the actual working state of the gravimeter, improving the adaptability and accuracy of drift correction.

[0014] Optionally, a drift rate function is established using the time interval between the outward measurement time and the outward measurement start time as the independent variable and the target instantaneous drift rate as the dependent variable. Specifically, this includes: calculating a second difference between the outward measurement time and the outward measurement start time to obtain the outward measurement time interval for each measuring point; determining the instantaneous drift rate difference between every two adjacent measuring points and the outward measurement time interval difference between every two adjacent measuring points; determining the drift acceleration based on the instantaneous drift rate difference and the outward measurement time interval difference; establishing an initial drift rate function using a preset numerical fitting method with the outward measurement time interval as the independent variable and the target instantaneous drift rate as the dependent variable; acquiring the gravimeter's operating state parameters and determining the gravimeter's abnormal drift periods based on these parameters; extracting the gravimeter's drift characteristics from the abnormal drift periods and adjusting the weights of the abnormal function portion within the abnormal drift periods in the initial drift rate function based on these drift characteristics to obtain the drift rate function.

[0015] By adopting the above technical solution, the outbound measurement time interval can be used as the independent variable to establish a unified time benchmark. The ratio of the instantaneous drift rate difference to the time interval difference can determine the drift acceleration, reflecting the rate of change of the drift rate. The preset numerical fitting method establishes an initial drift rate function based on the time interval and the target instantaneous drift rate, which can realize the conversion of discrete data into a continuous function. The drift anomaly periods identified by the working state parameters can mark the abnormal working period of the gravimeter, and the drift features extracted from the anomaly periods can quantify the degree of anomaly. By adjusting the weight of the anomaly function part in the initial drift rate function, the influence of abnormal data on the overall fitting effect can be reduced, so that the final drift rate function can more accurately reflect the normal drift law of the gravimeter, while not completely ignoring the information of the anomaly periods, achieving a balance between normal drift and abnormal drift, and improving the robustness of the function.

[0016] Optionally, starting from the measurement start time, a first definite integral operation is performed on the drift rate function along the time progression direction to obtain the target cumulative drift amount at each measurement point. Specifically, this includes: determining the first derivative of the drift rate function and comparing the absolute value of the first derivative with a preset gradient threshold to obtain a gradient comparison result; dividing the drift rate function into a stable period and a rapidly changing period based on the gradient comparison result; performing a first numerical integration calculation on the stable function part of the drift rate function located within the stable period using a first integration step size to obtain the stable period integral value; performing a second numerical integration calculation on the rapidly changing function part of the drift rate function located within the rapidly changing period using a second integration step size to obtain the rapidly changing period integral value, where the second integration step size is smaller than the first integration step size; and accumulating the stable period integral value and the rapidly changing period integral value according to the time sequence of each measurement point on the time axis to obtain the target cumulative drift amount.

[0017] By employing the above technical solution, the absolute value of the first derivative of the drift rate function can reflect the drastic change in the drift rate, and comparison with a preset gradient threshold can achieve a quantitative assessment of the function's changing characteristics. Based on the gradient comparison results, the function is divided into stationary periods and rapidly changing periods, allowing for the use of appropriate integration strategies for function parts with different changing characteristics. The first integration step size ensures computational efficiency for the stationary function part, while the second integration step size ensures integration accuracy for the rapidly changing function part through more refined integration calculations. The differentiated settings of the two step sizes achieve a balance between computational efficiency and accuracy. The integration values ​​for the stationary and rapidly changing periods are accumulated sequentially over time, ensuring the temporal continuity of the accumulated drift. This allows the target accumulated drift to accurately reflect the complete drift accumulation process from the start of the measurement to the current moment, achieving efficient integration calculation with adaptive accuracy.

[0018] Optionally, based on the proportion of outbound measurement time in the total outbound measurement time and the proportion of return measurement time in the total return measurement time, a time weighting coefficient for each measurement point is determined. Then, based on the target cumulative drift and the time weighting coefficient, a time drift correction amount for each measurement point is determined. Specifically, this includes: obtaining the outbound measurement end time and performing a third difference calculation between the outbound measurement end time and the outbound measurement start time to obtain the total outbound measurement time; obtaining the return measurement start time and return measurement end time and performing a fourth difference calculation between the return measurement end time and the return measurement start time to obtain the total return measurement time; and dividing the outbound measurement time interval by the outbound measurement time interval. The total measurement time is used to obtain the outbound measurement time ratio; the fifth difference between the end time and the return measurement time is calculated to obtain the return measurement time interval; the return measurement time interval is divided by the total return measurement time to obtain the return measurement time ratio; the cumulative target drift is used as the theoretical drift correction; the time weighting coefficient is determined based on the outbound and return measurement time ratios; the outbound and return theoretical corrections are determined based on the theoretical drift correction and the time weighting coefficient; the round-trip gravity values ​​are compared at multiple measurement points to obtain the round-trip gravity closure difference; the outbound and return theoretical corrections are corrected based on the round-trip gravity closure difference to obtain the time drift correction.

[0019] By adopting the above technical solution, the total time for outbound and return measurements can provide a benchmark for time normalization. The ratio of the outbound measurement time interval to the total time is the outbound measurement time ratio, and the ratio of the return measurement time interval to the total time is the return measurement time ratio. These two ratios can jointly determine the basis for calculating the time weighting coefficient. The target cumulative drift, as a theoretical drift correction, can reflect the theoretical drift influence based on the drift rate function model. The round-trip gravity closure error can quantify the inconsistency between round-trip measurements in actual measurements. The time weighting coefficient determines the weight allocation for outbound and return trips based on the relative positions of the measuring points in the measurement time sequence. The outbound and return theoretical corrections calculated based on the theoretical drift correction and the time weighting coefficient reflect the theoretical model's prediction of drift, while the round-trip gravity closure error reflects the deviation between the theoretical model and the actual measurement. By correcting the theoretical correction amounts for the outward and return journeys based on the round-trip gravity closure difference, the theoretical model predictions can be combined with actual measurement feedback. This ensures both the theoretical continuity and time consistency of drift correction, while also eliminating systematic biases in the theoretical model through feedback correction of the measured closure difference. This achieves synergistic optimization of theoretical and measured corrections, avoids excessive correction amounts caused by repeated corrections, and improves the accuracy and reliability of drift correction.

[0020] Optionally, the outbound and return gravity values ​​are corrected using a time drift correction to obtain the target corrected gravity value. Specifically, the method further includes: subtracting the outbound time drift correction from the outbound gravity value to obtain the outbound temporary corrected gravity value; and subtracting the return time drift correction from the return gravity value to obtain the return temporary corrected gravity value. Using the outbound and return temporary corrected gravity values ​​as iterative variables, the following iterative optimization operation is performed until the convergence characteristic value of the round-trip consistency deviation meets the preset convergence condition: calculating the sixth difference between the outbound and return temporary corrected gravity values ​​at the same measuring point to obtain the round-trip consistency deviation; determining the convergence characteristic value of the round-trip consistency deviation according to the preset convergence index; when... When the convergence eigenvalue does not meet the preset convergence condition, the fitting parameters of the drift rate function are adjusted according to the round-trip consistency deviation to obtain the adjusted drift rate function; the adjusted drift rate function is subjected to a second definite integral operation to obtain a new target cumulative drift amount; a new outbound time drift correction amount and a new return time drift correction amount are determined according to the new target cumulative drift amount and the time weight coefficient; the outbound gravity value is drift corrected using the new outbound time drift correction amount to obtain a new outbound temporary correction gravity value, and the return gravity value is drift corrected using the new return time drift correction amount to obtain a new return temporary correction gravity value; the outbound temporary correction gravity value and the return temporary correction gravity value that meet the preset convergence condition are weighted and averaged according to the time weight coefficient to obtain the target correction gravity value.

[0021] By employing the above technical solution, the difference between the outbound and return temporary correction gravity values ​​forms the round-trip consistency deviation, which directly reflects the quality of the current correction effect. A preset convergence index converts the deviation into a convergence characteristic value, providing a quantitative criterion for iteration termination. When the convergence condition is not met, the round-trip consistency deviation guides the adjustment direction and magnitude of the drift rate function fitting parameters. The adjusted drift rate function generates a new target cumulative drift amount through a second definite integral operation, which in turn updates the time drift correction amount. The new correction amount recalibrates the original gravity value, forming a new temporary correction value and constituting a complete iterative closed loop. This iterative optimization mechanism based on deviation feedback allows the correction process to adaptively adjust until the round-trip consistency meets the preset requirements. Finally, by using a time weighting coefficient to perform a weighted average of the round-trip correction values ​​that meet the convergence condition, the optimal target correction gravity value can be obtained by fully utilizing the round-trip measurement information.

[0022] In a second aspect, embodiments of this application provide an electronic device comprising: one or more processors and a memory; the memory is coupled to one or more processors and is used to store computer program code, the computer program code including computer instructions, wherein one or more processors invoke the computer instructions to cause the electronic device to perform the method described in the first aspect and any possible implementation thereof.

[0023] Thirdly, embodiments of this application provide a computer program product containing instructions that, when the computer program product is run on an electronic device, cause the electronic device to perform the method described in the first aspect and any possible implementation thereof.

[0024] Fourthly, embodiments of this application provide a computer-readable storage medium including instructions that, when executed on an electronic device, cause the electronic device to perform the method described in the first aspect and any possible implementation thereof. Attached Figure Description

[0025] Figure 1 This is a flowchart illustrating a gravimeter drift correction geodetic mapping method based on a drift rate function in an embodiment of this application.

[0026] Figure 2 This is a schematic diagram of the physical device structure of an electronic device in an embodiment of this application. Detailed Implementation

[0027] The terminology used in the following embodiments of this application is for the purpose of describing particular embodiments only and is not intended to be limiting of this application. As used in the specification and appended claims of this application, the singular expressions “a,” “an,” “the,” “the,” “the,” and “this” are intended to include the plural expressions as well, unless the context clearly indicates otherwise. It should also be understood that the term “and / or” as used in this application refers to any or all possible combinations including one or more of the listed items.

[0028] Hereinafter, the terms "first" and "second" are used for descriptive purposes only and should not be construed as implying or suggesting relative importance or implicitly indicating the number of indicated technical features. Thus, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature, and in the description of the embodiments of this application, unless otherwise stated, "multiple" means two or more.

[0029] This application provides a geodetic mapping method for gravimeter drift correction based on a drift rate function, see reference. Figure 1 , Figure 1This is a flowchart illustrating a gravimeter drift correction geodetic mapping method based on a drift rate function, as described in this application, including the following steps:

[0030] Step S101: Set up multiple measuring points on the preset measuring line, and use a gravimeter to measure the outgoing gravity of each measuring point in the outgoing order to obtain the outgoing gravity value and the corresponding outgoing measurement time of each measuring point. Also, use a gravimeter to measure the return gravity of each measuring point in the return order, which is the opposite of the outgoing order, to obtain the return gravity value and the corresponding return measurement time of each measuring point.

[0031] Step S102: Determine the target instantaneous drift rate of each measuring point based on the return gravity value, the outgoing gravity value, the outgoing measurement time, and the return measurement time. Establish a drift rate function with the time interval between the outgoing measurement time and the outgoing measurement start time as the independent variable and the target instantaneous drift rate as the dependent variable.

[0032] Step S103: Starting from the measurement start time, perform the first definite integral operation on the drift rate function along the time progression direction to obtain the target cumulative drift amount of each measuring point. The target cumulative drift amount represents the integral value of the drift amount of the measuring point from the measurement start time to the current measurement time along the measurement time axis.

[0033] Step S104: Determine the time weight coefficient of each measuring point based on the proportion of outbound measurement time in the total outbound measurement time and the proportion of return measurement time in the total return measurement time, and determine the time drift correction amount of each measuring point based on the target cumulative drift amount and the time weight coefficient.

[0034] Step S105: Use time drift correction to correct the outbound and return gravity values ​​to obtain the target corrected gravity value.

[0035] In the above embodiment, a CG-6 relative gravimeter was used to measure gravity in a mining area. This gravimeter has a reading resolution of 0.1 mGal (milligal, 1 mGal = 10^-5 m / s²). A pre-set survey line was laid out along the main exploration line of the mining area, with a total length of 12 kilometers and 25 measuring points, spaced approximately 500 meters apart. Each measuring point was marked with a concrete pier and precisely located using GNSS (Global Navigation Satellite System). The outbound measurement began at 8:00 AM on March 15th, with the gravimeter starting from the secondary gravity baseline A and sequentially reaching measuring points P1, P2...P25. At each measuring point, the gravimeter automatically performed continuous observation for 120 seconds, recording 60 readings. After automatically removing outliers, the average was taken to obtain the relative gravity reading for that measuring point. For example, the outbound gravity value of measuring point P1 is 2834.5 mGal (relative to the instrument's internal reference), and the corresponding outbound measurement time is 8:15:30; the outbound gravity value of measuring point P25 is 2891.2 mGal, and the outbound measurement time is 14:30:45. Return measurements begin at 15:00 on the same day, and are performed in the reverse order of P25, P24...P1, A. The return gravity value of measuring point P25 is 2891.3 mGal, and the return measurement time is 15:00:30; the return gravity value of measuring point P1 is 2834.6 mGal, and the return measurement time is 21:15:20.

[0036] In the above embodiment, the instantaneous drift rate of the target is calculated using round-trip measurement data. Taking measuring point P1 as an example, the difference between the return gravity value of 2834.6 mGal and the outbound gravity value of 2834.5 mGal is 0.1 mGal, and the round-trip time interval is 12 hours, 59 minutes, and 50 seconds (46790 seconds). The initial drift rate of this point is calculated to be 0.01 mGal / hour. Considering the gravimeter's travel time between measuring points (average 15 minutes) and the actual observation period at each measuring point (actual observation period = arrival time at the measuring point + dwell time), a continuous drift rate time function is established through cubic spline interpolation. The drift rate function can be fitted using a second-order polynomial: v(t) = a0 + a1 × t + a2 × t², where t is the time interval (in hours) relative to the start of the measurement. The fitting coefficients are a0 = 0.01 mGal / h, a1 = 0.001 mGal / h², and a2 = -0.0002 mGal / h³. This drift rate function reflects the nonlinear variation of the gravimeter's drift rate over time. The initial drift rate is small, gradually increasing with the increase in operating time, which conforms to the typical drift characteristics of the CG-6 gravimeter. The target cumulative drift is calculated by definite integral of the drift rate function: D(t) = ∫0 tv(τ)dτ=a0×t+(a1 / 2)×t²+(a2 / 3)×t³. For example, the cumulative drift of measuring point P10 at the outbound measurement time (3.5 hours relative to the starting time) is 0.04mGal, and the cumulative drift at the return measurement time (10.2 hours relative to the starting time) is 0.13mGal.

[0037] In the above embodiment, the time weighting coefficient is determined based on the measurement time ratio. The total measurement time for the outbound journey is 6 hours, 30 minutes, and 45 seconds, and the total measurement time for the return journey is 6 hours, 14 minutes, and 50 seconds. For measurement point P10, the outbound measurement time ratio is 3.5 / 6.51 = 0.54, and the return measurement time ratio is (6.25 - 3.7) / 6.25 = 0.41. The time weighting coefficient uses linear interpolation: the outbound weight is 0.45, and the return weight is 0.55, ensuring reasonable fusion of the round-trip measurement data. The time drift correction includes two parts: the main drift correction directly uses the target cumulative drift; the closure error correction is allocated according to the time ratio using the round-trip gravity closure error (the round-trip difference of 0.1mGal at base point A). The outward time drift correction for measuring point P10 is 0.04 + 0.1 × 0.54 = 0.09 mGal, and the return time drift correction is 0.13 + 0.1 × 0.41 = 0.17 mGal. Finally, the outward corrected gravity value for measuring point P10 is 2856.4 - 0.09 = 2856.31 mGal, and the return corrected gravity value is 2856.5 - 0.17 = 2856.33 mGal. The target corrected gravity value is obtained by weighted averaging: 2856.31 × 0.45 + 2856.33 × 0.55 = 2856.32 mGal. This value represents the gravity difference between measuring point P10 and the secondary gravity benchmark A, eliminating the influence of instrument drift and achieving a measurement accuracy of ±0.08 mGal, meeting the accuracy requirements for gravity measurements in mineral exploration. In practical applications, the required measurement accuracy is adjusted according to different measurement task needs: high-precision mode (0.06-0.09 mGal) is suitable for fine geological structure detection and oil and gas exploration; standard-precision mode (0.09-0.15 mGal) is suitable for general mineral exploration and regional gravity surveys; and rapid measurement mode (above 0.15 mGal) is suitable for preliminary exploration and large-scale gravity surveys. This embodiment uses high-precision mode as an example, improving measurement accuracy by extending the observation time and increasing the number of readings.

[0038] Through the above steps, the outbound and return gravity measurements are performed in reverse order along the same measuring line, ensuring a symmetrical temporal distribution of the outbound and return measurement times at each measuring point. This symmetry provides a reliable time reference for subsequently determining the instantaneous drift rate of the target. A drift rate function is established based on the difference between the return and outbound gravity values ​​and the corresponding time relationships of the measurement times. This function accurately reflects the drift variation of the gravimeter throughout the measurement process. The cumulative drift is obtained through definite integral calculation starting from the measurement initiation time, enabling continuous cumulative calculation of the drift and avoiding discontinuities that may arise from discrete point correction. The time weighting coefficient, combined with the ratio of outbound and return measurement times, allows for a reasonable allocation of drift influence in the time drift correction. Finally, the corrected target gravity value is obtained through comprehensive correction of the outbound and return gravity values, achieving precise compensation for drift errors. This solves the technical problem of low correction accuracy for gravimeter drift measurement in related technologies, achieving a significant improvement in the correction accuracy of gravimeter drift measurement.

[0039] The entity performing the above steps may be a system with gravimeter drift correction capability, or a device with gravimeter drift correction capability, or a controller or processor in the device or system, or a standalone controller or processor, or other processing devices or processing units with similar processing functions, but is not limited to these.

[0040] In an optional embodiment, the target instantaneous drift rate of each measuring point is determined based on the return gravity value, the outbound gravity value, the outbound measurement time, and the return measurement time. Specifically, this includes: subtracting the return gravity value from the outbound gravity value to obtain the gravity difference at each measuring point; calculating a first difference between the return measurement time and the outbound measurement time to obtain the round-trip time interval at each measuring point; dividing the gravity difference by the round-trip time interval to obtain the average drift rate at each measuring point; obtaining the movement time of the gravimeter between every two adjacent measuring points; and obtaining the dwell observation time of the gravimeter at each measuring point; and determining the drift rate based on the movement time and dwell observation time. The actual observation period for each measuring point is determined; the average drift rate of each measuring point is correlated with the midpoint of the corresponding actual observation period to construct a discrete drift rate dataset; interpolation is performed on the discrete drift rate dataset to obtain a continuous drift rate time function; the function value of the midpoint of the actual observation period for each measuring point is extracted from the drift rate time function, and the function value is used as the initial instantaneous drift rate of each measuring point; historical drift characteristic data of the gravimeter is obtained, and drift correction coefficients are determined based on the historical drift characteristic data; the initial instantaneous drift rate is corrected using the drift correction coefficients to obtain the target instantaneous drift rate.

[0041] In the above embodiment, a CG-6 relative gravimeter was used to conduct detailed gravity measurements on an iron ore area. First, the gravity difference at each measuring point was calculated. Taking measuring point P5 as an example, the return gravity value was 2847.9 mGal, and the outward gravity value was 2847.8 mGal. Subtracting the two yielded a gravity difference of 0.1 mGal. This gravity difference reflects the total drift of the instrument during the round-trip measurement, and is a relative change value relative to the instrument's internal quartz spring system reference. The round-trip time interval was calculated using precise time recording. The outward measurement time for measuring point P5 was 9:42:15, and the return measurement time was 19:28:40, with a time interval of 9 hours, 46 minutes, and 25 seconds (35185 seconds). Dividing the gravity difference of 0.1 mGal by the time interval of 9.773 hours yielded an average drift rate of 0.0102 mGal / hour. This average drift rate represents the average drift characteristics of this measuring point during the round-trip measurement. The movement time of the gravimeter between adjacent measuring points is accurately obtained through GPS (Global Positioning System) trajectory recording. For example, the movement time from measuring point P4 to P5 is 12 minutes and 30 seconds, and the movement time from P5 to P6 is 13 minutes and 45 seconds. The dwell observation time is the actual measurement duration of the instrument at each measuring point. The CG-6 gravimeter dwells at each measuring point for 120 seconds for automatic observation. The actual observation period refers to the complete time period from arriving at the measuring point and setting up the instrument to completing the observation and leaving. The actual observation period for measuring point P5 is from 9:40:00 to 9:43:00, with the midpoint time being 9:41:30.

[0042] In the above embodiment, the discrete drift rate dataset is constructed by pairing the average drift rate of 25 measurement points with the midpoint of their corresponding actual observation periods. For example, the midpoint time of measurement point P1 at 8:16:00 corresponds to a drift rate of 0.0094 mGal / h, the midpoint time of measurement point P5 at 9:41:30 corresponds to 0.0102 mGal / h, and the midpoint time of measurement point P10 at 11:23:00 corresponds to 0.0111 mGal / h. These discrete data points reflect the trend of drift rate change over time. Cubic spline interpolation can be used to ensure the second-order continuity of the drift rate time function. The continuous function obtained after interpolation is: v(t) = 0.0086 + 0.0014 × t + 0.000065 × t² - 0.0000032 × t³, where t is the time relative to the measurement start time (in hours). The function value at the midpoint of the actual observation period for each measuring point was extracted from the drift rate time function. The function value for measuring point P5 at t=1.692 hours was 0.0100 mGal / h, which was taken as the initial instantaneous drift rate for that point. Historical drift characteristic data were obtained from the long-term monitoring records of the gravimeter. Historical data from the CG-6 gravimeter showed that the drift rate was relatively stable, approximately 0.008-0.011 mGal / h, during the first 4 hours of continuous operation. After 4-8 hours of operation, due to the temperature effect of the quartz spring and material fatigue, the drift rate increased to 0.011-0.016 mGal / h. After more than 8 hours, the drift rate may reach 0.016-0.022 mGal / h. Based on these historical characteristics, a drift correction coefficient was determined.

[0043] In the above embodiments, the drift correction coefficient was determined considering the instrument's cumulative operating time and environmental factors. In this measurement, the measurement time at measuring point P5 was 2.5 hours after the instrument was powered on. Based on the drift characteristics under the same operating time in the built-in database, the correction coefficient K was determined to be 1.02. This means that the initial instantaneous drift rate needs to be corrected upwards by 2%. Therefore, the target instantaneous drift rate of measuring point P5 is 0.0100 × 1.02 = 0.0102 mGal / h. Through the above refinement process, each measuring point obtained a target instantaneous drift rate that takes into account instrument characteristics and time effects. For example, the target instantaneous drift rate of the early measuring point P1 is 0.0096 mGal / h, the target instantaneous drift rate of the mid-term measuring point P12 is 0.0118 mGal / h, and the target instantaneous drift rate of the late measuring point P25 is 0.0142 mGal / h, showing an increasing trend consistent with the drift characteristics of the CG-6 gravimeter.

[0044] In an optional embodiment, determining the drift correction coefficient based on historical drift characteristic data specifically includes: extracting the drift rate variation curves of the gravimeter under different working durations from the historical drift characteristic data; determining the drift acceleration characteristic value of the gravimeter based on the drift rate variation curves; obtaining the cumulative working duration of the current measurement task and determining the drift rate value corresponding to the cumulative working duration from the drift rate variation curves; comparing the drift rate value with a preset standard drift rate to obtain a drift rate ratio; and determining the drift correction coefficient K based on the drift acceleration characteristic value and the drift rate ratio using the following formula:

[0045] K = 1 + (R - 1) × (1 + α × T / T0)

[0046] Where K is the drift correction coefficient, α is the drift acceleration characteristic value, R is the drift rate ratio, T is the cumulative working time, and T0 is the preset standard working time.

[0047] In the above embodiment, a CG-6 relative gravimeter was used for gravity measurement in the mining area. This gravimeter is equipped with a built-in database and can continuously record 180 days of historical measurement data. Historical drift characteristic data refers to the set of data recording the drift rate changes over time in past measurement tasks, including drift performance under different temperature conditions, different working durations, and different measurement environments. Drift rate variation curves of the gravimeter under different working durations were extracted from the built-in database. These curves show that: during the first 0-2 hours of operation, the drift rate remained at 0.008-0.010 mGal / h; during 2-4 hours, the drift rate increased to 0.010-0.012 mGal / h; during 4-6 hours, the drift rate reached 0.012-0.015 mGal / h; during 6-8 hours, the drift rate rose to 0.015-0.018 mGal / h; and after more than 8 hours, the drift rate reached 0.018-0.022 mGal / h. This increasing trend is mainly caused by the temperature response and material elasticity changes of the quartz spring system. The drift acceleration characteristic value α is obtained by taking the second derivative of the drift rate change curve. In the specific calculation, the drift rate values ​​at adjacent time points on the curve are selected, and the rate of change of the rate is calculated. For example, in the working time interval of 2-4 hours, the drift rate increases from 0.010 mGal / h to 0.012 mGal / h, with a rate of change of 0.001 mGal / h²; in the 4-6 hour interval, the rate of change is 0.0015 mGal / h²; and in the 6-8 hour interval, the rate of change is 0.0015 mGal / h². The average value of these rates of change is taken and normalized to obtain the drift acceleration characteristic value α = 0.00133.

[0048] In the above embodiment, the cumulative working time of the current measurement task refers to the total time from the gravimeter's power-on to the current measurement time at the measuring point. Taking measuring point P15 as an example, the outbound measurement time of this point is 12:45:30, which is 5.25 hours since the power-on at 7:30 AM. From the historical drift rate change curve, the drift rate value corresponding to 5.25 hours is determined to be 0.0135 mGal / h through linear interpolation. The preset standard drift rate is a reference value determined according to the technical specifications of the CG-6 gravimeter, set at 0.011 mGal / h, representing the typical drift rate of the instrument under standard working conditions (temperature 20℃, working time 4 hours). Comparing the actual drift rate value of 0.0135 mGal / h with the standard value of 0.011 mGal / h, the drift rate ratio R = 0.0135 / 0.011 = 1.2273 is obtained. The preset standard working time T0 is set to 4 hours (based on the standard working conditions of the gravimeter). The current cumulative working time T = 5.25 hours. Substituting into the drift correction coefficient formula K = 1 + (R-1) × (1 + α × T / T0) = 1 + (1.2273 - 1) × (1 + 0.00133 × 5.25 / 4) = 1.2304, this correction coefficient K = 1.2304 indicates that the initial instantaneous drift rate needs to be corrected upward by 23.04% to compensate for the drift acceleration effect after the instrument has been working for a long time. In practical applications, the drift correction coefficient for different measuring points will be dynamically adjusted according to the cumulative working time at the time of measurement. The correction coefficient K = 1.018 for early measuring points such as P1 (cumulative working time of 0.5 hours), K = 1.195 for mid-term measuring points such as P12 (cumulative working time of 3.5 hours), and K = 1.312 for late measuring points such as P25 (cumulative working time of 7 hours). Through this dynamic correction mechanism, each measuring point can obtain a drift correction coefficient that matches its measurement time, ensuring the time-varying adaptability of drift correction.

[0049] In an optional embodiment, a drift rate function is established using the time interval between the outward measurement time and the outward measurement start time as the independent variable and the target instantaneous drift rate as the dependent variable. Specifically, this includes: calculating a second difference between the outward measurement time and the outward measurement start time to obtain the outward measurement time interval for each measuring point; determining the instantaneous drift rate difference between every two adjacent measuring points and the outward measurement time interval difference between every two adjacent measuring points; determining the drift acceleration based on the instantaneous drift rate difference and the outward measurement time interval difference; establishing an initial drift rate function using a preset numerical fitting method with the outward measurement time interval as the independent variable and the target instantaneous drift rate as the dependent variable; acquiring the working state parameters of the gravimeter and determining the abnormal drift periods of the gravimeter based on the working state parameters; extracting the drift characteristics of the gravimeter from the abnormal drift periods and adjusting the weights of the abnormal function portion within the abnormal drift periods in the initial drift rate function based on the drift characteristics to obtain the drift rate function.

[0050] In the above embodiment, a CG-6 relative gravimeter was used to perform gravity measurements in the copper mine area. First, the outward measurement time interval for each measuring point was calculated. Taking measuring point P8 as an example, the outward measurement time was 10:45:20, and the outward measurement start time was 8:00:00. The difference between the two was 2 hours, 45 minutes, and 20 seconds, or 2.756 hours. This time interval reflects the length of time elapsed from the start of the measurement to the outward measurement at that measuring point. The difference in instantaneous drift rate was determined by comparing the target instantaneous drift rates of adjacent measuring points. For example, the target instantaneous drift rate of measuring point P7 was 0.0104 mGal / h, and that of measuring point P8 was 0.0106 mGal / h, with a difference of 0.0002 mGal / h. Correspondingly, the outward measurement time interval for P7 was 2.423 hours, and for P8 it was 2.756 hours, with a time interval difference of 0.333 hours. The drift acceleration is calculated by dividing the difference in instantaneous drift rates by the difference in time intervals: 0.0002 / 0.333 = 0.0006 mGal / h². The preset numerical fitting method can be least squares polynomial fitting. Using the outward measurement time interval of the 25 measuring points as the independent variable t and the target instantaneous drift rate as the dependent variable v, a cubic polynomial function of the initial drift rate is obtained: v(t) = 0.0085 + 0.0013t + 0.000084t² - 0.0000049t³. This function describes the nonlinear change in drift rate over time, where the constant term 0.0085 mGal / h represents the initial drift rate, the coefficient of the first term 0.0013 reflects the linear growth trend, and the coefficients of the second and third terms describe the nonlinear characteristics.

[0051] In the above embodiments, the operating status parameters include key indicators such as the instrument's internal temperature, tilt compensation status, and battery voltage. The internal temperature sensor of the CG-6 gravimeter detected that during the 4.5 to 5.2 hour measurement period, due to rapid changes in ambient temperature (from 18°C ​​to 26°C), the instrument's internal temperature control system entered an adjustment state, with temperature fluctuations reaching ±0.3°C. The tilt compensation system recorded that at the 3.8 hour, due to ground vibration, the instantaneous change in the X-axis tilt angle exceeded 15 arcseconds. Battery voltage monitoring showed a stable range of 12.4-12.6V throughout the measurement, without any abnormalities. Abnormal drift periods were determined based on the abnormal values ​​of the operating status parameters. The temperature adjustment period (4.5-5.2 hours) was marked as the first abnormal period, during which the drift rate exhibited significant fluctuations, deviating from the normal trend by 0.0022 mGal / h. The tilt disturbance moment (around 3.8 hours ±0.1 hours) was marked as the second abnormal period, during which the drift rate experienced brief jumps. The measurement points corresponding to these abnormal periods include P11, P12, P13 (temperature anomalies) and P9, P10 (tilt disturbances). Drift characteristics were extracted from the data analysis of these abnormal periods. The drift characteristics during the temperature anomaly periods exhibited periodic fluctuations, with the fluctuation amplitude directly proportional to the rate of temperature change, and a correlation coefficient of 0.82. The drift characteristics during the tilt disturbance periods exhibited pulse-like jumps, with the jump amplitude correlated with the change in tilt angle. Based on these characteristics, the initial drift rate function was corrected.

[0052] In the above embodiments, the weight adjustment can employ a segmented weighting method. For periods of abnormal temperature (4.5-5.2 hours), the weight coefficient of the function for this period is set to 0.6, meaning that the reliability of the data during this period is reduced by 40%. For periods of tilted disturbance (3.7-3.9 hours), the weight coefficient is set to 0.7. The weight coefficient for normal periods remains at 1.0. The adjusted drift rate function is expressed as: v(t) = w(t) × [0.0085 + 0.0013t + 0.000084t² - 0.0000049t³], where the weight function w(t) is defined as: w(t) = 1.0 when t ∈ [0, 3.7) ∪ (3.9, 4.5) ∪ (5.2, 8]; w(t) = 0.7 when t ∈ [3.7, 3.9]; w(t) = 0.6 when t ∈ [4.5, 5.2]. Through this weight adjustment mechanism, the influence of abnormal periods on the overall drift rate function is reasonably suppressed, while preserving its trend information. The adjusted drift rate function more realistically reflects the drift characteristics of the gravimeter under normal operating conditions, achieving accurate modeling of drift characteristics under complex measurement environments.

[0053] In an optional embodiment, starting from the measurement start time, a first definite integral operation is performed on the drift rate function along the time progression direction to obtain the target cumulative drift amount at each measurement point. Specifically, this includes: determining the first derivative of the drift rate function and comparing the absolute value of the first derivative with a preset gradient threshold to obtain a gradient comparison result; dividing the drift rate function into a stable period and a rapidly changing period based on the gradient comparison result; performing a first numerical integration calculation on the stable function portion of the drift rate function located within the stable period using a first integration step size to obtain the stable period integral value; performing a second numerical integration calculation on the rapidly changing function portion of the drift rate function located within the rapidly changing period using a second integration step size to obtain the rapidly changing period integral value, wherein the second integration step size is smaller than the first integration step size; and accumulating the stable period integral value and the rapidly changing period integral value according to the time sequence of each measurement point on the time axis to obtain the target cumulative drift amount.

[0054] In the above embodiment, a CG-6 relative gravimeter was used to measure gravity in the mining area, and the established drift rate function v(t) = 0.0124 + 0.0019t + 0.000124t² - 0.0000072t³ was integrated. First, the first derivative of this function was calculated, yielding v'(t) = 0.0019 + 0.000248t - 0.0000216t², which characterizes the rate of change of the drift rate. The preset gradient threshold was set to 0.00032 mGal / h², which was determined based on the typical range of drift rate changes of the CG-6 gravimeter under normal operating conditions. The gradient comparison results were obtained by calculating the absolute value of the derivative |v'(t)| point by point and comparing it with the threshold. For example, at t=1.5 hours, |v'(1.5)|=|0.0019+0.000372-0.0000486|≈0.002223, which is greater than the threshold; at t=4.8 hours, |v'(4.8)|=|0.0019+0.001190-0.000498|≈0.002592, which is greater than the threshold; at t=6.2 hours, |v'(6.2)|=|0.0019+0.001538-0.000832|≈0.002606, which is greater than the threshold. Based on the gradient comparison results, the 8-hour measurement period was divided into two phases: a stable phase, including hours [0, 1.2] and [6.8, 8], during which the absolute values ​​of the derivatives were relatively small; and a rapidly changing phase, including hours [1.2, 6.8], which corresponds to the instrument's main operating period and shows significant drift rate changes. The stable phase accounted for 30% of the total duration, while the rapidly changing phase accounted for 70%.

[0055] In the above embodiment, the first integration step size is set to 0.05 hours (3 minutes), suitable for integration calculations during stationary periods. Using the trapezoidal rule, the stationary function portion for the [0, 1.2] hour period is calculated: D1 = ∫0¹·²v(t)dt ≈ Σᵢ₌0²³[v(tᵢ) + v(tᵢ₊1)] / 2 × 0.05, yielding the first stationary period integration value D1 = 0.2 mGal. Similarly, for the [6.8, 8] hour period, the second stationary period integration value D2 = 0.3 mGal is obtained. The second integration step size is set to 0.01 hours (36 seconds), used for fine integration during rapidly varying periods. For the rapidly varying function portion for the [1.2, 6.8] hour period, Simpson's integral method is used: D3 = ∫1.2 6 · 8 v(t)dt≈Σ j ₌0 559 [v(t 2j )+4v(t 2j+1 )+v(t 2j+2 )] / 6×0.01, the integral value D3=1.1mGal is calculated for the rapidly changing period. The use of Simpson's integral method improves the integration accuracy in the rapidly changing drift rate range. The cumulative drift of the target is obtained by accumulating the values ​​sequentially over time. For measuring point P10 (measurement time is 2.8 hours), its cumulative drift is the integral value before that time: C(2.8)=D1+∫1.2²· 8 v(t)dt=0.2+0.3=0.5mGal. For measuring point P15 (measurement time is 4.2 hours), the integral value of the stationary period [0,1.2] and the partial integral value of the rapidly changing period [1.2,4.2] need to be accumulated: C(4.2)=D1+∫1.2 4 ·²v(t)dt=0.2+0.6=0.8mGal.

[0056] In the above embodiment, the calculation of cumulative drift needs to be segmented for measurement points spanning different time periods. The cumulative drift of measurement point P20 (measured at hour 6.5) is: C(6.5) = D1 + ∫1.2 6 · 5 v(t)dt=0.2+1.0=1.2mGal. The cumulative drift at measuring point P25 (measurement time is 7.8 hours) is: C(7.8)=D1+D3+∫6.8 7 · 8v(t)dt = 0.2 + 1.1 + 0.3 = 1.6 mGal. During the integration calculation, the data acquisition system of the CG-6 gravimeter recorded the raw data at a frequency of 1 Hz (Hertz), providing ample data support for numerical integration. An adaptive step-size integration strategy optimized computational efficiency while maintaining accuracy. The final cumulative drift values ​​for the 25 measurement points ranged from 0.1 to 1.6 mGal, with an average of approximately 0.9 mGal. These values ​​accurately reflect the cumulative drift characteristics of the gravimeter throughout the measurement process, providing a reliable quantitative basis for subsequent drift correction.

[0057] In an optional embodiment, the time weighting coefficient for each measurement point is determined based on the proportion of outbound measurement time in the total outbound measurement duration and the proportion of return measurement time in the total return measurement duration. The time drift correction amount for each measurement point is then determined based on the target cumulative drift and the time weighting coefficient. Specifically, this includes: obtaining the outbound measurement end time and performing a third difference calculation between the outbound measurement end time and the outbound measurement start time to obtain the total outbound measurement duration; obtaining the return measurement start time and the return measurement end time and performing a fourth difference calculation between the return measurement end time and the return measurement start time to obtain the total return measurement duration; and dividing the outbound measurement time interval by... The total outbound measurement time is used to obtain the outbound measurement time ratio; the fifth difference between the end time and the return measurement time is calculated to obtain the return measurement time interval; the return measurement time interval is divided by the total return measurement time to obtain the return measurement time ratio; the cumulative target drift is used as the theoretical drift correction; the time weighting coefficient is determined based on the outbound and return measurement time ratios; the outbound and return theoretical corrections are determined based on the theoretical drift correction and the time weighting coefficient; the round-trip gravity values ​​are compared at multiple measurement points to obtain the round-trip gravity closure difference; the outbound and return theoretical corrections are corrected based on the round-trip gravity closure difference to obtain the time drift correction.

[0058] In the above embodiment, a CG-6 relative gravimeter was used to conduct gravity measurements in an iron ore area. The outbound measurement started at 7:30:00 and ended at 11:45:30, with the difference resulting in a total outbound measurement time of 4.258 hours. The return measurement started at 12:15:00 (including a 30-minute rest period) and ended at 16:28:45, with the difference resulting in a total return measurement time of 4.229 hours. Taking measuring point P12 as an example, its outbound measurement time was 9:42:18, and the time interval from the outbound start time was 2.206 hours. Dividing this by the total outbound measurement time of 4.258 hours yields an outbound measurement time ratio of 0.518. The return measurement time for this measuring point is 14:16:27. The difference between the return measurement end time (16:28:45) and the return measurement time is 2.204 hours. Dividing this return measurement time interval by the total return measurement duration (4.229 hours) yields a return measurement time ratio of 0.521. The calculation of the time ratio reflects the relative position of the measuring point throughout the entire measurement process. The theoretical drift correction amount is directly adopted from the previously calculated target cumulative drift amount. The time interval between the outward measurement time and the outward start time for measuring point P12 is 2.206 hours. Integrating the drift rate function v(t) = 0.0085 + 0.0013t + 0.000084t² - 0.0000049t³ over the interval [0, 2.206] yields an outward target cumulative drift amount of 0.62mGal. The time interval between the return measurement time and the return start time at this measuring point is 2.204 hours. Integrating the drift rate function over the interval [0, 2.204] yields a cumulative drift of 0.61 mGal. These two cumulative drift values ​​reflect the cumulative drift effect of the gravimeter's quartz spring system under different operating durations.

[0059] In the above embodiments, the time weighting coefficient is determined based on the similarity of the outbound and return measurement time ratios. For measurement points with a time ratio difference of less than 0.05, the average weighting method is used: outbound weight w f =0.5, return trip weight w b =0.5. For measuring points with differences between 0.05 and 0.15, the linear adjustment method is used: w f =0.5-0.2×(|P f -P b | / 0.15), w b =1-w f , where P f P represents the proportion of the outbound measurement time. b This represents the proportion of the return measurement time. The proportion difference at measuring point P12 is |0.518-0.521|=0.003, which is less than 0.05, therefore w f =w b=0.5. The theoretical correction amounts for the outward and return journeys are determined based on the target cumulative drift. In this embodiment, the theoretical correction amount is directly equal to the target cumulative drift amount, without involving adjustments to the time weighting coefficient. The theoretical correction amount for the outward journey at measuring point P12 is Cf,12 = 0.62 mGal, and the theoretical correction amount for the return journey is Cb,12 = 0.61 mGal. These two theoretical correction amounts reflect the theoretical prediction of the drift influence at the outward and return measurement times of this measuring point based on the drift rate function model. The round-trip gravity closure error is obtained by comparing the relative gravity values ​​for the outward and return journeys at the same measuring point. The relative gravity value for the outward journey at measuring point P12 is 1842.4 mGal, and the relative gravity value for the return journey is 1842.9 mGal, with a difference of 0.5 mGal. After applying the theoretical corrections, the temporary correction gravity value for the outbound journey is 1842.4 - 0.62 = 1841.78 mGal, and the temporary correction gravity value for the return journey is 1842.9 - 0.61 = 1842.29 mGal. The round-trip consistency deviation is |1841.78 - 1842.29| = 0.51 mGal. The round-trip consistency deviation ranges from 0.28 to 0.85 mGal for the 25 measuring points along the entire survey line, with a root mean square error (RMSE) of 0.62 mGal. The round-trip consistency deviation is mainly caused by three factors: the nonlinear characteristics of instrument drift, the influence of ambient temperature changes, and the interference of ground micro-vibrations. The round-trip consistency deviation reflects the discrepancy between the theoretical drift model and the actual measurement, requiring correction of the theoretical correction values ​​for both the outbound and return journeys.

[0060] In the above embodiments, the specific process of correcting the theoretical correction amounts for the outbound and return journeys based on the round-trip consistency deviation is based on the least squares principle. The goal of the correction is to make the corrected outbound and return journey correction gravity values ​​as close as possible to the true gravity values, while minimizing the deviation of the correction amounts from the theoretical model. The correction model is established as follows: Let the outbound correction amount for measuring point i be ΔCf,i, and the return correction amount be ΔCb,i. Then the corrected outbound correction gravity value is Gf,i-Cf,i-ΔCf,i, and the return correction gravity value is Gb,i-Cb,i-ΔCb,i, where Gf,i and Gb,i are the relative gravity values ​​for the outbound and return journeys, respectively, and Cf,i and Cb,i are the theoretical correction amounts for the outbound and return journeys, respectively. The modified constraints include: 1) Round trip consistency constraint: The corrected gravity values ​​for the outbound and return trips should be equal, i.e. (Gf,i-Cf,i-ΔCf,i)=(Gb,i-Cb,i-ΔCb,i), which can be transformed to obtain ΔCf,i-ΔCb,i=(Gf,i-Cf,i)-(Gb,i-Cb,i). Define the round-trip consistency deviation as Δi=(Gf,i-Cf,i)-(Gb,i-Cb,i), then the constraint condition simplifies to ΔCf,i-ΔCb,i=Δi; 2) Minimize correction amount constraint: The correction amount should be as small as possible to maintain the effectiveness of the theoretical model. The objective function is minΣ(ΔCf,i²+ΔCb,i²); 3) Time weight constraint: The allocation of correction amount should consider the measurement time sequence and measurement point location. The outbound measurement time is earlier, the theoretical model prediction accuracy is higher, and the correction amount should be smaller; the return measurement time is later, the uncertainty of cumulative drift increases, and the correction amount should be larger.

[0061] In the above embodiments, based on the aforementioned constraints, the weighted least squares method is used to solve for the optimal correction amount. A time weighting coefficient w is introduced. f and w b As a weighting factor for the adjustment allocation, a weighted Lagrangian function is constructed: L=Σ(w f ·ΔCf,i²+w b Taking the partial derivatives of ΔCf,i and ΔCb,i with respect to each and setting them to zero, we get: ∂L / ∂ΔCf,i = 2w f ·ΔCf,i+λ=0, ∂L / ∂ΔCb,i=2w b ·ΔCb,i-λ=0. Solving the system of equations simultaneously, we get: ΔCf,i=-λ / (2w f ), ΔCb,i=λ / (2w b Substituting this into the round-trip consistency constraint ΔCf,i-ΔCb,i=Δi, we get: -λ / (2w f )-λ / (2w b )=Δi, that is, λ=-Δi / [1 / (2w f)+1 / (2w b )]=-Δi·w f ·w b / (w f +w b Therefore, the optimal correction is: ΔCf,i = Δi·w b / (w f +w b ), ΔCb,i=-Δi·w f / (w f +w b This result indicates that, in the weighted least squares sense, round-trip consistency bias should be allocated inversely proportional to the time weighting coefficient. When w f =w b When = 0.5, the formula simplifies to: ΔCf,i = Δi / 2, ΔCb,i = -Δi / 2, i.e., equal distribution. When w f >w b When the weight of the outbound journey is larger, the correction amount of the outbound journey is smaller (ΔCf,i<Δi / 2), while the correction amount of the return journey is larger (|ΔCb,i|>Δi / 2). This is consistent with the physical intuition that "the outbound journey measurement is more reliable and the correction amount should be smaller".

[0062] In the above embodiment, for measuring point P12, the round-trip consistency deviation Δ12 = (1842.4 - 0.62) - (1842.9 - 0.61) = 1841.78 - 1842.29 = -0.51mGal (a negative value indicates that the outbound correction gravity value is less than the return correction gravity value). Time weighting coefficient w f =w b=0.5. According to the weighted least squares correction model: the outbound correction ΔCf,12 = -0.51 × 0.5 / (0.5 + 0.5) = -0.51 / 2 = -0.255mGal; the return correction ΔCb,12 = -(-0.51) × 0.5 / (0.5 + 0.5) = 0.51 / 2 = 0.255mGal. The time drift correction is obtained by correcting the theoretical correction for round-trip consistency deviation. For outbound measurement: outbound time drift correction = outbound theoretical correction + outbound correction, for measuring point P12 outbound: 0.62 + (-0.255) = 0.365mGal. For return measurement: return time drift correction = return theoretical correction + return correction, for measuring point P12 return: 0.61 + 0.255 = 0.865mGal. The corrected outbound gravity value is 1842.4 - 0.365 = 1842.035 mGal, and the corrected return gravity value is 1842.9 - 0.865 = 1842.035 mGal. The two values ​​are completely consistent, reducing the round-trip consistency deviation to 0 mGal. This verifies the effectiveness of the weighted least squares correction method. The corrected outbound and return gravity values ​​are weighted and averaged according to the time weighting coefficient to obtain the final corrected gravity value for this measuring point: Target corrected gravity value = w f ×1842.035+w b ×1842.035=0.5×1842.035+0.5×1842.035=1842.035mGal. Since the round-trip consistency deviation is zero after correction, the weighted average result is the same as the single value.

[0063] In the above embodiment, the same correction method was applied to all 25 measuring points along the entire measurement line. Before correction, the root mean square error (RMSE0) of the round-trip consistency deviation was 0.62 mGal; after correction, the RMSE1 of the round-trip consistency deviation was 0.08 mGal, a reduction of 87%. The corrected round-trip consistency deviation mainly originated from random errors and environmental disturbances during the measurement process, and is now close to the reading resolution of the CG-6 gravimeter (0.1 mGal). Through this correction method based on the weighted least squares principle, each measuring point obtained a drift correction value that matched its measurement sequence and location. The correction process has a clear mathematical model and physical meaning: the weighted least squares principle ensures the optimality of the correction amount, the time weight coefficient reflects the characteristics of drift accumulation over time and the differences in measurement reliability, and the combination of the two ensures both the theoretical continuity of drift correction and eliminates the systematic bias of the theoretical model through feedback correction of the measured round-trip consistency deviation, ensuring the accuracy of relative gravity measurement conversion.

[0064] In an optional embodiment, the outbound and return gravity values ​​are corrected using a time drift correction to obtain a target corrected gravity value. The specific method further includes: subtracting the outbound time drift correction from the outbound gravity value to obtain a temporary outbound corrected gravity value; and subtracting the return time drift correction from the return gravity value to obtain a temporary return gravity value. Using the outbound and return temporary corrected gravity values ​​as iterative variables, the following iterative optimization operation is performed until the convergence characteristic value of the round-trip consistency deviation meets a preset convergence condition: calculating the sixth difference between the outbound and return temporary corrected gravity values ​​at the same measurement point to obtain the round-trip consistency deviation; and determining the convergence characteristic of the round-trip consistency deviation according to a preset convergence index. When the convergence characteristic value does not meet the preset convergence condition, the fitting parameters of the drift rate function are adjusted according to the round-trip consistency deviation to obtain the adjusted drift rate function; the adjusted drift rate function is subjected to a second definite integral operation to obtain a new target cumulative drift amount; a new outbound time drift correction amount and a new return time drift correction amount are determined according to the new target cumulative drift amount and the time weight coefficient; the outbound gravity value is drift corrected using the new outbound time drift correction amount to obtain a new outbound temporary correction gravity value, and the return gravity value is drift corrected using the new return time drift correction amount to obtain a new return temporary correction gravity value; the outbound temporary correction gravity value and the return temporary correction gravity value that meet the preset convergence condition are weighted and averaged according to the time weight coefficient to obtain the target correction gravity value.

[0065] In the above embodiment, a CG-6 relative gravimeter was used to conduct gravity measurements in a coal mining area. The measurement line consisted of 25 measuring points, using a round-trip observation method. The relative gravity value at measuring point P18 was 2156.8 mGal on the outward journey and 2157.7 mGal on the return journey. First, the initial time drift correction was calculated. The time interval between the outward measurement time and the outward start time at measuring point P18 was 2.8 hours. A definite integral was performed on the initial drift rate function v(t) = 0.0085 + 0.0013t + 0.000084t² - 0.0000049t³ in the interval [0, 2.8], yielding a cumulative outward drift of 0.92 mGal. The time interval between the return measurement time and the return start time at the same measuring point was 2.6 hours. A definite integral was performed on the initial drift rate function in the interval [0, 2.6], yielding a cumulative return drift of 0.88 mGal. The time weighting coefficient is determined based on the ratio of outbound and return measurement time. For measuring point P18, the outbound measurement time ratio is 0.52, and the return measurement time ratio is 0.54, with a ratio difference of 0.02, which is less than 0.05. Therefore, the average weighting method is used: w f =w b=0.5. According to the technical solution of claim 6, the theoretical correction amount for the outbound journey is equal to the cumulative drift amount of the outbound target (0.92 mGal), and the theoretical correction amount for the return journey is equal to the cumulative drift amount of the return target (0.88 mGal). After applying the round-trip consistency deviation correction, the outbound time drift correction amount for measuring point P18 is 0.90 mGal, and the return time drift correction amount is 0.90 mGal. Subtracting the outbound time drift correction amount from the outbound relative gravity value yields a temporary corrected gravity value of 2156.8 - 0.90 = 2155.90 mGal. Subtracting the return time drift correction amount from the return relative gravity value yields a temporary corrected gravity value of 2157.7 - 0.90 = 2156.80 mGal. The difference between these two temporary correction values ​​reflects the initial deviation of the drift model.

[0066] In the above embodiment, the round-trip consistency deviation is calculated using the sixth difference. The round-trip consistency deviation of measuring point P18 is |2155.90-2156.80|=0.90mGal. The round-trip consistency deviation of the 25 measuring points along the entire measuring line ranges from 0.2 to 1.2mGal, reflecting the imperfections of the initial drift correction. The preset convergence index uses the root mean square error, calculated as: RMSE=√(Σ(Δᵢ²) / n), where Δᵢ is the round-trip consistency deviation of the i-th measuring point, and n is the total number of measuring points. The convergence characteristic value before iteration is RMSE0=0.78mGal. The preset convergence condition is set to the root mean square error of the round-trip consistency deviation <0.2mGal. This 0.2mGal threshold can be determined based on the 0.1mGal reading resolution of the CG-6 gravimeter and the actual measurement accuracy requirements. Since RMSE0=0.78mGal is greater than the convergence condition of 0.2mGal, parameter adjustment is required. The original form of the drift rate function is v(t) = 0.0085 + 0.0013t + 0.000084t² - 0.0000049t³. Based on the distribution characteristics of the round-trip consistency deviation, the gradient descent method is used to adjust the fitting parameters. The distribution characteristics of the round-trip consistency deviation are analyzed as follows: 1) Sign distribution of the deviation: Among the 25 measuring points, the return correction gravity value of 18 measuring points is greater than the outgoing correction gravity value (positive deviation), and the return correction gravity value of 7 measuring points is less than the outgoing correction gravity value (negative deviation). The positive deviation accounts for 72%, indicating that the drift rate function underestimates the actual drift as a whole; 2) Spatial distribution of the deviation: The average deviation of the first section of the measuring line (measuring points 1-8) is 0.45mGal, the average deviation of the middle section of the measuring line (measuring points 9-17) is 0.92mGal, and the average deviation of the last section of the measuring line (measuring points 18-25) is 0.68mGal. The middle section has the largest deviation, indicating that the quadratic term coefficient needs to be increased; 3) Temporal trend of the deviation: The deviation shows a trend of first increasing and then decreasing with the increase of measurement time, indicating that the cubic term coefficient needs to be adjusted.

[0067] In the above embodiments, the parameter adjustment strategy is based on the sensitivity analysis of each parameter to the round-trip consistency deviation. The sensitivity analysis is performed by calculating the partial derivatives of the deviation with respect to each parameter: ∂Δᵢ / ∂a0≈1, ∂Δᵢ / ∂a1≈tᵢ, ∂Δᵢ / ∂a2≈tᵢ², ∂Δᵢ / ∂a3≈tᵢ³. Based on the gradient descent method, the parameter adjustment amount is: Δa j =-α j ×Σ(Δᵢ×∂Δᵢ / ∂a j ) / Σ(∂Δᵢ / ∂a j )², where α j Here is the learning rate parameter. The specific calculations are as follows: Constant term adjustment: Δa0 = -α0 × Σ(Δᵢ × 1) / Σ(1²) = -0.8 × 18.5 / 25 = -0.592 / 25 = -0.00024. Where Σ(Δᵢ × 1) = 18.5mGal is the algebraic sum of the round-trip consistency deviations of all measurement points (considering the sign), and the learning rate α0 = 0.8. Linear term coefficient adjustment: Δa1 = -α1 × Σ(Δᵢ × tᵢ) / Σ(tᵢ²) = -0.6 × 52.3 / 1458 = -0.0000215. Where Σ(Δᵢ × tᵢ) = 52.3mGal·h is the weighted sum of deviation and time, Σ(tᵢ²) = 1458h² is the sum of squares of time, and the learning rate α1 = 0.6. Adjustment of quadratic term coefficient: Δa² = -α² × Σ(Δᵢ × tᵢ²) / Σ(tᵢ) 4 = -0.4 × (-128.6) / 68420 = 0.00000075. Where Σ(Δᵢ×tᵢ²) = -128.6mGal·h², Σ(tᵢ 4 )=68420h 4 The learning rate α² = 0.4. A negative value indicates a large mid-range bias, requiring an increase in the quadratic term coefficient. Cubic term coefficient adjustment: Δa³ = -α³ × Σ(Δᵢ × tᵢ³) / Σ(tᵢ) 6 )=-0.2×(-485.2) / 3125000=0.000000031. Among them, Σ(Δᵢ×tᵢ³)=-485.2mGal·h³, Σ(tᵢ 6 )=3125000h 6 The learning rate α3 = 0.2. The learning rate parameters α0, α1, α2, and α3 are set to 0.8, 0.6, 0.4, and 0.2, respectively, following a decreasing principle to ensure more careful adjustment of higher-order terms and avoid overfitting. The adjusted drift rate function is: v⁽¹⁾(t) = (0.0085 - 0.00024) + (0.0013 - 0.0000215)t + (0.000084 + 0.00000075)t² + (-0.0000049 + 0.000000031)t³ = 0.00826 + 0.001278t + 0.000085t² - 0.0000049t³.

[0068] In the above embodiment, a second definite integral operation is performed on the adjusted drift rate function using Simpson's integral method to obtain a new target cumulative drift amount. The time interval between the outbound measurement time and the start time of measurement point P18 is 2.8 hours. A definite integral is performed on the adjusted drift rate function in the interval [0, 2.8], yielding a new outbound cumulative drift amount of 0.87 mGal, a decrease of 0.05 mGal compared to the original value of 0.92 mGal. The time interval between the return measurement time and the return start time of the same measurement point is 2.6 hours. A definite integral is performed on the adjusted drift rate function in the interval [0, 2.6], yielding a new return cumulative drift amount of 0.83 mGal, a decrease of 0.05 mGal compared to the original value of 0.88 mGal. A new time drift correction amount is determined based on the new cumulative drift amount. During the iterative optimization process, the time drift correction amount is directly equal to the cumulative drift amount, without involving adjustments to the time weight coefficients or corrections for round-trip consistency deviations, thus maintaining the purity of the model optimization. The new outbound time drift correction is equal to the new outbound cumulative drift of 0.87 mGal; the new return time drift correction is equal to the new return cumulative drift of 0.83 mGal. This direct use of the cumulative drift as the correction ensures that the iterative process is entirely based on the optimization of the drift rate function model, avoiding interference from human factors in adjusting model parameters. After applying the new correction, the new outbound temporary correction gravity value for measuring point P18 is 2156.8 - 0.87 = 2155.93 mGal, and the new return temporary correction gravity value is 2157.7 - 0.83 = 2156.87 mGal. The new round-trip consistency deviation is reduced to |2155.93 - 2156.87| = 0.94 mGal. Although the deviation of a single measuring point increases slightly (from 0.90 to 0.94), the root mean square error of the entire measuring line decreases significantly. After the first iteration, RMSE1 = 0.51mGal, a 35% decrease compared to the initial value of 0.78mGal. This is because the parameter adjustment mainly optimized the middle section measurement points with larger deviations. The deviation of the middle section measurement points decreased from an average of 0.92mGal to 0.58mGal, a reduction of 37%. The deviations of the front and rear section measurement points increased slightly, but the overall RMSE decreased significantly.

[0069] In the above embodiment, the iterative process continues, and each iteration adjusts the drift rate function parameters based on the current round-trip consistency deviation distribution, recalculating the cumulative drift and correction. The second iteration uses the same gradient descent strategy, with parameter adjustments of: Δa0 = -0.00018, Δa1 = -0.0000162, Δa2 = 0.00000058, Δa3 = 0.000000024. The adjusted drift rate function is v⁽²⁾(t) = 0.00808 + 0.001262t + 0.000086t² - 0.0000048t³. After applying the correction from the second iteration, the temporary correction gravity value for the outbound journey at measuring point P18 is 2156.05 mGal, the temporary correction gravity value for the return journey is 2156.92 mGal, and the round-trip consistency deviation is 0.87 mGal. After the second iteration, RMSE2 = 0.32 mGal, a decrease of 37% compared to the first iteration. The parameter adjustments for the third iteration are: Δa0 = -0.00012, Δa1 = -0.0000108, Δa2 = 0.00000042, Δa3 = 0.000000018. The adjusted drift rate function is v⁽³⁾(t) = 0.00796 + 0.001251t + 0.000086t² - 0.0000048t³. After applying the corrections from the third iteration, the temporary corrected gravity value for the outbound journey at measuring point P18 is 2156.12 mGal, and the temporary corrected gravity value for the return journey is 2156.95 mGal, with a round-trip consistency deviation of 0.83 mGal. After the third iteration, RMSE3 = 0.19 mGal, a decrease of 41% compared to the second iteration. The parameter adjustments for the fourth iteration are: Δa0 = -0.00008, Δa1 = -0.0000072, Δa2 = 0.00000028, Δa3 = 0.000000012. The adjusted drift rate function is v⁽ 4⁾(t) = 0.00788 + 0.001244t + 0.000087t² - 0.0000047t³. After applying the correction value from the fourth iteration, the temporary correction gravity value for the outbound journey at measuring point P18 is 2156.16 mGal, and the temporary correction gravity value for the return journey is 2156.97 mGal, with a round-trip consistency deviation of 0.81 mGal. The RMSE4 after the fourth iteration is 0.15 mGal, a decrease of 21% compared to the third iteration, satisfying the preset convergence condition RMSE < 0.2 mGal. Throughout the iteration process, the parameters of the drift rate function are gradually optimized, and the final form is: v_final(t) = 0.00788 + 0.001244t + 0.000087t² - 0.0000047t³. The convergence of the iterative process is verified by the monotonically decreasing RMSE: RMSE0=0.78mGal→RMSE1=0.51mGal→RMSE2=0.32mGal→RMSE3=0.19mGal→RMSE4=0.15mGal, with a total decrease of 81%.

[0070] In the above embodiment, after the convergence condition is met, the final cumulative drift amount of each measuring point is calculated based on the final optimized drift rate function. The time interval between the outward measurement time and the start time of measuring point P18 is 2.8 hours. A definite integral is performed on the final drift rate function in the interval [0, 2.8], yielding a final outward cumulative drift amount of 0.82 mGal. The time interval between the return measurement time and the return start time of this measuring point is 2.6 hours. A definite integral is performed on the final drift rate function in the interval [0, 2.6], yielding a final return cumulative drift amount of 0.78 mGal. The original gravity value is corrected using the final cumulative drift amount: the converged outward temporary corrected gravity value = 2156.8 - 0.82 = 2155.98 mGal; the converged return temporary corrected gravity value = 2157.7 - 0.78 = 2156.92 mGal. At this point, the round-trip consistency deviation is |2155.98-2156.92|=0.94mGal. Although the round-trip consistency deviation of a single measuring point P18 slightly increased from the initial value of 0.90mGal to 0.94mGal, the root mean square error of the entire measuring line significantly decreased from 0.78mGal to 0.15mGal, a reduction of 81%, indicating that the iterative optimization method has a significant effect on the drift correction of the overall measuring line. The target corrected gravity value is obtained by weighting the outbound and return temporary correction gravity values ​​that meet the preset convergence conditions according to the time weight coefficient. The time weight coefficient for measuring point P18 is wf=wb=0.5, therefore the target corrected gravity value = 0.5×2155.98+0.5×2156.92=2156.45mGal. The role of the time weight coefficient is to combine the outbound and return measurement results after iterative convergence to obtain the final corrected gravity value. For measurement points with significant time ratio differences, the time weighting coefficient will tilt towards the side with more reliable measurement timing, thereby improving the accuracy of the final result. This target corrected gravity value is a relative gravity value compared to the starting baseline of the measurement line (i.e., the secondary gravity baseline). To obtain the absolute gravity value, the absolute gravity value of the starting baseline of the measurement line needs to be added. Assuming the absolute gravity value of the starting baseline of the measurement line, measured by an FG5 absolute gravimeter (accuracy better than 2 μGal), is 979254.3 mGal, then the absolute gravity value of measurement point P18 is 979254.3 + 2156.45 = 981410.75 mGal. In this embodiment, using the standard accuracy mode as an example, high consistency in round-trip measurements is achieved through iterative optimization of the drift rate function parameters.During the iteration process, the parameters of the drift rate function were gradually adjusted from the initial value v(t) = 0.0085 + 0.0013t + 0.000084t² - 0.0000049t³ to the optimal value v_final(t) = 0.00788 + 0.001244t + 0.000087t² - 0.0000047t³. The root mean square error of the round-trip consistency deviation decreased from 0.78 mGal to 0.15 mGal, a reduction of 81%, verifying the effectiveness of the iterative optimization method. Finally, a weighted average of the corrected gravity values ​​for the outbound and return journeys was performed using a time weighting coefficient to ensure the accurate measurement of the relative gravity value.

[0071] It should also be noted that the examples of actual values ​​for the above-mentioned parameters are merely exemplary embodiments, and the examples of actual values ​​for the above-mentioned parameters are not limited to the examples mentioned above.

[0072] In this embodiment, the outbound and return gravity measurements are performed in reverse order along the same measuring line, ensuring a symmetrical temporal distribution between the outbound and return measurement times at each measuring point. This symmetry provides a reliable time reference for determining the instantaneous drift rate of the target. A drift rate function is established based on the difference between the return and outbound gravity values ​​and the corresponding temporal relationship of the measurement times. This function accurately reflects the drift variation of the gravimeter throughout the measurement process. The cumulative drift is obtained through definite integral calculation starting from the measurement initiation time, enabling continuous accumulation calculation of the drift and avoiding discontinuities that may arise from discrete point correction. The time weighting coefficient, combined with the ratio of outbound and return measurement times, allows for a reasonable allocation of drift correction effects. Finally, the corrected target gravity value is obtained through comprehensive correction of the outbound and return gravity values, achieving precise compensation for drift errors.

[0073] The electronic device in the embodiments of this invention is described below from the perspective of hardware processing. (See attached document.) Figure 2 , Figure 2 This is a schematic diagram of the physical device structure of an electronic device in an embodiment of this application.

[0074] It should be noted that, Figure 2 The structure of the electronic device shown is merely an example and should not impose any limitation on the functionality and scope of use of the embodiments of the present invention.

[0075] like Figure 2As shown, the electronic device includes a Central Processing Unit (CPU) 201, which can perform various appropriate actions and processes according to a program stored in Read-Only Memory (ROM) 202 or a program loaded from storage portion 208 into Random Access Memory (RAM) 203, such as performing the methods described in the above embodiments. The RAM 203 also stores...

[0076] It contains various programs and data required for system operation. CPU 201, ROM 202, and RAM 203 are interconnected via bus 204. Input / output (I / O) interface 205 is also connected to bus 204.

[0077] The following components are connected to I / O interface 205: input section 206 including audio input devices, push-button switches, etc.; output section 207 including a liquid crystal display (LCD) and audio output devices, indicator lights, etc.; storage section 208 including a hard disk, etc.; and communication section 209 including a network interface card such as a LAN (Local Area Network) card, modem, etc. Communication section 209 performs communication processing via a network such as the Internet. Drive 210 is also connected to I / O interface 205 as needed. Removable media 211, such as a disk, optical disk, magneto-optical disk, semiconductor memory, etc., are installed on drive 210 as needed so that computer programs read from them can be installed into storage section 208 as needed.

[0078] In particular, according to embodiments of the present invention, the processes described above with reference to the flowcharts can be implemented as computer software programs. For example, embodiments of the present invention include a computer program product comprising a computer program carried on a computer-readable medium, the computer program containing computer programs for performing the methods shown in the flowcharts. In such embodiments, the computer program can be downloaded and installed from a network via communication section 209, and / or installed from removable medium 211. When the computer program is executed by central processing unit (CPU) 201, it performs the various functions defined in the present invention.

[0079] It should be noted that specific examples of computer-readable storage media may include, but are not limited to: electrical connections having one or more wires, portable computer disks, hard disks, random access memory (RAM), read-only memory (ROM), erasable programmable read-only memory (EPROM), flash memory, optical fiber, portable compact disc read-only memory (CD-ROM), optical storage devices, magnetic storage devices, or any suitable combination thereof. In this invention, a computer-readable storage medium can be any tangible medium containing or storing a program that can be used by or in conjunction with an instruction execution system, apparatus, or device.

[0080] The flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present invention. Each block in a flowchart or block diagram may represent a module, program segment, or portion of code, which contains one or more executable instructions for implementing a specified logical function. It should also be noted that in some alternative implementations, the functions indicated in the blocks may occur in a different order than those shown in the drawings.

[0081] Specifically, the electronic device in this embodiment includes a processor and a memory. The memory stores a computer program. When the computer program is executed by the processor, it implements the gravimeter drift correction geodetic mapping method based on the drift rate function provided in the above embodiment.

[0082] In another aspect, the present invention also provides a computer-readable storage medium, which may be included in the electronic device described in the above embodiments; or it may exist independently and not assembled into the electronic device. The storage medium carries one or more computer programs that, when executed by a processor of the electronic device, cause the electronic device to implement the gravimeter drift correction geodetic mapping method based on the drift rate function provided in the above embodiments.

[0083] The above-described embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit it. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of this application.

[0084] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. This program can be stored in a computer-readable storage medium, and when executed, it can include the processes described in the above method embodiments. The aforementioned storage medium includes various media capable of storing program code, such as ROM or random access memory (RAM), magnetic disks, or optical disks.

Claims

1. A method for drift correction of gravimeter in geodetic surveying based on drift rate function, characterized by, The method comprises the following steps: a plurality of measuring points are set on a preset measuring line, and a gravity meter is used to measure the gravity of each measuring point in a departure sequence to obtain the gravity value and the corresponding measuring time of each measuring point, and the gravity meter is used to measure the gravity of each measuring point in a return sequence opposite to the departure sequence to obtain the gravity value and the corresponding measuring time of each measuring point; a target instantaneous drift rate of each measuring point is determined according to the return gravity value, the departure gravity value, the departure measuring time and the return measuring time, a drift rate function is established by taking the time interval between the departure measuring time and the departure measuring starting time as the independent variable and the target instantaneous drift rate as the dependent variable; a first definite integral operation is performed on the drift rate function in a time advancing direction from the departure measuring starting time to obtain a target cumulative drift amount of each measuring point, which represents the integral value of the measuring point period drift amount along the measuring time axis from the departure measuring starting time to the current measuring time; a time weight coefficient of each measuring point is determined according to the departure measuring time proportion of the departure measuring time in the total departure measuring time and the return measuring time proportion of the return measuring time in the total return measuring time, and a time drift correction amount of each measuring point is determined according to the target cumulative drift amount and the time weight coefficient; the time drift correction amount is used to correct the departure gravity value and the return gravity value to obtain a target corrected gravity value; the method for correcting the departure gravity value and the return gravity value by using the time drift correction amount to obtain the target corrected gravity value further comprises: a departure temporary corrected gravity value is obtained by subtracting the departure time drift correction amount in the corresponding time drift correction amount from the departure gravity value, and a return temporary corrected gravity value is obtained by subtracting the return time drift correction amount in the corresponding time drift correction amount from the return gravity value; the departure temporary corrected gravity value and the return temporary corrected gravity value are taken as iterative variables, and the following iterative optimization operation is performed until the convergence characteristic value of the round-trip consistency deviation meets a preset convergence condition: a sixth difference calculation is performed on the departure temporary corrected gravity value and the return temporary corrected gravity value of the same measuring point to obtain the round-trip consistency deviation; the convergence characteristic value of the round-trip consistency deviation is determined according to a preset convergence index; when the convergence characteristic value does not meet the preset convergence condition, the fitting parameters of the drift rate function are adjusted according to the round-trip consistency deviation to obtain an adjusted drift rate function; a second definite integral operation is performed on the adjusted drift rate function to obtain a new target cumulative drift amount; new departure time drift correction amounts and new return time drift correction amounts are determined according to the new target cumulative drift amount and the time weight coefficient; The new in-transit time drift correction quantity is used to correct the in-transit gravity value to obtain a new in-transit temporarily corrected gravity value, and the new return time drift correction quantity is used to correct the return gravity value to obtain a new return temporarily corrected gravity value; The in-transit temporarily corrected gravity value and the return temporarily corrected gravity value that meet the preset convergence condition are weighted and averaged according to the time weight coefficient to obtain the target corrected gravity value.

2. The method of claim 1, wherein, The target instantaneous drift rate of each measuring point is determined according to the return gravity value, the in-transit gravity value, the in-transit measuring time and the return measuring time, and specifically includes: The return gravity value is subtracted from the in-transit gravity value to obtain a gravity difference value of each measuring point; A first difference calculation is performed on the return measuring time and the in-transit measuring time to obtain a round-trip time interval of each measuring point; The gravity difference value is divided by the round-trip time interval to obtain an average drift rate of each measuring point; The moving time of the gravimeter between each two adjacent measuring points in the plurality of measuring points is obtained, and the stay observation time of the gravimeter at each measuring point is obtained; The actual observation period of each measuring point is determined according to the moving time and the stay observation time; The average drift rate of each measuring point is associated with the midpoint time of the corresponding actual observation period to construct a discrete drift rate data set; An interpolation calculation is performed on the discrete drift rate data set to obtain a continuous drift rate time function; The function value of the midpoint time of the actual observation period of each measuring point is extracted from the drift rate time function, and the function value is taken as the initial instantaneous drift rate of each measuring point; The historical drift characteristic data of the gravimeter is obtained, and a drift correction coefficient is determined according to the historical drift characteristic data; The initial instantaneous drift rate is corrected by using the drift correction coefficient to obtain the target instantaneous drift rate.

3. The method of claim 2, wherein, The drift correction coefficient is determined according to the historical drift characteristic data, and specifically includes: A drift rate change curve of the gravimeter under different working durations is extracted from the historical drift characteristic data; A drift acceleration characteristic value of the gravimeter is determined according to the drift rate change curve; An accumulated working duration of a current measurement task is obtained, and a drift rate value corresponding to the accumulated working duration is determined from the drift rate change curve; The drift rate value is compared with a preset standard drift rate to obtain a drift rate ratio; The drift correction coefficient K is determined according to the drift acceleration characteristic value and the drift rate ratio through the following formula: K = 1 + (R - 1) × (1 + α × T / T0) Wherein, K is the drift correction coefficient, α is the drift acceleration characteristic value, R is the drift rate ratio, T is the accumulated working duration, and T0 is a preset standard working duration.

4. The method of claim 1, wherein, The in-transit measuring time and the in-transit measuring start time are taken as independent variables, and the target instantaneous drift rate is taken as dependent variable to establish a drift rate function, and specifically includes: performing a second difference calculation on the out-of-run measurement time and the out-of-run measurement starting time to obtain an out-of-run measurement time interval of each measurement point; determining an instantaneous drift rate difference between each two adjacent measurement points in the plurality of measurement points, and determining an out-of-run measurement time interval difference between each two adjacent measurement points in the plurality of measurement points; determining a drift acceleration according to the instantaneous drift rate difference and the out-of-run measurement time interval difference; taking the out-of-run measurement time interval as the independent variable and the target instantaneous drift rate as the dependent variable, and establishing an initial drift rate function by using a preset numerical fitting method; obtaining a working state parameter of the gravimeter, and determining a drift abnormal period of the gravimeter according to the working state parameter; extracting a drift feature of the gravimeter from the drift abnormal period, and performing weight adjustment on an abnormal function part in the initial drift rate function located in the drift abnormal period according to the drift feature to obtain the drift rate function.

5. The method of claim 1, wherein, performing a first definite integral operation on the drift rate function in a time advancing direction starting from the measurement starting time to obtain a target cumulative drift amount of each measurement point, specifically including: determining a first derivative of the drift rate function, and comparing an absolute value of the first derivative with a preset gradient threshold to obtain a gradient comparison result; dividing the drift rate function into a stable period and a fast-changing period according to the gradient comparison result; performing a first numerical integral calculation on a stable function part in the drift rate function located in the stable period by using a first integral step to obtain a stable period integral value; performing a second numerical integral calculation on a fast-changing function part in the drift rate function located in the fast-changing period by using a second integral step to obtain a fast-changing period integral value, the second integral step being smaller than the first integral step; adding the stable period integral value and the fast-changing period integral value in a time sequence of the measurement points on the time axis to obtain the target cumulative drift amount.

6. The method of claim 1, wherein, determining a time weight coefficient of each measurement point according to a proportion of the out-of-run measurement time in the total out-of-run measurement time and a proportion of the return measurement time in the total return measurement time, and determining a time drift correction amount of each measurement point according to the target cumulative drift amount and the time weight coefficient, specifically including: obtaining an out-of-run measurement ending time, and performing a third difference calculation on the out-of-run measurement ending time and the out-of-run measurement starting time to obtain the total out-of-run measurement time; obtaining a return measurement starting time and a return measurement ending time, and performing a fourth difference calculation on the return measurement ending time and the return measurement starting time to obtain the total return measurement time; dividing the out-of-run measurement time interval by the total out-of-run measurement time to obtain the proportion of the out-of-run measurement time; performing a fifth difference calculation on the return measurement ending time and the return measurement time to obtain a return measurement time interval; dividing the return measurement time interval by the total return measurement time to obtain the proportion of the return measurement time; The target accumulated drift amount is taken as a theoretical drift correction amount; A time weight coefficient is determined according to the outbound measurement time proportion and the inbound measurement time proportion; Outbound and inbound theoretical correction amounts are determined according to the theoretical drift correction amount and the time weight coefficient; Round-trip gravity values of the multiple measuring points are compared to obtain a round-trip gravity closure error; The outbound and inbound theoretical correction amounts are corrected according to the round-trip gravity closure error to obtain the time drift correction amount.

7. An electronic device, comprising: The electronic device comprises one or more processors and a memory; the memory is coupled with the one or more processors; the memory is configured to store computer program codes; the computer program codes comprise computer instructions; the one or more processors invoke the computer instructions to enable the electronic device to perform the method according to any one of claims 1-6.

8. A computer-readable storage medium comprising program instructions, characterized in that, When the program instructions run on the electronic device, the electronic device is enabled to perform the method according to any one of claims 1-6.

9. A computer program product, characterised in that, When the computer program product runs on the electronic device, the electronic device is enabled to perform the method according to any one of claims 1-6.

Images