Gravimeter drift correction geodetic surveying and mapping method and device based on drift rate function, medium and product

Abstract

The invention provides a gravimeter drift correction geodetic surveying and mapping method and device based on a drift rate function, a medium and a product, and relates to the technical field of geophysical exploration, and the method comprises the steps: setting a plurality of measuring points on a preset measuring line, and carrying out the going and returning gravity measurement through a gravimeter, respectively obtaining an outbound gravity value, a return gravity value and a corresponding measurement moment of each measurement point; determining a target instantaneous drift rate of each measuring point according to a time difference between the forward and backward gravity values, and establishing a drift rate function by taking the measuring time interval as an independent variable and the target instantaneous drift rate as a dependent variable; performing definite integral operation on the drift rate function to obtain a target accumulated drift distance of each measuring point; determining a time weight coefficient by combining the measurement time proportion of the forward stroke and the backward stroke, and multiplying the target accumulated drift amount by the time weight coefficient to obtain a time drift correction amount; and performing drift correction on the forward and backward gravity values by using the time drift correction value to obtain a target correction gravity value.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of geophysical exploration, and particularly relates to a gravity meter drift correction geodetic surveying method, device, medium and product based on a drift rate function. BACKGROUND

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

[0003] In the related art, a correction method based on a linear time drift assumption is usually used. In specific implementation, a surveyor first performs initial gravity measurement at the starting point of a survey line using a gravity meter, records the initial gravity value and the measurement time; then performs single-pass gravity measurement on each survey point along the survey line in turn to obtain the observed gravity value of each survey point; finally, returns to the starting point of the survey line to perform repeated measurement, calculates the average drift rate according to the gravity value difference and the time interval of the two measurements at the starting point, and calculates the drift correction value of each survey point according to the time difference between the measurement time of each survey point and the initial time by using linear interpolation. Although the basic drift correction can be completed, the entire process is heavily dependent on the linear drift assumption.

[0004] However, the above drift correction method estimates the unified drift rate of the entire survey line only according to the two measurement results at the starting point, which is based on the ideal assumption that the instrument drift strictly follows the linear law. In fact, the drift of the gravity meter is often nonlinear due to the influence of temperature change and material fatigue of the quartz spring system of the gravity meter. However, the drift correction method ignores the difference in actual drift of each survey point, resulting in inherent deviation between the linear model and the real drift curve. More importantly, this deviation will be systematically accumulated with the passage of time, so that the farther the survey point is from the starting point, the greater the correction error will be, thereby leading to low correction accuracy of the gravity meter measurement drift in the related art. SUMMARY

[0005] The present application provides a gravity meter drift correction geodetic surveying method, device, medium and product based on a drift rate function, which is used to improve the correction accuracy of the gravity meter measurement drift.

[0006] In a first aspect, the present application provides a drift correction method for gravimeter drift based on a drift rate function, applied to the electronic device, the method comprising: setting a plurality of measuring points on a preset measuring line, measuring the out-bound gravity of each measuring point in an out-bound order using the gravimeter to obtain the out-bound gravity value and the corresponding out-bound measuring time of each measuring point, and measuring the return-bound gravity of each measuring point in a return-bound order opposite to the out-bound order using the gravimeter to obtain the return-bound gravity value and the corresponding return-bound measuring time of each measuring point; determining the target instantaneous drift rate of each measuring point according to the return-bound gravity value, the out-bound gravity value, the out-bound measuring time and the return-bound measuring time, and establishing a drift rate function with the time interval between the out-bound measuring time and the out-bound measuring starting time as the independent variable and the target instantaneous drift rate as the dependent variable; performing a first definite integral operation on the drift rate function in the time advancing direction from the measuring starting time to obtain the target cumulative drift of each measuring point, the target cumulative drift representing the integral value of the measuring point period drift along the measuring time axis from the measuring starting time to the current measuring time; determining the time weight coefficient of each measuring point according to the out-bound measuring time proportion of the out-bound measuring time in the total out-bound measuring time and the return-bound measuring time proportion of the return-bound measuring time in the total return-bound measuring time, and determining the time drift correction value of each measuring point according to the target cumulative drift and the time weight coefficient; and performing drift correction on the out-bound gravity value and the return-bound gravity value using the time drift correction value to obtain the target corrected gravity value.

[0007] By adopting the above technical solution, the out-bound gravity measurement and the return-bound gravity measurement are performed in opposite orders on the same measuring line, so that the out-bound measuring time and the return-bound measuring time of each measuring point can form a symmetrical distribution relationship in time, and this symmetry can provide a reliable time reference for subsequent determination of the target instantaneous drift rate. The drift rate function is established based on the difference between the return-bound gravity value and the out-bound gravity value and the time relationship of the corresponding measuring time, so that the function can accurately reflect the drift change rule of the gravimeter in the entire measurement process. The target cumulative drift is obtained by the definite integral operation from the measuring starting time, which can realize continuous cumulative calculation of the drift, and avoid the discontinuity problem caused by discrete point correction. The time weight coefficient combines the measurement time proportion of the out-bound and return-bound, so that the time drift correction value can reasonably allocate the drift influence, and finally the target corrected gravity value is obtained by comprehensive correction of the out-bound gravity value and the return-bound gravity value, which can realize accurate compensation of the drift error. Further, the technical problem of low correction accuracy of gravimeter measurement drift in related technologies is solved, and the technical effect of improving the correction accuracy of gravimeter measurement drift is achieved.

[0008] Optionally, the target instantaneous drift rate of each measuring point is determined according to the return gravity value, the outbound gravity value, the outbound measuring time and the return measuring time, specifically including: subtracting the return gravity value from the outbound gravity value to obtain a gravity difference value of each measuring point; performing a first difference calculation on the return measuring time and the outbound measuring time to obtain a round-trip time interval of each measuring point; dividing the gravity difference value by the round-trip time interval to obtain an average drift rate of each measuring point; obtaining a moving time of the gravimeter between each two adjacent measuring points in the plurality of measuring points, and obtaining a stay observation time of the gravimeter at each measuring point; determining an actual observation period of each measuring point according to the moving time and the stay observation time; associating the average drift rate of each measuring point with a midpoint time of the corresponding actual observation period to construct a discrete drift rate data set; performing interpolation calculation on the discrete drift rate data set to obtain a continuous drift rate time function; extracting a function value of the midpoint time of the actual observation period of each measuring point from the drift rate time function, and taking the function value as an initial instantaneous drift rate of each measuring point; obtaining historical drift characteristic data of the gravimeter, and determining a drift correction coefficient according to the historical drift characteristic data; correcting the initial instantaneous drift rate by using the drift correction coefficient to obtain the target instantaneous drift rate.

[0009] By adopting the above technical solution, the difference between the return gravity value and the outbound gravity value is divided by the round-trip time interval to obtain the average drift rate, which can provide a basic drift estimation value for each measuring point. The combination of the moving time and the stay observation time can determine the actual observation period, so that the average drift rate can be accurately associated with the midpoint time of the corresponding observation period to form the discrete drift rate data set. The discrete data is converted into the continuous drift rate time function through interpolation calculation, which can eliminate the sudden change of the drift rate between the measuring points and ensure the continuity of the drift change. The initial instantaneous drift rate extracted from the continuous function can reflect the real-time drift characteristics of each measuring point, and the drift correction coefficient determined according to the historical drift characteristic data is used to correct the initial instantaneous drift rate, which can fully consider the influence of the historical working state of the gravimeter on the current drift, so that the target instantaneous drift rate can be more consistent with the actual drift characteristics of the gravimeter, and the accuracy of the drift rate estimation is improved.

[0010] Optionally, the drift correction coefficient is determined according to the historical drift characteristic data, specifically including: extracting a drift rate change curve of the gravimeter under different working durations from the historical drift characteristic data; determining a drift acceleration characteristic value of the gravimeter according to the drift rate change curve; obtaining an accumulated working duration of the current measurement task, and determining a drift rate value corresponding to the accumulated working duration from the drift rate change curve; comparing the drift rate value with a preset standard drift rate to obtain a drift rate ratio; determining the drift correction coefficient K 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 a drift correction coefficient, a is a drift acceleration characteristic value, R is a drift rate ratio, T is a cumulative working time length, and T0 is a preset standard working time length.

[0011] By adopting the above technical solution, the drift rate change curve extracted from the historical drift characteristic data can reflect the drift evolution law of the gravimeter under different working time lengths, and the drift acceleration characteristic value can quantify the change trend of the drift rate. The ratio of the drift rate value corresponding to the cumulative working time length to the preset standard drift rate can reflect the deviation degree of the current drift state relative to 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 a and the time ratio T / T0 can reflect the time-varying characteristics of the drift. The two interact through the form of (R-1)×(1+α×T / T0), so that the drift correction coefficient K can dynamically adapt to the drift characteristic changes of the gravimeter in different working stages. This correction method which comprehensively considers the basic deviation and time-varying characteristics can ensure that the drift correction coefficient accurately reflects the actual working state of the gravimeter, and improves the adaptability and accuracy of drift correction.

[0012] Optionally, a drift rate function is established with the time interval between the out-of-travel measurement time and the out-of-travel measurement starting time as the independent variable and the target instantaneous drift rate as the dependent variable, specifically including: performing a second difference calculation on the out-of-travel measurement time and the out-of-travel measurement starting time to obtain the out-of-travel measurement time interval of each measurement point; determining the instantaneous drift rate difference between each two adjacent measurement points in the plurality of measurement points, and determining the out-of-travel measurement time interval difference between each two adjacent measurement points in the plurality of measurement points; determining the drift acceleration according to the instantaneous drift rate difference and the out-of-travel measurement time interval difference; establishing an initial drift rate function by using a preset numerical fitting method with the out-of-travel measurement time interval as the independent variable and the target instantaneous drift rate as the dependent variable; obtaining the working state parameters of the gravimeter, and determining the drift abnormal period of the gravimeter according to the working state parameters; extracting the drift characteristics of the gravimeter from the drift abnormal period, and adjusting the weight of the abnormal function part in the initial drift rate function located in the drift abnormal period according to the drift characteristics to obtain the drift rate function.

[0013] By adopting the technical solution, the time interval of the approach measurement can be taken as an independent variable to establish a unified time reference, and the ratio of the difference in the instantaneous drift rate to the difference in the time interval can determine the drift acceleration and reflect the change speed of the drift rate. The preset numerical fitting mode establishes an initial drift rate function based on the time interval and the target instantaneous drift rate, and can realize the conversion of discrete data to continuous functions. The drift anomaly period identified by the working state parameter identification can mark the abnormal working period of the gravimeter, and the drift features extracted from the anomaly period can quantify the abnormality degree. By adjusting the weight of the abnormal 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 rule of the gravimeter, while not completely ignoring the information of the anomaly period, realizing the balanced processing of normal drift and abnormal drift, and improving the robustness of the function.

[0014] Optionally, the first numerical integral calculation is performed on the drift rate function in the time advancing direction from the measurement starting moment, to obtain the target cumulative drift amount of each measurement point, specifically including: 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 fast-changing period according to the gradient comparison result; performing first numerical integral calculation on the stable function part in the stable period of the drift rate function by using a first integral step, to obtain a stable period integral value; performing second numerical integral calculation on the fast-changing function part in the fast-changing period of the drift rate function 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; and accumulating the stable period integral value and the fast-changing period integral value in the time sequence of each measurement point on the time axis to obtain the target cumulative drift amount.

[0015] By adopting the technical solution, the absolute value of the first derivative of the drift rate function can reflect the change intensity of the drift rate, and the comparison with the preset gradient threshold can realize the quantitative evaluation of the change characteristics of the function. Based on the gradient comparison result, the function is divided into a stable period and a fast-changing period, so that the function parts with different change characteristics can adopt the corresponding integral strategy. The numerical integral of the first integral step on the stable function part can ensure the calculation efficiency, and the more precise integral calculation of the second integral step on the fast-changing function part can ensure the integral accuracy, and the differential setting of the two steps can realize the balance between the calculation efficiency and the accuracy. The stable period integral value and the fast-changing period integral value are accumulated in the time sequence, which can ensure the time sequence continuity of the cumulative drift amount, so that the target cumulative drift amount can accurately reflect the complete drift accumulation process from the measurement starting moment to the current moment, realizing the efficient integral calculation with adaptive accuracy.

[0016] Optionally, the time weight coefficient of each measuring point is determined according to the time proportion of the outbound measurement at the outbound measurement time in the total outbound measurement time and the time proportion of the return measurement at the return measurement time in the total return measurement time, and the time drift correction of each measuring point is determined according to the target cumulative drift and the time weight coefficient, specifically including: obtaining the outbound measurement end time, and performing third difference calculation on 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 the return measurement end time, and performing fourth difference calculation on the return measurement end time and the return measurement start time to obtain the total return measurement time; dividing the outbound measurement time interval by the total outbound measurement time to obtain the outbound measurement time proportion; performing fifth difference calculation on the return measurement end time and the return measurement time to obtain the return measurement time interval; dividing the return measurement time interval by the total return measurement time to obtain the return measurement time proportion; taking the target cumulative drift as the theoretical drift correction; determining the time weight coefficient according to the outbound measurement time proportion and the return measurement time proportion; determining the outbound theoretical correction and the return theoretical correction according to the theoretical drift correction and the time weight coefficient; comparing the round-trip gravity values of the plurality of measuring points to obtain the round-trip gravity closure error; and modifying the outbound theoretical correction and the return theoretical correction according to the round-trip gravity closure error to obtain the time drift correction.

[0017] By adopting the above technical solution, the total outbound measurement time and the total return measurement time can provide the basis for time normalization, the ratio of the outbound measurement time interval to the total time obtains the outbound measurement time proportion, the ratio of the return measurement time interval to the total time obtains the return measurement time proportion, and the two proportion values can jointly determine the calculation basis of the time weight coefficient. The target cumulative drift as the theoretical drift correction can reflect the influence of the theoretical drift based on the drift rate function model, and the round-trip gravity closure error can quantify the inconsistency of the round-trip measurement in the actual measurement. The time weight coefficient determines the weight distribution of the outbound and return according to the relative position of the measuring point in the measurement time sequence, and the outbound theoretical correction and the return theoretical correction calculated based on the theoretical drift correction and the time weight coefficient reflect the prediction of the theoretical model on the drift, while the round-trip gravity closure error reflects the deviation between the theoretical model and the actual measurement. By modifying the outbound theoretical correction and the return theoretical correction according to the round-trip gravity closure error, the theoretical model prediction and the actual measurement feedback can be combined, which can not only ensure the theoretical continuity and time consistency of the drift correction, but also eliminate the systematic deviation of the theoretical model through the feedback correction of the actual measurement closure error, realize the cooperative optimization of the theoretical correction and the actual measurement correction, avoid the problem of excessive correction caused by repeated correction, and improve the accuracy and reliability of the drift correction.

[0018] Optionally, the time drift correction amount is used to correct the outbound gravity value and the return gravity value to obtain a target corrected gravity value, and the specific method further comprises: subtracting the outbound time drift correction amount in the corresponding time drift correction amount from the outbound gravity value to obtain an outbound temporary corrected gravity value, and subtracting the return time drift correction amount in the corresponding time drift correction amount from the return gravity value to obtain a return temporary corrected gravity value; taking the outbound temporary corrected gravity value and the return temporary corrected gravity value as iteration variables, and performing the following iteration optimization operation until the convergence eigenvalue of the round-trip consistency deviation meets the preset convergence condition: performing sixth difference calculation on the outbound temporary corrected gravity value and the return temporary corrected gravity value of the same measurement point to obtain the round-trip consistency deviation; determining the convergence eigenvalue of the round-trip consistency deviation according to the preset convergence index; when the convergence eigenvalue does not meet the preset convergence condition, adjusting the fitting parameters of the drift rate function according to the round-trip consistency deviation to obtain an adjusted drift rate function; performing second definite integral operation on the adjusted drift rate function to obtain a new target cumulative drift amount; determining new outbound time drift correction amount and new return time drift correction amount according to the new target cumulative drift amount and the time weight coefficient; using the new outbound time drift correction amount to correct the outbound gravity value to obtain a new outbound temporary corrected gravity value, and using the new return time drift correction amount to correct the return gravity value to obtain a new return temporary corrected gravity value; and performing weighted average on the outbound temporary corrected gravity value and the return temporary corrected gravity value that meet the preset convergence condition according to the time weight coefficient to obtain the target corrected gravity value.

[0019] By adopting the above technical solution, the difference between the outbound temporary corrected gravity value and the return temporary corrected gravity value forms the round-trip consistency deviation, which can directly reflect the pros and cons of the current correction effect. The preset convergence index converts the deviation into the convergence eigenvalue, which can provide a quantitative criterion for iteration termination. When the convergence condition is not met, the round-trip consistency deviation can guide the adjustment direction and amplitude of the fitting parameters of the drift rate function, and the adjusted drift rate function generates a new target cumulative drift amount through the second definite integral operation, which can then update the time drift correction amount. The new correction amount re-corrects the original gravity value to form a new temporary corrected value, forming a complete iterative closed loop. This iteration optimization mechanism based on deviation feedback enables the correction process to be self-adaptive, until the round-trip consistency meets the preset requirement. Finally, the weighted average of the round-trip corrected values that meet the convergence condition according to the time weight coefficient can fully utilize the round-trip measurement information to obtain the optimal target corrected gravity value.

[0020] In a second aspect, an electronic device is provided, which includes one or more processors and a memory; the memory is coupled to the one or more processors, and is configured to store computer program codes including computer instructions; the 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 manner of the first aspect.

[0021] In a third aspect, a computer program product including instructions is provided, which, when executed on an electronic device, causes the electronic device to perform the method described in the first aspect and any possible implementation manner of the first aspect.

[0022] In a fourth aspect, a computer-readable storage medium is provided, which includes instructions, which, when executed on an electronic device, causes the electronic device to perform the method described in the first aspect and any possible implementation manner of the first aspect. BRIEF DESCRIPTION OF DRAWINGS

[0023] Figure 1 is a flowchart of a gravity meter drift correction geodetic survey method based on a drift rate function in embodiments of the present application; Figure 2 is a schematic diagram of an entity device structure of an electronic device in embodiments of the present application. DETAILED DESCRIPTION

[0024] The terms used in the following embodiments of the present application are only for the purpose of describing specific embodiments and are not intended to be limiting of the present application. As used in the specification and the appended claims of the present application, the singular forms "a," "an" and "the" are intended to include the plural forms as well, unless the context clearly indicates otherwise. It will be further understood that the terms "and / or," as used in the present application, refers to any or all possible combinations of one or more of the associated listed items.

[0025] Hereinafter, the terms "first" and "second" are only for the purpose of description, and cannot be understood as implying or suggesting relative importance or implicitly indicating the number of the indicated technical features. Therefore, the features defined with "first" and "second" can explicitly or implicitly include one or more of the features, and in the description of the embodiments of the present application, the meaning of "a plurality of" is two or more, unless otherwise specified.

[0026] The present application provides a gravity meter drift correction geodetic survey method based on a drift rate function, as shown in Figure 1 , Figure 1is a flowchart of a gravity meter drift correction geodetic surveying method based on a drift rate function in the embodiment of the present application, comprising the following steps: Step S101, a plurality of survey points are set on a preset survey line, the gravity meter is used to measure the gravity of each survey point in a departure sequence to obtain the departure gravity value and the corresponding departure measurement time of each survey point, and the gravity meter is used to measure the gravity of each survey point in a return sequence opposite to the departure sequence to obtain the return gravity value and the corresponding return measurement time of each survey point; Step S102, the target instantaneous drift rate of each survey point is determined according to the return gravity value, the departure gravity value, the departure measurement time and the return measurement time, and the drift rate function is established with the time interval between the departure measurement time and the departure measurement starting time as the independent variable and the target instantaneous drift rate as the dependent variable; Step S103, the drift rate function is subjected to a first definite integral operation in the time advancing direction from the measurement starting time to obtain the target cumulative drift of each survey point, and the target cumulative drift represents the integral value of the survey point period drift along the measurement time axis from the measurement starting time to the current measurement time; Step S104, the time weight coefficient of each survey point is determined according to the departure measurement time proportion of the departure measurement time in the total departure measurement time and the return measurement time proportion of the return measurement time in the total return measurement time, and the time drift correction amount of each survey point is determined according to the target cumulative drift and the time weight coefficient; Step S105, the departure gravity value and the return gravity value are subjected to drift correction by using the time drift correction amount to obtain the target corrected gravity value.

[0027] In the above embodiment, a CG-6 relative gravimeter is used for gravity measurement in a certain mining area, and the reading resolution of the gravimeter reaches 0.1 mGal (milligal, 1 mGal = 10-5m / s2). A preset survey line is arranged along the main exploration line of the mining area, with a total length of 12 kilometers and 25 survey points arranged at intervals of about 500 meters. Each survey point is marked with a concrete monument and precisely positioned by GNSS (Global Navigation Satellite System). The outbound measurement starts at 8:00 on March 15, and the gravimeter starts from the gravity secondary base point A and reaches survey points P1, P2...P25 in turn. At each survey point, the gravimeter automatically performs 120 seconds of continuous observation, records 60 readings, automatically eliminates abnormal values, and then takes the average to obtain the relative gravity reading of the survey point. For example, the outbound gravity value of survey point P1 is 2834.5 mGal (relative value with respect to the internal reference of the instrument), and the outbound measurement time is 8:15:30; the outbound gravity value of survey point P25 is 2891.2 mGal, and the outbound measurement time is 14:30:45. The return measurement starts at 15:00 on the same day, and the measurement is performed in the reverse order of P25, P24...P1, A. The return gravity value of survey point P25 is 2891.3 mGal, and the return measurement time is 15:00:30; the return gravity value of survey point P1 is 2834.6 mGal, and the return measurement time is 21:15:20.

[0028] In the above embodiment, the target instantaneous drift rate is obtained by calculating the data of the round trip measurement. Taking survey point P1 as an example, the difference between the return gravity value 2834.6 mGal and the outbound gravity value 2834.5 mGal is 0.1 mGal, and the round trip time interval is 12 hours 59 minutes and 50 seconds (46790 seconds), and the preliminary drift rate of this point is calculated to be 0.01 mGal / hour. Considering the moving time of the gravimeter between survey points (average 15 minutes) and the actual observation period of each survey point (actual observation period = arrival time at survey point + observation time), a continuous drift rate time function is established by cubic spline interpolation. The drift rate function can be fitted by a second-order polynomial: v(t) = a0 + a1 x t + a2 x t2, where t is the time interval with respect to the start time of the measurement (unit: hour), the fitting coefficients a0 = 0.01 mGal / h, a1 = 0.001 mGal / h2, and a2 = -0.0002 mGal / h3. The drift rate function reflects the nonlinear variation characteristics of the drift rate of the gravimeter with time, and the initial drift rate is small, which gradually increases with the working time, which is consistent with the typical drift characteristics of the CG-6 gravimeter. The target cumulative drift amount is calculated by definite integral of the drift rate function: D(t) = ∫0 tv(τ)dτ = a0x t + (a1 / 2)x t2 + (a2 / 3)x t3. For example, the cumulative drift of measuring point P10 is 0.04 mGal at the departure measurement time (3.5 hours relative to the starting time), and 0.13 mGal at the return measurement time (10.2 hours relative to the starting time).

[0029] In the above embodiment, the time weight coefficient is determined according to the measurement time ratio. The total duration of the departure measurement is 6 hours 30 minutes 45 seconds, and the total duration of the return measurement is 6 hours 14 minutes 50 seconds. For measuring point P10, the departure 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 weight coefficient is linearly interpolated: the departure weight is 0.45, and the return weight is 0.55, ensuring the reasonable fusion of the data of the departure and return measurements. The time drift correction amount includes two parts: the main body drift correction amount directly uses the target cumulative drift amount; and the closure correction amount is distributed according to the time ratio through the departure and return gravity closure difference (the departure and return difference of base point A is 0.1 mGal). The departure time drift correction amount of measuring point P10 is 0.04 + 0.1x 0.54 = 0.09 mGal, and the return time drift correction amount is 0.13 + 0.1x 0.41 = 0.17 mGal. Finally, the corrected gravity value of measuring point P10 in the departure is 2856.4-0.09 = 2856.31 mGal, and the corrected gravity value in the return is 2856.5-0.17 = 2856.33 mGal. The target corrected gravity value is obtained by weighted average: 2856.31x 0.45 + 2856.33x 0.55 = 2856.32 mGal. This value represents the gravity difference of measuring point P10 relative to gravity secondary base point A, eliminates the influence of instrument drift, and the measurement accuracy reaches ±0.08 mGal, meeting the accuracy requirements of mineral exploration gravity measurement. In actual application, the required measurement accuracy is adjusted according to different measurement task requirements: high-precision mode (0.06-0.09 mGal), suitable for fine geological structure detection and oil and gas exploration; standard precision mode (0.09-0.15 mGal), suitable for general mineral exploration and regional gravity survey; and fast measurement mode (0.15 mGal and above), suitable for preliminary exploration and large-scale gravity survey. In this embodiment, the high-precision mode is taken as an example, and the measurement accuracy is improved by prolonging the observation time and increasing the number of readings.

[0030] By the above steps, the outbound gravity measurement and the return gravity measurement are performed in the same measurement line in reverse order, so that the outbound measurement time and the return measurement time of each measurement point can form a symmetrical distribution relationship in time. This symmetry can provide a reliable time reference for subsequent determination of the target instantaneous drift rate. A drift rate function is established based on the difference between the return gravity value and the outbound gravity value and the time relationship of the corresponding measurement time, so that the function can accurately reflect the drift change rule of the gravimeter during the entire measurement process. The target cumulative drift amount is obtained by a definite integral operation from the start time of the measurement, which can realize continuous cumulative calculation of the drift amount and avoid the discontinuity problem caused by discrete point correction. The time weight coefficient combines the measurement time proportion of the outbound and return, so that the time drift correction amount can reasonably distribute the drift influence. Finally, the target corrected gravity value is obtained by comprehensive correction of the outbound gravity value and the return gravity value, which can realize accurate compensation of the drift error. Further, the technical problem of low correction accuracy of the gravimeter measurement drift in the related art is solved, and the technical effect of improving the correction accuracy of the gravimeter measurement drift is achieved.

[0031] The execution subject of the above steps can 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 separate controller or processor, or other processing devices or processing units with similar processing functions, but is not limited thereto.

[0032] In an optional embodiment, the target instantaneous drift rate of each measurement point is determined according to the return gravity value, the outbound gravity value, the outbound measurement time and the return measurement time, specifically including: subtracting the return gravity value from the outbound gravity value to obtain the gravity difference value of each measurement point; performing a first difference calculation on the return measurement time and the outbound measurement time to obtain the round-trip time interval of each measurement point; dividing the gravity difference value by the round-trip time interval to obtain the average drift rate of each measurement point; obtaining the moving time of the gravimeter between each two adjacent measurement points in the plurality of measurement points, and obtaining the stay observation time of the gravimeter at each measurement point; determining the actual observation period of each measurement point according to the moving time and the stay observation time; associating the average drift rate of each measurement point with the midpoint time of the corresponding actual observation period to construct a discrete drift rate data set; performing interpolation calculation on the discrete drift rate data set to obtain a continuous drift rate time function; extracting the function value of the midpoint time of the actual observation period of each measurement point from the drift rate time function, and taking the function value as the initial instantaneous drift rate of each measurement point; obtaining the historical drift characteristic data of the gravimeter, and determining a drift correction coefficient according to the historical drift characteristic data; correcting the initial instantaneous drift rate by using the drift correction coefficient to obtain the target instantaneous drift rate.

[0033] In the above embodiment, a CG-6 relative gravimeter is used for detailed gravity measurement in a certain iron mine area. First, the gravity difference of each measuring point is calculated. Take measuring point P5 as an example. The return gravity value of this point is 2847.9 mGal, and the outbound gravity value is 2847.8 mGal. The difference between the two is 0.1 mGal. This gravity difference reflects the total drift of the instrument during the round-trip measurement, which is the relative change value relative to the internal quartz spring system reference of the instrument. The round-trip time interval is obtained by accurate time recording calculation. The outbound measurement time of measuring point P5 is 9:42:15, and the return measurement time is 19:28:40. The time interval is 9 hours 46 minutes and 25 seconds (35185 seconds). Divide the gravity difference 0.1 mGal by the time interval 9.773 hours to obtain the average drift rate of this point, which is 0.0102 mGal / hour. This average drift rate represents the average drift characteristics of this measuring point during the round-trip measurement. The moving time of the gravimeter between adjacent measuring points is accurately obtained by GPS (Global Positioning System) trajectory recording. For example, the moving time from measuring point P4 to P5 is 12 minutes and 30 seconds, and the moving time from P5 to P6 is 13 minutes and 45 seconds. The stay observation time is the actual measurement time of the instrument at each measuring point. The CG-6 gravimeter stays at each measuring point for 120 seconds for automatic observation. The actual observation period refers to the complete period from setting the instrument when arriving at the measuring point to completing the observation and leaving. The actual observation period of measuring point P5 is from 9:40:00 to 9:43:00, and the midpoint time is 9:41:30.

[0034] In the above embodiment, the construction of the discrete drift rate dataset is achieved by pairing the average drift rate of the 25 measurement points with their corresponding midpoint time in the actual observation period. For example, the midpoint time of 8:16:00 for measurement point P1 corresponds to a drift rate of 0.0094 mGal / h, the midpoint time of 9:41:30 for measurement point P5 corresponds to a drift rate of 0.0102 mGal / h, and the midpoint time of 11:23:00 for measurement point P10 corresponds to a drift rate of 0.0111 mGal / h. These discrete data points reflect the trend of drift rate change over time. A cubic spline interpolation method 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 × t2- 0.0000032 × t3, where t is the time relative to the start of the measurement (unit: hours). The function value at the midpoint time of the actual observation period for each measurement point is extracted from this drift rate time function. The function value at t = 1.692 hours for measurement point P5 is 0.0100 mGal / h, which is used as the initial instantaneous drift rate for this point. The historical drift characteristic data is obtained from the long-term monitoring records of the gravimeter. The historical data of the CG-6 gravimeter shows that within the first 4 hours of continuous operation, the drift rate is relatively stable, about 0.008-0.011 mGal / h; after 4-8 hours of operation, due to the temperature effect and material fatigue of the quartz spring, the drift rate increases to 0.011-0.016 mGal / h; and after more than 8 hours, the drift rate can reach 0.016-0.022 mGal / h. Based on these historical characteristics, the drift correction coefficient is determined.

[0035] In the above embodiment, the determination of the drift correction coefficient takes into account the cumulative working time of the instrument and environmental factors. In this measurement, the measurement time of measurement point P5 is 2.5 hours away from the start of the instrument, and according to the drift characteristics under the same working time in the built-in database, the correction coefficient K = 1.02 is determined. This means that the initial instantaneous drift rate needs to be corrected upward by 2%. Therefore, the target instantaneous drift rate of measurement point P5 is 0.0100 × 1.02 = 0.0102 mGal / h. Through the above refinement, each measurement point obtains a target instantaneous drift rate that takes into account the characteristics of the instrument and the time effect. For example, the target instantaneous drift rate of early measurement point P1 is 0.0096 mGal / h, the target instantaneous drift rate of middle measurement point P12 is 0.0118 mGal / h, and the target instantaneous drift rate of late measurement point P25 is 0.0142 mGal / h, showing an increasing trend consistent with the drift characteristics of the CG-6 gravimeter.

[0036] In an optional embodiment, the drift correction coefficient is determined according to the historical drift characteristic data, specifically comprising: extracting the drift rate change curve of the gravimeter under different working durations from the historical drift characteristic data; determining the drift acceleration characteristic value of the gravimeter according to the drift rate change curve; 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 change curve; comparing the drift rate value with the preset standard drift rate to obtain a drift rate ratio; and determining the drift correction coefficient K according to the drift acceleration characteristic value and the drift rate ratio through the following formula: K = 1 + (R - 1) x (1 + a x T / T0) wherein K is the drift correction coefficient, a is the drift acceleration characteristic value, R is the drift rate ratio, T is the cumulative working duration, and T0 is the preset standard working duration.

[0037] In the above embodiment, a CG-6 relative gravimeter is used for mine area gravity measurement, which is equipped with an internal database to continuously record historical measurement data for 180 days. The historical drift characteristic data refers to the data set of the drift rate of the gravimeter changing with time recorded in the past measurement tasks, including the drift performance under different temperature conditions, different working durations and different measurement environments. The drift rate change curve of the gravimeter under different working durations is extracted from the internal database. The curve shows that within 0-2 hours of work, the drift rate remains at 0.008-0.010 mGal / h; within 2-4 hours of work, the drift rate increases to 0.010-0.012 mGal / h; within 4-6 hours of work, the drift rate reaches 0.012-0.015 mGal / h; within 6-8 hours of work, the drift rate rises to 0.015-0.018 mGal / h; and after more than 8 hours of work, the drift rate can reach 0.018-0.022 mGal / h. This increasing trend is mainly caused by the temperature response and material elastic change of the quartz spring system. The drift acceleration characteristic value a is obtained by taking the second derivative of the drift rate change curve. In the specific calculation, the drift rate values of adjacent time points on the curve are selected to calculate the rate of change. For example, within the 2-4 hour working duration interval, the drift rate increases from 0.010 mGal / h to 0.012 mGal / h, with a change rate of 0.001 mGal / h²; within the 4-6 hour interval, the change rate is 0.0015 mGal / h²; and within the 6-8 hour interval, the change rate is 0.0015 mGal / h². The average of these change rates is taken and normalized to obtain the drift acceleration characteristic value a = 0.00133.

[0038] In the above embodiment, the cumulative working time length of the current measurement task refers to the total time length from the start of the gravimeter to the current measurement time of the measurement point. Taking measurement point P15 as an example, the travel time of this point is 12:45:30, which has been 5.25 hours since the start of the morning 7:30. From the historical drift rate change curve, the drift rate value corresponding to 5.25 hours is determined by linear interpolation to be 0.0135 mGal / h. The preset standard drift rate is a reference value determined according to the technical specifications of the CG-6 gravimeter, and is set to 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 0.0135 mGal / h with the standard value 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), and the current cumulative working time T = 5.25 hours, substituting the drift correction coefficient formula K = 1 + (R-1) x (1+αxT / T0) = 1 + (1.2273-1) x (1+0.00133x5.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 of the instrument after long time working. In practical application, the drift correction coefficient of different measurement points will be dynamically adjusted according to the cumulative working time length of its measurement time. The correction coefficient of early measurement points such as P1 (cumulative working time 0.5 hours) is K = 1.018, the correction coefficient of middle measurement points such as P12 (cumulative working time 3.5 hours) is K = 1.195, and the correction coefficient of late measurement points such as P25 (cumulative working time 7 hours) is K = 1.312. Through this dynamic correction mechanism, each measurement point can obtain a drift correction coefficient matched with its measurement time, ensuring the time-varying adaptability of drift correction.

[0039] In an optional embodiment, a drift rate function is established with the time interval between the time of the transit measurement and the starting time of the transit measurement as the independent variable and the target instantaneous drift rate as the dependent variable, specifically including: performing a second difference calculation on the time of the transit measurement and the starting time of the transit measurement to obtain the transit measurement time interval of each measuring point; determining the instantaneous drift rate difference between each two adjacent measuring points in the plurality of measuring points, and determining the transit measurement time interval difference between each two adjacent measuring points in the plurality of measuring points; determining the drift acceleration according to the instantaneous drift rate difference and the transit measurement time interval difference; establishing an initial drift rate function by using a preset numerical fitting method with the transit measurement time interval as the independent variable and the target instantaneous drift rate as the dependent variable; obtaining the working state parameter of the gravimeter, and determining the drift abnormal period of the gravimeter according to the working state parameter; extracting the drift feature of the gravimeter from the drift abnormal period, and performing weight adjustment on the 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.

[0040] In the above embodiment, the CG-6 relative gravimeter is used for copper mine area gravity measurement. First, the transit measurement time interval of each measuring point is calculated. Take measuring point P8 as an example. The time of the transit measurement of this point is 10:45:20, the starting time of the transit measurement is 8:00:00, and the difference between the two is 2 hours 45 minutes and 20 seconds, i.e. 2.756 hours. This time interval reflects the length of time experienced from the start of the measurement to the transit measurement of the measuring point. The determination of the instantaneous drift rate difference is achieved by comparing the target instantaneous drift rates of adjacent measuring points. For example, the target instantaneous drift rate of measuring point P7 is 0.0104 mGal / h, and that of measuring point P8 is 0.0106 mGal / h, with a difference of 0.0002 mGal / h. Correspondingly, the transit measurement time interval of P7 is 2.423 hours, that of P8 is 2.756 hours, and the time interval difference is 0.333 hours. The drift acceleration is calculated by dividing the instantaneous drift rate difference by the time interval difference: 0.0002 / 0.333=0.0006 mGal / h². The preset numerical fitting method can use the least square polynomial fitting. The transit measurement time interval of the 25 measuring points is taken as the independent variable t, and the target instantaneous drift rate is taken as the dependent variable v, and a cubic polynomial form of the initial drift rate function is fitted: v(t)=0.0085+0.0013t+0.000084t²-0.0000049t³. The function can describe the nonlinear variation law of the drift rate with time, wherein the constant term 0.0085 mGal / h represents the initial drift rate, the linear term coefficient 0.0013 reflects the linear growth trend, and the quadratic and cubic term coefficients describe the nonlinear characteristics.

[0041] In the above embodiment, the working state parameters include the internal temperature of the instrument, the tilt compensation state, the battery voltage, and other key indicators. The internal temperature sensor of the CG-6 gravimeter monitored that, during the 4.5th to 5.2nd hour of measurement, due to the rapid change in ambient temperature (from 18°C to 26°C), the internal temperature control system of the instrument entered the adjustment state, and the temperature fluctuation reached ±0.3°C. The tilt compensation system record showed that, at the 3.8th hour, the X-axis tilt angle changed instantaneously by more than 15 arcseconds due to ground vibration. The battery voltage monitoring showed that it remained stable in the range of 12.4-12.6V throughout the measurement, and no abnormalities occurred. The abnormal period of drift was determined according to the abnormal values of the working state parameters. The temperature adjustment period (4.5-5.2 hours) was marked as the first abnormal period, and the drift rate during this period showed obvious volatility, deviating from the normal trend by 0.0022 mGal / h. The tilt disturbance time (around 3.8th hour ±0.1 hour) was marked as the second abnormal period, and the drift rate during this period showed a short jump. The measurement points corresponding to these abnormal periods include P11, P12, P13 (temperature anomaly) and P9, P10 (tilt disturbance). The drift characteristics were extracted from the data analysis of the abnormal periods. The drift characteristics of the temperature abnormal period showed periodic fluctuations, and the fluctuation amplitude was proportional to the temperature change rate, with a correlation coefficient of 0.82. The drift characteristics of the tilt disturbance period showed pulse-like jumps, and the jump amplitude was related to the change in tilt angle. Based on these characteristics, the initial drift rate function was modified.

[0042] In the above embodiment, the weight adjustment can adopt a segmented weighting method. For the temperature abnormal period (4.5-5.2 hours), the weight coefficient of the period function is set to 0.6, meaning that the reliability of the data in this period is reduced by 40%. For the tilt disturbance period (3.7-3.9 hours), the weight coefficient is set to 0.7. The normal period maintains the weight coefficient 1.0. The adjusted drift rate function expression is: 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 retaining its trend information. The adjusted drift rate function more truly reflects the drift characteristics of the gravimeter under normal working conditions, achieving accurate modeling of the drift characteristics under complex measurement environments.

[0043] In an optional embodiment, the first definite integral operation is performed on the drift rate function in the time advancing direction from the measurement starting moment, to obtain the target cumulative drift amount of each measuring point, specifically including: 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 fast-changing period according to the gradient comparison result; performing a first numerical integral calculation on the stable function part in the stable period of the drift rate function by using a first integral step, to obtain a stable period integral value; performing a second numerical integral calculation on the fast-changing function part in the fast-changing period of the drift rate function 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; and accumulating the stable period integral value and the fast-changing period integral value in the time sequence of each measuring point on the time axis, to obtain the target cumulative drift amount.

[0044] In the above embodiment, the CG-6 relative gravimeter is used for mine area gravity measurement, and the established drift rate function v(t) = 0.0124 + 0.0019t + 0.000124t2- 0.0000072t3 is subjected to integral operation. First, the first derivative of the function is calculated to obtain v'(t) = 0.0019 + 0.000248t - 0.0000216t2, which represents the change rate of the drift rate. The preset gradient threshold is set to 0.00032 mGal / h2, which is determined based on the typical range of the drift rate change of the CG-6 gravimeter under normal working conditions. The gradient comparison result is 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. According to the gradient comparison result, the 8-hour measurement period is divided into: the stable period including [0, 1.2] hours and [6.8, 8] hours, in which the absolute value of the derivative is relatively small; and the fast-changing period of [1.2, 6.8] hours, which corresponds to the main working period of the instrument, and the drift rate changes significantly. The stable period accounts for 30% of the total time, and the fast-changing period accounts for 70%.

[0045] In the above embodiment, the first integral step is set to 0.05 hours (3 minutes), which is suitable for integral calculation in the stable period. The trapezoidal integral method is used to calculate the stable function part of [0, 1.2] hours: D1 = ∫0¹·²v(t)dt ≈ Σᵢ=0 23[v(t +1 )] / 2×0.05, the first stable period integral value D1=0.2mGal is calculated. Similarly, the second stable period integral value D2=0.3mGal is calculated for the [6.8, 8] hour interval. The second integral step is set to 0.01 hour (36 seconds) for fine integration of the fast varying interval. For the fast varying function part of [1.2, 6.8] hours, Simpson's integration 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 fast varying interval integral value D3=1.1mGal is calculated. The use of Simpson's integration method improves the integration accuracy in the interval where the drift rate changes rapidly. The target cumulative drift is obtained by time-sequential accumulation. For the measurement point P10 (measurement time is 2.8 hours), its cumulative drift is the integral value before this time: C(2.8)=D1+∫1.2 8 v(t)dt=0.2+0.3=0.5mGal. For the measurement point P15 (measurement time is 4.2 hours), the integral value of the stable interval [0, 1.2] and the partial integral value of the fast varying interval [1.2, 4.2] need to be accumulated: C(4.2)=D1+∫1.2 4 ·²v(t)dt=0.2+0.6=0.8mGal.

[0046] In the above embodiment, for measurement points spanning different time intervals, the calculation of the cumulative drift requires segmented processing. The cumulative drift of measurement point P20 (measurement time is 6.5 hours) is: C(6.5)=D1+∫1.2 6 · 5 v(t)dt=0.2+1.0=1.2mGal. The cumulative drift of measurement 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 process, the data acquisition system of CG-6 gravimeter recorded the original data at a frequency of 1 Hz, which provided sufficient data support for numerical integration. Through the adaptive step integration strategy, the calculation efficiency was optimized while ensuring the calculation accuracy. The final target cumulative drift of 25 measuring points ranges from 0.1 to 1.6 mGal, with an average value of about 0.9 mGal. These values accurately reflect the cumulative drift characteristics of the gravimeter during the entire measurement process, providing a reliable quantitative basis for subsequent drift correction.

[0047] In an optional embodiment, the time weight coefficients of the measuring points are determined according to the outbound measurement time proportion of the outbound measurement time in the total outbound measurement time and the return measurement time proportion of the return measurement time in the total return measurement time, and the time drift correction amount of each measuring point is determined according to the target cumulative drift and the time weight coefficient, which specifically includes: obtaining the outbound measurement end time, and performing a third difference calculation on 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 the return measurement end time, and performing a fourth difference calculation on the return measurement end time and the return measurement start time to obtain the total return measurement time; dividing the outbound measurement time interval by the total outbound measurement time to obtain the outbound measurement time proportion; performing a fifth difference calculation on the return measurement end time and the return measurement time to obtain the return measurement time interval; dividing the return measurement time interval by the total return measurement time to obtain the return measurement time proportion; taking the target cumulative drift as the theoretical drift correction amount; determining the time weight coefficient according to the outbound measurement time proportion and the return measurement time proportion; determining the outbound theoretical correction amount and the return theoretical correction amount according to the theoretical drift correction amount and the time weight coefficient; comparing the forward and return gravity values of the plurality of measuring points to obtain the forward and return gravity closure error; and correcting the outbound theoretical correction amount and the return theoretical correction amount according to the forward and return gravity closure error to obtain the time drift correction amount.

[0048] In the above embodiment, the CG-6 relative gravimeter is used for gravity measurement in a certain iron mine area. The starting time of the outbound measurement is 7:30:00, and the ending time of the outbound measurement is 11:45:30. The total time of the outbound measurement is calculated by the difference between the two times, which is 4.258 hours. The starting time of the return measurement is 12:15:00 (including 30 minutes of rest adjustment time), and the ending time of the return measurement is 16:28:45. The total time of the return measurement is calculated by the difference between the two times, which is 4.229 hours. Taking the measurement point P12 as an example, the outbound measurement time is 9:42:18, the time interval from the starting time of the outbound measurement is 2.206 hours, and the outbound measurement time ratio is obtained by dividing the time interval by the total time of the outbound measurement, which is 4.258 hours, and the result is 0.518. The return measurement time of the measurement point is 14:16:27, and the time interval between the ending time of the return measurement 16:28:45 and the starting time of the return measurement is 2.204 hours. The return measurement time ratio is obtained by dividing the time interval by the total time of the return measurement, which is 4.229 hours, and the result is 0.521. The calculation of the time ratio reflects the relative position of the measurement point in the entire measurement process. The target cumulative drift amount is directly obtained by the foregoing calculation. The time interval from the starting time of the outbound measurement to the outbound measurement time of the measurement point P12 is 2.206 hours. The integral of the drift rate function v(t) = 0.0085 + 0.0013t + 0.000084t²- 0.0000049t³ in the interval [0, 2.206] is calculated, and the result is that the target cumulative drift amount of the outbound measurement is 0.62 mGal. The time interval from the starting time of the return measurement to the return measurement time of the measurement point is 2.204 hours. The integral of the drift rate function in the interval [0, 2.204] is calculated, and the result is that the target cumulative drift amount of the return measurement is 0.61 mGal. The two cumulative drift amounts reflect the cumulative drift effect of the quartz spring system of the gravimeter under different working times.

[0049] In the above embodiment, the time weight coefficient is determined according to the proximity of the outbound and return measurement time ratios. For the measurement points with a time ratio difference less than 0.05, the average weight method is used: the outbound weight w f =0.5, and the return weight w b =0.5. For the measurement points with a time ratio difference 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 is the outbound measurement time ratio, and P b is the return measurement time ratio. The ratio difference of the measurement point P12 is |0.518-0.521|=0.003, which is less than 0.05, so w f =w b= 0.5. The forward and return theoretical corrections are determined according to the target cumulative drift. In this embodiment, the theoretical corrections are directly equal to the target cumulative drift, without involving the adjustment of the time weighting coefficient. The forward theoretical correction Cf,12of the station P12is 0.62 mGal, and the return theoretical correction Cb,12is 0.61 mGal. These two theoretical corrections embody the theoretical prediction of the drift influence on the station P12at the forward and return measurement time based on the drift rate function model. The forward and return consistency deviation is obtained by comparing the relative gravity values of the same station at the forward and return times. The relative gravity value of the station P12at the forward time is 1842.4 mGal, and the relative gravity value at the return time is 1842.9 mGal, with a difference of 0.5 mGal. After applying the theoretical corrections, the temporary corrected gravity value at the forward time is 1842.4 - 0.62 = 1841.78 mGal, and the temporary corrected gravity value at the return time is 1842.9 - 0.61 = 1842.29 mGal, with a forward and return consistency deviation of |1841.78 - 1842.29| = 0.51 mGal. The forward and return consistency deviations of the 25 stations along the whole profile range from 0.28 mGal to 0.85 mGal, with a root mean square error RMSE = 0.62 mGal. The existence of the forward and return consistency deviation is mainly caused by three factors: the nonlinear characteristics of the instrument drift, the influence of the ambient temperature change, and the interference of the ground micro-vibration. The forward and return consistency deviation reflects the deviation between the theoretical drift model and the actual measurement, and the forward and return theoretical corrections need to be corrected.

[0050] In the above embodiments, the specific process of modifying the outbound theoretical correction and the return theoretical correction according to the round-trip consistency deviation is based on the least square principle. The goal of the modification is to make the modified outbound and return corrected gravity values as close to the true gravity values as possible, while minimizing the deviation of the correction from the theoretical model. The modified model is established as follows: let the outbound correction of the measuring point i be ΔCf,i, and the return correction be ΔCb,i, then the modified outbound corrected gravity value is Gf,i-Cf,i-ΔCf,i, and the return corrected gravity value is Gb,i-Cb,i-ΔCb,i, wherein Gf,i and Gb,i are the outbound and return relative gravity values respectively, and Cf,i and Cb,i are the outbound and return theoretical correction values respectively. The constraint conditions of the modification include: 1) round-trip consistency constraint: the modified outbound and return corrected gravity values should be equal, that is, (Gf,i-Cf,i-ΔCf,i)=(Gb,i-Cb,i-ΔCb,i), which is transformed to Δ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 is simplified as ΔCf,i-ΔCb,i=Δi; 2) correction minimum constraint: the correction should be as small as possible to maintain the effectiveness of the theoretical model, and the objective function is minΣ(ΔCf,i²+ΔCb,i²); 3) time weight constraint: the distribution of the correction should consider the measurement timing and the position of the measuring point, the outbound measurement time is earlier, the prediction accuracy of the theoretical model is higher, and the correction should be smaller; the return measurement time is later, the uncertainty of the cumulative drift increases, and the correction should be larger.

[0051] In the above embodiments, based on the above constraint conditions, the weighted least square method is used to solve the optimal correction. The time weight coefficients w f and w b are introduced as the weight factors of the correction distribution, and the weighted Lagrange function is constructed: L=Σ(w f ·ΔCf,i²+w b ·ΔCb,i²)+λ·Σ[ΔCf,i-ΔCb,i-Δi]. The partial derivatives of ΔCf,i and ΔCb,i are taken respectively and set to zero, and the following is obtained: ∂L / ∂ΔCf,i=2w f ·ΔCf,i+λ=0, ∂L / ∂ΔCb,i=2w b ·ΔCb,i-λ=0. The simultaneous equations are solved to obtain: ΔCf,i=-λ / (2w f ), ΔCb,i=λ / (2w b ). Substituting them into the round-trip consistency constraint ΔCf,i-ΔCb,i=Δi, the following is obtained: -λ / (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 for the outbound journey is smaller (ΔCf,i<Δi / 2), while the correction amount for the return journey is larger (|ΔCb,i|>Δi / 2). This is consistent with the physical intuition that "the outbound measurement is more reliable and the correction amount should be smaller".

[0052] 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 departure correction amount AFCf,12 = -0.51 x 0.5 / (0.5 + 0.5) = -0.51 / 2 = -0.255 mGal; the return correction amount AFCb,12 = -(-0.51) x 0.5 / (0.5 + 0.5) = 0.51 / 2 = 0.255 mGal. The time drift correction amount is obtained by correcting the theoretical correction amount with the round-trip consistency deviation. For the departure measurement: the departure time drift correction amount = the departure theoretical correction amount + the departure correction amount, the measurement point P12 departure: 0.62 + (-0.255) = 0.365 mGal. For the return measurement: the return time drift correction amount = the return theoretical correction amount + the return correction amount, the measurement point P12 return: 0.61 + 0.255 = 0.865 mGal. The corrected departure correction gravity value is 1842.4 - 0.365 = 1842.035 mGal, and the return correction gravity value is 1842.9 - 0.865 = 1842.035 mGal, which are completely consistent, and the round-trip consistency deviation is reduced to 0 mGal. This verifies the effectiveness of the weighted least squares correction method. According to the time weight coefficient, the corrected departure and return correction gravity values are weighted and averaged to obtain the final correction gravity value of the measurement point: the target correction gravity value = w f x 1842.035 + w b x 1842.035 = 0.5 x 1842.035 + 0.5 x 1842.035 = 1842.035 mGal. Since the corrected round-trip consistency deviation is zero, the weighted average result is the same as the single value.

[0053] In the above embodiment, the same correction method is applied to 25 measurement points of the entire measurement line. Before correction, the root mean square error RMSE0 of the round-trip consistency deviation is 0.62 mGal; after correction, the root mean square error RMSE1 of the round-trip consistency deviation is 0.08 mGal, with a reduction of 87%. The corrected round-trip consistency deviation is mainly derived from random errors and environmental disturbances in the measurement process, which is close to the reading resolution 0.1 mGal of the CG-6 type gravity meter. Through this correction method based on the weighted least squares principle, each measurement point obtains a drift correction value matched with its measurement time sequence and position. The correction process has a clear mathematical model and physical meaning: the weighted least squares principle ensures the optimality of the correction amount, and the time weight coefficient reflects the characteristics of the drift accumulation with time and the difference in measurement reliability, which, in combination, not only guarantees the theoretical continuity of the drift correction, but also eliminates the systematic deviation of the theoretical model through the feedback correction of the measured round-trip consistency deviation, ensuring the accuracy of the relative gravity measurement conversion.

[0054] In an optional embodiment, the en route gravity value and the return gravity value are drift-corrected by the time drift correction amount to obtain the target corrected gravity value, and the method further comprises: subtracting the en route time drift correction amount in the corresponding time drift correction amount from the en route gravity value to obtain an en route temporary corrected gravity value, and subtracting the return time drift correction amount in the corresponding time drift correction amount from the return gravity value to obtain a return temporary corrected gravity value; taking the en route temporary corrected gravity value and the return temporary corrected gravity value as iteration variables, and performing the following iteration optimization operation until the convergence eigenvalue of the round-trip consistency deviation meets the preset convergence condition: performing a sixth difference calculation on the en route temporary corrected gravity value and the return temporary corrected gravity value of the same measuring point to obtain the round-trip consistency deviation; determining the convergence eigenvalue of the round-trip consistency deviation according to the preset convergence index; when the convergence eigenvalue does not meet the preset convergence condition, adjusting the fitting parameters of the drift rate function according to the round-trip consistency deviation to obtain an adjusted drift rate function; performing a second definite integral operation on the adjusted drift rate function to obtain a new target cumulative drift amount; determining a new en route time drift correction amount and a new return time drift correction amount according to the new target cumulative drift amount and the time weight coefficient; drift-correcting the en route gravity value by the new en route time drift correction amount to obtain a new en route temporary corrected gravity value, and drift-correcting the return gravity value by the new return time drift correction amount to obtain a new return temporary corrected gravity value; and performing a weighted average on the en route temporary corrected gravity value and the return temporary corrected gravity value that meet the preset convergence condition according to the time weight coefficient to obtain the target corrected gravity value.

[0055] In the above embodiment, a CG-6 relative gravimeter is used for gravity measurement in a certain coal mine area. There are 25 measuring points in the measuring line, and the round-trip observation method is used. The en route relative gravity value of measuring point P18 is 2156.8 mGal, and the return relative gravity value is 2157.7 mGal. First, the initial time drift correction amount is calculated. The time interval of the en route measurement time of measuring point P18 from the en route starting time is 2.8 hours, and the initial drift rate function v(t) = 0.0085 + 0.0013t + 0.000084t2- 0.0000049t3 is integrated in the interval [0, 2.8] to obtain the en route target cumulative drift amount of 0.92 mGal. The time interval of the return measurement time of this measuring point from the return starting time is 2.6 hours, and the initial drift rate function is integrated in the interval [0, 2.6] to obtain the return target cumulative drift amount of 0.88 mGal. The time weight coefficient is determined according to the en route and return measurement time proportion. The en route measurement time proportion of measuring point P18 is 0.52, the return measurement time proportion is 0.54, and the proportion difference is 0.02, which is less than 0.05, so the average weight method is used: f b ​= 0.5. According to the technical solution of claim 6, the theoretical correction amount of the outgoing pass is equal to the target cumulative drift amount 0.92 mGal of the outgoing pass, and the theoretical correction amount of the return pass is equal to the target cumulative drift amount 0.88 mGal of the return pass. After applying the round-trip consistency deviation correction, the time drift correction amount of the outgoing pass of the measuring point P18 is 0.90 mGal, and the time drift correction amount of the return pass is 0.90 mGal. The temporary corrected gravity value of the outgoing pass is obtained by subtracting the time drift correction amount of the outgoing pass from the relative gravity value of the outgoing pass, that is, 2156.8-0.90 = 2155.90 mGal. The temporary corrected gravity value of the return pass is obtained by subtracting the time drift correction amount of the return pass from the relative gravity value of the return pass, that is, 2157.7-0.90 = 2156.80 mGal. The difference between the two temporary corrected values reflects the initial deviation of the drift model.

[0056] In the above embodiment, the round-trip consistency deviation is obtained by the sixth difference calculation. The round-trip consistency deviation of the measuring point P18 is |2155.90-2156.80| = 0.90 mGal. The range of the round-trip consistency deviation of the 25 measuring points of the whole measuring line is 0.2-1.2 mGal, and these deviation values reflect the imperfection of the initial drift correction. The preset convergence index adopts the root mean square error, and the calculation formula is: RMSE = V(∑(Δi2) / n), where Δi 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.78 mGal. The preset convergence condition is set to be that the root mean square error of the round-trip consistency deviation is <0.2 mGal, and this 0.2 mGal threshold value can be determined based on the reading resolution 0.1 mGal of the CG-6 type gravity meter and the actual measurement accuracy requirement. Since RMSE0 = 0.78 mGal is greater than the convergence condition 0.2 mGal, the parameter adjustment is needed. The original form of the drift rate function is v(t) = 0.0085 + 0.0013t + 0.000084t2-0.0000049t3. According to 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) the sign distribution of the deviation: among the 25 measuring points, the return corrected gravity value of 18 measuring points is greater than the outgoing corrected gravity value (positive deviation), and the return corrected gravity value of 7 measuring points is less than the outgoing corrected gravity value (negative deviation), and the positive deviation accounts for 72%, which indicates that the drift rate function overall underestimates the actual drift; 2) the spatial distribution of the deviation: the average deviation of the front section (measuring points 1-8) of the measuring line is 0.45 mGal, the average deviation of the middle section (measuring points 9-17) of the measuring line is 0.92 mGal, and the average deviation of the rear section (measuring points 18-25) of the measuring line is 0.68 mGal, and the deviation of the middle section is the largest, which indicates that the quadratic term coefficient needs to be increased; 3) the time trend of the deviation: the deviation increases first and then decreases with the increase of the measurement time, which indicates that the cubic term coefficient needs to be adjusted.

[0057] In the above embodiments, the parameter adjustment strategy is based on the sensitivity analysis of the round-trip consistency bias to each parameter. The sensitivity analysis is performed by calculating the partial derivatives of the bias to each parameter: ∂Δi / ∂a0≈1, ∂Δi / ∂a1≈ti, ∂Δi / ∂a2≈ti2, ∂Δi / ∂a3≈ti3. Based on the gradient descent method, the parameter adjustment amount is: Δa j =-α j ×Σ(Δᵢ×∂Δᵢ / ∂a j ) / Σ(∂Δᵢ / ∂a j )², where α j is the learning rate parameter. The specific calculation is 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 bias of all measurement points (considering the sign), and the learning rate α0=0.8. First-order coefficient adjustment: Δa1=-α1×Σ(Δᵢ×ti) / Σ(ti²)=-0.6×52.3 / 1458=-0.0000215. Where Σ(Δᵢ×ti)=52.3mGal·h is the weighted sum of the bias and time, Σ(ti²)=1458h² is the sum of squares of time, and the learning rate α1=0.6. Second-order coefficient adjustment: Δa2=-α2×Σ(Δᵢ×ti²) / Σ(ti 4 )=-0.4×(-128.6) / 68420=0.00000075. Where Σ(Δᵢ×ti²)=-128.6mGal·h², Σ(ti 4 )=68420h 4 , and the learning rate α2=0.4. The negative value indicates that the middle section bias is large, and the second-order coefficient needs to be increased. Third-order coefficient adjustment: Δa3=-α3×Σ(Δᵢ×ti³) / Σ(ti 6 )=-0.2×(-485.2) / 3125000=0.000000031. Where Σ(Δᵢ×ti³)=-485.2mGal·h³, Σ(ti 6 )=3125000h 6 , and the learning rate α3=0.2. The learning rate parameters α0, α1, α2, α3 are set to 0.8, 0.6, 0.4, 0.2, following the decreasing principle to ensure that the adjustment of high-order terms is more cautious and avoids overfitting. The adjusted drift rate function is: v （ ¹ ）(t) = (0.0085 - 0.00024) + (0.0013 - 0.0000215)t + (0.000084 + 0.00000075)t2+ (-0.0000049 + 0.000000031)t3= 0.00826 + 0.001278t + 0.000085t2- 0.0000049t3.

[0058] In the above embodiment, the second integral operation is performed on the adjusted drift rate function, Simpson's rule is used, and a new target cumulative drift is obtained. The time interval of the outbound measurement time of the measuring point P18 from the starting time is 2.8 hours. The integral of the adjusted drift rate function is performed in the interval [0, 2.8], and the new outbound cumulative drift is 0.87 mGal, which is reduced by 0.05 mGal compared with the original value of 0.92 mGal. The time interval of the return measurement time of the measuring point from the starting time of the return is 2.6 hours. The integral of the adjusted drift rate function is performed in the interval [0, 2.6], and the new return cumulative drift is 0.83 mGal, which is reduced by 0.05 mGal compared with the original value of 0.88 mGal. The new time drift correction is determined based on the new cumulative drift. In the iterative optimization process, the time drift correction is directly equal to the cumulative drift, and there is no adjustment of the time weight coefficient and no correction of the round-trip consistency deviation to maintain the purity of the model optimization. The new outbound time drift correction = the new outbound cumulative drift = 0.87 mGal; the new return time drift correction = the new return cumulative drift = 0.83 mGal. This direct use of the cumulative drift as the correction ensures that the iterative process is completely based on the optimization of the drift rate function model, avoiding the interference of human factors on the adjustment of the model parameters. After applying the new correction, the new outbound temporary corrected gravity value of the measuring point P18 is 2156.8-0.87 = 2155.93 mGal, and the new return temporary corrected gravity value is 2157.7-0.83 = 2156.87 mGal, and 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 significantly decreases. The RMSE1 after the first iteration is 0.51 mGal, which is reduced by 35% compared with the initial value of 0.78 mGal. This is because the parameter adjustment is mainly optimized for the measuring points with large deviations in the middle section. The deviation of the measuring points in the middle section is reduced from an average of 0.92 mGal to 0.58 mGal, with a decrease of 37%, while the deviations of the measuring points in the front and rear sections increase slightly, but the overall RMSE decreases significantly.

[0059] In the above example, the iteration process continues, and each iteration adjusts the drift rate function parameters based on the current round-trip consistency bias distribution, recalculates the cumulative drift and the correction. The second iteration uses the same gradient descent strategy, and the parameter adjustment is: Δ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.000086t2-0.0000048t3. After applying the correction of the second iteration, the outbound temporary corrected gravity value at point P18 is 2156.05 mGal, the inbound temporary corrected gravity value is 2156.92 mGal, and the round-trip consistency bias is 0.87 mGal. The RMSE2=0.32 mGal after the second iteration, which is 37% lower than the first iteration. The parameter adjustment of the third iteration is: Δ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.000086t2-0.0000048t3. After applying the correction of the third iteration, the outbound temporary corrected gravity value at point P18 is 2156.12 mGal, the inbound temporary corrected gravity value is 2156.95 mGal, and the round-trip consistency bias is 0.83 mGal. The RMSE3=0.19 mGal after the third iteration, which is 41% lower than the second iteration. The parameter adjustment of the fourth iteration is: Δ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.000087t2- 0.0000047t3. After applying the correction of the fourth iteration, the temporary corrected gravity values of P18 are 2156.16 mGal for the outbound and 2156.97 mGal for the inbound, with a round-trip consistency deviation of 0.81 mGal. The RMSE4after the fourth iteration is 0.15 mGal, which is 21% lower than 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.000087t2- 0.0000047t3. The convergence of the iteration process is verified by the monotonic decrease of RMSE: RMSE0= 0.78 mGal→ RMSE1= 0.51 mGal→ RMSE2= 0.32 mGal→ RMSE3= 0.19 mGal→ RMSE4= 0.15 mGal, with a total decrease of 81%.

[0060] In the above embodiment, after the convergence condition is met, the final cumulative drift of each measuring point is calculated based on the final optimized drift rate function. The time interval between the departure measurement time of measuring point P18 and the starting time is 2.8 hours, and the final drift rate function is integrated in the interval [0, 2.8] to obtain the final departure cumulative drift of 0.82 mGal. The time interval between the return measurement time of the measuring point and the starting time of the return is 2.6 hours, and the final drift rate function is integrated in the interval [0, 2.6] to obtain the final return cumulative drift of 0.78 mGal. The original gravity value is corrected using the final cumulative drift: the temporarily corrected gravity value after convergence of the departure is 2156.8-0.82=2155.98 mGal; the temporarily corrected gravity value after convergence of the return is 2157.7-0.78=2156.92 mGal. At this time, the one-way consistency deviation is |2155.98-2156.92|=0.94 mGal. Although the one-way consistency deviation of a single measuring point P18 increases slightly from the initial value of 0.90 mGal to 0.94 mGal, the root mean square error of the entire measuring line decreases significantly from 0.78 mGal to 0.15 mGal, with a decrease of 81%, indicating that the iterative optimization method has obvious effect on drift correction of the entire measuring line. The target corrected gravity value is obtained by weighted averaging the departure temporarily corrected gravity value and the return temporarily corrected gravity value according to the time weight coefficient. The time weight coefficient of measuring point P18 is wf=wb=0.5, so the target corrected gravity value is 0.5×2155.98+0.5×2156.92=2156.45 mGal. The role of the time weight coefficient is to combine the measurement results of the departure and the return after iterative convergence to obtain the final corrected gravity value. For measuring points with large time proportion difference, the time weight coefficient will tilt to the side with more reliable time sequence, thereby improving the accuracy of the final result. This target corrected gravity value is the relative gravity value relative to the starting point of the measuring line (i.e. the gravity secondary base point). If you want to get the absolute gravity value, you need to add the absolute gravity value of the starting point of the measuring line. Assuming that the absolute gravity value of the starting point of the measuring line measured by the FG5 type absolute gravity meter (accuracy better than 2 μGal) is 979254.3 mGal, then the absolute gravity value of measuring point P18 is 979254.3+2156.45=981410.75 mGal. In this embodiment, the standard accuracy mode is taken as an example, and the iterative optimization of the drift rate function parameters is realized to achieve high consistency of the one-way measurement.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.

[0061] 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.

[0062] 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.

[0063] 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.

[0064] 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.

[0065] like Figure 2As shown, the electronic device includes a central processing unit (CPU) 201 which can perform various appropriate actions and processes, such as executing the methods described in the above embodiments, in accordance with a program stored in a read-only memory (ROM) 202 or a program loaded from a storage section 208 into a random access memory (RAM) 203. In the RAM 203, various data required for the system operation is also stored. The CPU 201, the ROM 202, and the RAM 203 are connected to each other through a bus 204. An input / output (I / O) interface 205 is also connected to the bus 204.

[0066] The following components are connected to the I / O interface 205: an input section 206 including an audio input device, a button switch, and the like; an output section 207 including a liquid crystal display (LCD), an audio output device, an indicator, and the like; a storage section 208 including a hard disk and the like; and a communication section 209 including a network interface card such as a LAN (Local Area Network) card, a modem, and the like. The communication section 209 performs communication processing via a network such as the Internet. A drive 210 is also connected to the I / O interface 205 as necessary. A removable media 211 such as a magnetic disk, an optical disk, a magneto-optical disk, a semiconductor memory, and the like is attached to the drive 210 as necessary, so that a computer program read therefrom is installed into the storage section 208 as necessary.

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

[0068] Note that specific examples of computer-readable storage media can include but are not limited to an electrical connection having one or more wires, a portable computer diskette, a hard disk, a random access memory (RAM), a read-only memory (ROM), an erasable programmable read-only memory (EPROM or Flash memory), an optical fiber, a portable compact disc read-only memory (CD-ROM), an optical storage device, a magnetic storage device, or any suitable combination of the foregoing. In the present disclosure, computer-readable storage media can be any tangible medium that can contain, or store a program for use by or in connection with an instruction execution system, apparatus, or device.

[0069] The flow diagrams and the block diagrams in the drawings are illustrations of architectures, functional processes, and operational processes, according to various embodiments of the present application. It will be understood that each block of the flow diagrams and / or block diagrams, and combinations of blocks in the flow diagrams and / or the block diagrams, can be implemented by computer readable program instructions such as program code. Such computer readable program instructions can be provided to a processor of a computer, or other programmable data processing apparatus, to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions / acts specified in the flow diagrams and / or block diagrams. These computer readable program instructions can also be stored in a computer readable storage medium that can direct a computer, a programmable data processing apparatus, and the other

[0070] Specifically, the electronic device of the embodiment includes a processor and a memory, and the memory stores a computer program. When the computer program is executed by the processor, the method for gravimeter drift correction geodetic survey based on a drift rate function is implemented.

[0071] As another aspect, the present application also provides a computer readable storage medium. The storage medium can be included in the electronic device described in the above embodiments, or can exist separately and not be assembled into the electronic device. The storage medium carries one or more computer programs. When the one or more computer programs are executed by a processor of the electronic device, the electronic device implements the method for gravimeter drift correction geodetic survey based on a drift rate function provided in the above embodiments.

[0072] The above embodiments are only used to illustrate the technical solutions of the present application, but not limit the present application; even though the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: they can still modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacements to some technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application.

[0073] Those skilled in the art can understand that all or part of the processes in the above-mentioned method embodiments can be implemented by a computer program instructing relevant hardware to complete, the program can be stored in a computer readable storage medium, and the program can include the processes of the above-mentioned method embodiments when executed. The aforementioned storage medium includes ROM or random storage memory RAM, magnetic disc or optical disc and various storage code medium.

Claims

1. A geodetic mapping method for gravimeter drift correction based on a drift rate function, characterized in that, include: Multiple measuring points are set up on a preset measuring line. Gravity meters are used 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. In the return order, the gravimeter is used to measure the return gravity of each measuring point in the reverse order to obtain the return gravity value and the corresponding return measurement time of each measuring point. The target instantaneous drift rate of each measuring point is determined based on the return gravity value, the outgoing gravity value, the outgoing measurement time, and the return measurement time. A drift rate function is established 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. Starting from the measurement start time, the drift rate function is integrated along the time progression direction to obtain the target cumulative drift amount of each measurement point. The target cumulative drift amount represents the integral value of the drift amount of the measurement point from the measurement start time to the current measurement time along the measurement time axis. 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, the time weight coefficient of each measurement point is determined, and the time drift correction amount of each measurement 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.

2. The method according to claim 1, characterized in that, The determination of 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 specifically includes: Subtracting the return gravity value from the outbound gravity value yields the gravity difference at each measuring point. The first difference between the return measurement time and the outbound measurement time is calculated to obtain the round-trip time interval of each measurement point. Divide the gravity difference by the round-trip time interval to obtain the average drift rate of each measuring point; The movement time of the gravimeter between every two adjacent measuring points in the plurality of measuring points is obtained, and the dwell time of the gravimeter at each measuring point is obtained; The actual observation period for each measuring point is determined based on the movement time and the dwell observation time. The average drift rate of each measuring point is correlated with the midpoint time of the corresponding actual observation period to construct a discrete drift rate dataset; Interpolation calculations are performed on the discrete drift rate dataset to obtain a continuous drift rate time function; Extract the function value of the midpoint of the actual observation period of each measuring point from the drift rate time function, and use the function value as the initial instantaneous drift rate of each measuring point; Acquire historical drift characteristic data of the gravimeter, and determine the drift correction coefficient based on the historical drift characteristic data; The initial instantaneous drift rate is corrected using the drift correction coefficient to obtain the target instantaneous drift rate.

3. The method according to claim 2, characterized in that, The step of determining the drift correction coefficient based on the historical drift characteristic data specifically includes: Extract the drift rate variation curves of the gravimeter under different operating durations from the historical drift characteristic data; The drift acceleration characteristic value of the gravimeter is determined based on the drift rate change curve. Obtain the cumulative working time of the current measurement task, and determine the drift rate value corresponding to the cumulative working time from the drift rate change curve; The drift rate value is compared with a preset standard drift rate to obtain a drift rate ratio. Based on the drift acceleration characteristic value and the drift rate ratio, the drift correction coefficient K is determined using 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 cumulative working time, and T0 is the preset standard working time.

4. The method according to claim 1, characterized in that, The establishment of a drift rate function, 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 includes: A second difference is calculated between the outbound measurement time and the outbound measurement start time to obtain the outbound measurement time interval for each measurement point. Determine the instantaneous drift rate difference between any two adjacent measuring points in the plurality of measuring points, and determine the outbound measurement time interval difference between any two adjacent measuring points in the plurality of measuring points; The drift acceleration is determined based on the difference between the instantaneous drift rate and the difference between the outward measurement time interval. Using the outbound measurement time interval as the independent variable and the target instantaneous drift rate as the dependent variable, an initial drift rate function is established using a preset numerical fitting method; The working status parameters of the gravimeter are obtained, and the period of abnormal drift of the gravimeter is determined based on the working status parameters. The drift characteristics of the gravimeter are extracted from the drift anomaly period, and the abnormal function portion of the initial drift rate function located in the drift anomaly period is weighted according to the drift characteristics to obtain the drift rate function.

5. The method according to claim 1, characterized in that, The step of performing a first definite integral operation on the drift rate function along the time progression direction starting from the measurement start time to obtain the target cumulative drift amount at each measurement point specifically includes: The first derivative of the drift rate function is determined, and the absolute value of the first derivative is compared with a preset gradient threshold to obtain the gradient comparison result. Based on the gradient comparison results, the drift rate function is divided into a stable period and a rapidly changing period; The first numerical integration is performed on the stationary part of the drift rate function located in the stationary period using the first integration step size to obtain the stationary period integral value. A second numerical integration is performed on the fast-changing function portion of the drift rate function located within the fast-changing time period using a second integration step size to obtain the fast-changing time period integral value. The second integration step size is smaller than the first integration step size. The integral values ​​of the stable period and the integral values ​​of the rapidly changing period are accumulated according to the time sequence of each measuring point on the time axis to obtain the cumulative drift of the target.

6. The method according to claim 1, characterized in that, The step of determining the time weighting coefficient for each measurement 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 determining the time drift correction amount for each measurement point based on the target cumulative drift amount and the time weighting coefficient, specifically includes: Obtain the end time of the outbound measurement, and calculate the third difference between the end time of the outbound measurement and the start time of the outbound measurement to obtain the total duration of the outbound measurement; Obtain the start time and end time of the return trip measurement, and perform a fourth difference calculation between the end time and the start time to obtain the total duration of the return trip measurement; Divide the outbound measurement time interval by the total outbound measurement duration to obtain the outbound measurement time ratio; The fifth difference is calculated between the end time of the return measurement and the time of the return measurement to obtain the return measurement time interval; Divide the return trip measurement time interval by the total return trip measurement duration to obtain the return trip measurement time ratio; The target cumulative drift amount is used as the theoretical drift correction amount; The time weighting coefficient is determined based on the outbound measurement time ratio and the return measurement time ratio. The theoretical drift correction amount and the theoretical return correction amount are determined based on the theoretical drift correction amount and the time weighting coefficient. By comparing the forward and backward gravity values ​​at the multiple measuring points, the forward and backward gravity closure difference is obtained. The time drift correction is obtained by correcting the outbound theoretical correction and the return theoretical correction based on the round-trip gravity closure difference.

7. The method according to claim 1, characterized in that, The method of using the time drift correction amount to correct the drift of the outbound gravity value and the return gravity value to obtain the target corrected gravity value further includes: Subtract the outbound time drift correction amount from the outbound gravity value to obtain the outbound temporary correction gravity value; and subtract the return time drift correction amount from the return gravity value to obtain the return temporary correction gravity value. Using the outbound temporary correction gravity value and the return temporary correction gravity value as iterative variables, perform the following iterative optimization operation until the convergence characteristic value of the round-trip consistency deviation satisfies the preset convergence condition: The sixth difference is calculated by comparing the outbound temporary correction gravity value and the return temporary correction gravity value at 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 feature 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. Perform a second definite integral operation on the adjusted drift rate function to obtain a new target cumulative drift amount; The new outbound time drift correction and the new return time drift correction are determined based on the new target cumulative drift amount and the time weighting coefficient. The outbound gravity value is corrected by the new outbound time drift correction amount to obtain a new outbound temporary correction gravity value. The return gravity value is corrected by the new return time drift correction amount to obtain a new return temporary correction gravity value. The target corrected gravity value is obtained by weighting the outbound temporary correction gravity value and the return temporary correction gravity value that meet the preset convergence condition according to the time weighting coefficient.

8. An electronic device, characterized in that, The electronic device includes: one or more processors and a memory; the memory is coupled to the one or more processors, the memory is used to store computer program code, the computer program code including computer instructions, and the one or more processors call the computer instructions to cause the electronic device to perform the method as described in any one of claims 1-7.

9. A computer-readable storage medium comprising instructions, characterized in that, When the instructions are executed on an electronic device, the electronic device causes the electronic device to perform the method as described in any one of claims 1-7.

10. A computer program product, characterized in that, When the computer program product is run on an electronic device, it causes the electronic device to perform the method as described in any one of claims 1-7.

Images