Gravity anomaly correction method and system for high-altitude tunnel distance measurement

Abstract

The invention discloses a gravity anomaly correction method and system for high-altitude tunnel distance measurement, and relates to the technical field of engineering surveying and physical geodetic surveying. The method comprises the following steps: firstly, judging the type of a measurement section, and constructing potential constraints by using a global background gravity field model and a digital elevation model; performing gravitational field forward modeling based on physical model parameters, generating a shaft axis gravity or tunnel axis vertical line deviation sequence, and performing physical reduction on ranging data; actual measurement gravity is collected, a residual error sequence is constructed, the residual error sequence is mapped into a ranging correction error through a sensitivity matrix, and a through error projection value is obtained through propagation; and finally, taking the minimum residual error as a target to invert and update the parameter until the through error meets a threshold value. According to the method, the closed-loop iteration of measurement section-strategy-index is established, so that the data calibers of the vertical shaft and the adit are effectively unified, and the influence of gravity anomaly on the penetration precision is inhibited.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of engineering surveying and physical geodesy, and in particular to a gravity anomaly correction method and system for high-altitude tunnel distance measurement. BACKGROUND

[0002] In high-altitude and gravity anomaly obvious mountain tunnel engineering, vertical shaft and horizontal tunnel often exist at the same time: vertical shaft is used for vertical control transmission, and horizontal tunnel is used for horizontal control extension. Due to the combined action of altitude gradient, terrain undulation and local density anomaly, the spatial variation of gravity field is significant, which leads to systematic deviation of distance measurement observation in the process of converting from geometric observation to unified reference quantity required by engineering control, and the deviation is propagated and accumulated in the traverse network, and finally the through error is out of limit.

[0003] The common practice in existing engineering includes using static Bouguer correction or using empirical correction with fixed parameters to compensate for distance measurement, but these methods often have difficulty in simultaneously considering the sensitivity of vertical shaft measurement section to vertical gravity gradient, the sensitivity of horizontal tunnel measurement section to vertical deflection (deviation of plumb line from reference normal), and the comprehensive constraint relationship of through error as the final engineering index to various error terms, thereby leading to lack of pertinence and iterative convergence of correction strategy. SUMMARY

[0004] The present application aims to provide a technical solution to solve the above problems in the prior art. Specifically, the present application is realized by the following technical solution: A gravity anomaly correction method for high-altitude tunnel distance measurement, comprising: S1, obtaining the original distance measurement observation data and the spatial position information of the measurement station of the to-be-processed measurement section, determining whether the measurement section is a vertical shaft transmission measurement section or a horizontal tunnel extension measurement section, calculating the potential of each measurement station based on a global background gravity field model and a digital elevation model of the measurement area, and constructing the to-be-corrected input data with potential difference constraint or potential elevation constraint; S2, obtaining initial physical model parameters, performing gravity field forward calculation on the measurement area: generating vertical shaft axis gravity distribution sequence for vertical shaft transmission measurement section, generating vertical line deflection distribution sequence along tunnel axis for horizontal tunnel extension measurement section, and performing vertical depth reduction or length reduction on the to-be-corrected input data according to the sequence to obtain physical reduction distance measurement data; S3, collecting the measured gravity observation value at the measurement station, and calculating the gravity residual sequence by subtracting the same coordinate base theoretical gravity observation value calculated from the physical model parameters from the measured gravity observation value, calculating the distance correction residual error from the gravity residual sequence, and obtaining the through error projection value through error propagation of the traverse network; S4, judging whether the through error projection value is less than a preset threshold value; if yes, outputting physical reduction ranging data; if no and the iteration termination condition is not met, updating the physical model parameters by inversion with the objective of minimizing gravity residual, and returning to S2 until the iteration termination condition is met, and outputting physical reduction ranging data containing an alarm identifier.

[0005] Further, the constructing of the input data to be corrected comprises: Based on the global background gravity field model, an earth gravity potential reference frame is constructed, and potential values of each station are calculated; when the survey section is a vertical shaft transfer survey section, a potential difference constraint is formed based on the potential values of the upper and lower stations, and the potential difference constraint is associated with the vertical shaft geometric depth observation to generate vertical input data to be corrected; when the survey section is a flat tunnel extension survey section, physical elevations of the stations relative to a preset reduction reference surface are calculated based on the potential values of the stations, and a reduction scale constraint is established by using the physical elevations on the flat tunnel geometric length observation to generate flat tunnel input data to be corrected.

[0006] Further, the obtaining of the initial physical model parameters comprises gravity field forward calculation of the survey area: a vertical shaft axis gravity distribution sequence is generated for the vertical shaft transfer survey section, and a vertical deflection distribution sequence along the tunnel axis is generated for the flat tunnel extension survey section, comprising: The physical model parameters comprise local lithology density parameters, regional background density parameters, and vertical gradient correction parameters; when the survey section is a vertical shaft transfer survey section, gravity acceleration of each depth point of the vertical shaft axis is calculated based on the vertical gradient correction parameters and in combination with the local lithology density parameters to form the vertical shaft axis gravity distribution sequence; when the survey section is a flat tunnel extension survey section, near-zone topographic mass gravity contribution is calculated based on the local lithology density parameters, far-zone Bouguer plate mass gravity contribution is calculated based on the regional background density parameters, and the gravity contribution is subjected to integral modeling to obtain a vertical deflection component at each survey point along the tunnel axis to form the vertical deflection distribution sequence.

[0007] Further, the performing of vertical depth reduction or length reduction on the input data to be corrected according to the sequence to obtain physical reduction ranging data comprises: When the survey section is a vertical shaft transfer survey section, gravity path integration is performed on the vertical shaft axis based on the vertical shaft axis gravity distribution sequence to obtain a potential difference theoretical value between stations, and the vertical shaft geometric depth observation is reduced under the potential difference constraint to obtain vertical physical reduction ranging data; When the survey section is a flat tunnel extension survey section, an inclination correction term and a reduction correction term are calculated based on the vertical deflection distribution sequence, the flat tunnel geometric length observation is corrected and reduced to a preset reference surface to obtain flat tunnel physical reduction ranging data.

[0008] Further, the measured gravity observation value includes at least one of a gravity acceleration scalar value and a gravity vector component observation value; the theoretical gravity observation value is a corresponding scalar value or a corresponding component value in the same coordinate basis as the measured gravity observation value; and the gravity residual sequence is constructed by a difference between the measured gravity observation value and the theoretical gravity observation value.

[0009] Further, the ranging correction residual error is calculated from the gravity residual sequence, and a through error projection value is obtained by means of error propagation of the traverse network, including: The gravity residual sequence is decomposed to obtain a high-frequency fluctuation component and a low-frequency trend component; a mapping sensitivity matrix of residual samples to ranging observation edges is constructed to map the high-frequency fluctuation component as a random disturbance item and the low-frequency trend component as a systematic disturbance item, and a ranging correction residual error affecting the physical reduction ranging data is calculated; the ranging correction residual error is taken as an input vector to substitute into a traverse network adjustment or error propagation function to calculate a transverse position deviation projection value and a longitudinal position deviation projection value at a through surface, and to determine the through error projection value according to the transverse position deviation projection value and the longitudinal position deviation projection value.

[0010] Further, when the measuring section is a flat tunnel extension measuring section, the calculation domain of the integral modeling is mutually exclusive divided into a near-zone integral domain and a far-zone integral domain, and an adaptive boundary radius is introduced; in the closed loop iteration process, the adaptive boundary radius is updated according to the spectral characteristics of the gravity residual sequence to suppress repeated counting or missing of the near-zone topographic mass gravity contribution and the far-zone Bouguer plate mass gravity contribution.

[0011] Further, the method further includes: Robustly weighting the gravity residual sequence and determining an abnormal section: the gravity residual is standardized based on the observation standard deviation of the gravity observation value, and a segmented weight function is used to obtain a robust weight; when the robust weight of consecutive stations is lower than an abnormal threshold, the corresponding interval is marked as an abnormal section; when updating the physical model parameters in the inversion, the objective function is weighted based on the robust weight, and segmented parameterized updating is performed on the local lithology density parameters corresponding to the abnormal section.

[0012] Further, the iteration termination condition includes: The number of iterations reaches a preset maximum value, the physical model parameter update amount of adjacent two iterations is less than a convergence threshold, or the through error projection value of adjacent two iterations decreases by less than a preset proportion threshold.

[0013] Further, the method further includes: Determining the measuring section as a vertical shaft transfer measuring section or a flat tunnel extension measuring section by using an average slope threshold of edges, including: When the average slope of the edges is greater than a preset slope threshold, the measuring section is determined as the vertical shaft transfer measuring section, otherwise, the measuring section is determined as the flat tunnel extension measuring section. in For the station Slope distance observation between , for The elevation difference between the two stations.

[0014] Furthermore, when the measurement section is a vertical shaft transfer measurement section, the theoretical value of the inter-station potential difference is obtained by performing gravity path integration on the vertical shaft axis based on the gravity distribution sequence of the vertical shaft axis, and the vertical shaft geometric depth observation is reduced under the potential difference constraint to obtain vertical physical reduced distance data, including: The theoretical value of the inter-station potential difference is obtained by performing gravity path integration on the shaft axis based on the gravity distribution sequence. : in The number of discrete depth points along the shaft axis. , For the first Each depth point; in the vertical shaft survey section, according to the upper station → Next stop Order definition and order and Take the positive value accumulated in the same direction; if the order of the endpoints does not meet the above convention, use... and Involved in ratio calculation; Under potential difference constraints, geometric depth observations The vertical physical distance measurement data was obtained by performing the reduction calculation. : in For adjacent stations Potential difference constraint.

[0015] Furthermore, when the measured section is an extension of a horizontal tunnel, the tilt correction term and the reduction correction term are calculated based on the vertical deviation distribution sequence. The observed geometric length of the horizontal tunnel is corrected and reduced to a preset reference surface to obtain the physical reduction distance data of the horizontal tunnel, including: The tilt correction term and reduction correction term are calculated using the vertical deviation distribution sequence. The observed geometric length of the horizontal tunnel is corrected and reduced to a preset reference surface. The projection angle of the vertical deviation on the side azimuth is calculated. in , ; Define the equivalent zenith angle under the normal basis: Get the geometric horizontal length under the normal base: According to the average value of the physical elevation of the two end stations Calculate the length to the preset reference surface: wherein is the equivalent curvature radius corresponding to the preset reference surface or the constant of the average radius of the earth.

[0016] A gravity anomaly correction system for high-altitude tunnel ranging, applying the gravity anomaly correction method for high-altitude tunnel ranging, comprising: a data initialization module, a strategy generation and reduction module, a residual and accuracy evaluation module, a closed-loop control and output module, and a data processing module; The data initialization module, the strategy generation and reduction module, the residual and accuracy evaluation module, and the closed-loop control and output module are respectively connected with the data processing module; The data initialization module is used to obtain original ranging observation data and station information, determine the type of the section, and construct the input data to be corrected under the constraint of potential difference or potential elevation; The strategy generation and reduction module is used to perform gravity field forward calculation based on physical model parameters to generate a vertical shaft axis gravity distribution sequence or a tunnel axis vertical deviation distribution sequence, and to perform physical reduction on geometric observation data to obtain physical reduction ranging data; The residual and accuracy evaluation module is used to construct a gravity residual sequence and calculate ranging correction residual error and through error projection value; The closed-loop control and output module is used to determine whether the through error projection value is less than a preset threshold; if yes, the physical reduction ranging data is output; if no and the iteration termination condition is not met, the physical model parameters are updated and a new round of calculation is triggered; if no and the iteration termination condition is met, the physical reduction ranging data containing an alarm identifier is output.

[0017] Preferably, if no and the iteration termination condition is not met, the physical model parameters are updated, comprising: Selecting an update object according to a dominant component: When the ranging correction residual error dominant through error projection value corresponding to the high-frequency fluctuation component is over-limit, updating the local lithology density parameter corresponding to the spatial position ; When the ranging correction residual error dominant over-limit corresponding to the low-frequency trend component of the tunnel extension section is over-limit, updating the regional background density parameter ; When the residual error of the low-frequency trend component of the shaft transfer section dominates the overrun, update the vertical gradient correction parameter .

[0018] When the termination condition is met but the residual error is still , output the physical reduction range data with residual error alarm mark, The through error projection threshold.

[0019] Preferably, the alarm mark at least includes the overrun type, the overrun amplitude and the corresponding section identifier.

[0020] Compared with the prior art, the present application has the following advantages and beneficial effects: the data caliber of the shaft and the tunnel is unified by the potential constraint; the closed-loop iteration is established by the section type-correction strategy-through index, the parameter update has a clear objective function and termination condition; the influence of gravity anomaly can be controllably transferred and suppressed from the observation level to the through index level, thereby improving the through control measurement accuracy and stability. BRIEF DESCRIPTION OF DRAWINGS

[0021] The drawings described herein are used to provide further understanding of the embodiments of the present application, form a part of the present application, and do not constitute a limitation on the embodiments of the present application. In the drawings: Figure 1 It is a flowchart of a gravity anomaly correction method for high-altitude tunnel ranging. DETAILED DESCRIPTION

[0022] In order to make the purpose, technical scheme and advantages of the present application clearer, further detailed description will be made below in combination with embodiments and drawings, the illustrative embodiments of the present application and the description thereof are only used to explain the present application, and do not constitute a limitation on the present application.

[0023] Embodiment 1 As Figure 1 shown, a gravity anomaly correction method for high-altitude tunnel ranging includes: S1, obtaining the original ranging observation data and the spatial position information of the measuring station of the section to be processed, determining whether the section is a shaft transfer section or a tunnel extension section, calculating the potential number of each measuring station based on the global background gravity field model and the digital elevation model of the survey area, and constructing the input data to be corrected based on the potential difference constraint or the potential elevation constraint; S2, obtaining the initial physical model parameters, performing gravity field forward calculation on the survey area: generating shaft axis gravity distribution sequence for shaft transfer section, generating vertical deflection distribution sequence along the tunnel axis for tunnel extension section, and performing vertical depth reduction or length reduction on the input data to be corrected according to the sequence to obtain physical reduction ranging data; S3. Collect the measured gravity observation values ​​at the station, and calculate the difference between them and the theoretical gravity observation values ​​calculated from the physical model parameters to obtain the gravity residual sequence. Calculate the distance measurement correction residual error from the gravity residual sequence, and obtain the penetration error projection value through the traverse network error propagation. S4, determine whether the projection value of the penetration error is less than the preset threshold; if yes, output the physical reduction distance data; if no and the iteration termination condition is not met, then update the physical model parameters with the goal of minimizing the gravity residual and return to S2, until the iteration termination condition is met and output the physical reduction distance data containing the alarm flag.

[0024] Specifically, in this embodiment, the data processing terminal can be an industrial control computer, a laptop, or a server, and it includes at least a processor and a memory, in which program instructions and data structures for executing this embodiment are stored.

[0025] The spatial location information of the station includes at least the three-dimensional coordinates of the station. ,in For the station The plane coordinates, where For the station Eastward coordinates For the station The northward coordinates (local tangent plane coordinate system). For the station The elevation can be used equivalently to calculate the height difference; the original distance measurement observation data should include at least the station pairs. Slope distance observation between and the azimuth angle used for projection calculation With zenith angle Zenith In this specification, the angle between the line of sight and the vertical downward direction is defined, and its value range is [value missing]. .

[0026] When the instrument outputs the zenith angle with vertical reference (downward viewing may occur). When converting to and with It participates in the calculation of vertical and horizontal components. The elevation difference between the two stations is denoted as... For each observation edge ,Record and endpoint station index .

[0027] The data processing terminal acquires the raw distance measurement data and station spatial location information of the section to be processed, and determines whether the section is a vertical shaft transfer section or a horizontal tunnel extension section. The determination rule can be based on tunnel axis segmentation marking, design mileage section identification, or the average slope threshold of the sides: when... When greater than the preset slope threshold, it is determined as a shaft transmission section, otherwise it is determined as a flat tunnel extension section.

[0028] The potential of each station is calculated based on a global background gravity field model and a digital elevation model of the survey area. The geodetic coordinates of the station (respectively, geodetic latitude, longitude and distance from the center of the earth) and the maximum order , based on spherical harmonic expansion, the , that is, the background potential provided by the global background gravity field model. Wherein Can take 360, 720 or 2160 and other engineering feasible values; the digital elevation model DEM can adopt regular grid form In this model, each grid cell (indexed by u, v) is constructed as a vertical prism cell for subsequent gravitational integral modeling; each grid cell is constructed as a vertical prism cell ; Or use the triangular net (TIN) form, extrude the triangular patches to form polyhedral cells. Cell geometric parameters include cell top and bottom surface elevation, plane range and volume When the calculation efficiency requirement is higher, the DEM in the far field can be down-sampled or layered. In order to avoid the ambiguity of the potential reference, the preset reference potential is defined as the station potential : Wherein, The gravitational potential of the station , unit ; The reference potential The reference potential constant corresponding to the engineering adopted reduction reference surface (preset equipotential surface) is given by the given vertical reference or design reference; In one implementation way, The background potential provided by the global background gravity field model And the topographic mass disturbance potential calculated by the digital elevation model Superimposed to get: Wherein, in order to avoid the repeated calculation of topographic effect between the background model and the terrain integral term, the background potential Adopt long wave / background caliber (for example, only keep the potential components of the background model within the preset bandwidth, or use the reference potential framework after removing the terrain), the topographic mass disturbance potential Residual / fine potential components calculated from digital elevation model in local range; in another implementation, remove-restore process can be used: first get reference potential from background model, then superimpose local residual potential calculated from DEM cell integration.

[0029] Terrain quality disturbing potential Can be calculated by cell integration: where, is the gravitational constant, is the density value of the th cell, is the volume of the cell, is the distance from the station to the cell infinitesimal. is the number of cells used for integration modeling, which can be prismatic or polyhedral elements generated from digital elevation model. In this embodiment, is the density value of the th cell, is the volume of the cell, is the distance from the station to the cell infinitesimal. is the number of cells used for integration modeling, which can be prismatic or polyhedral elements generated from digital elevation model. In this embodiment, is taken as the reference density used to build potential constraints , which is given by the initial geologic density model or lithology table and remains unchanged in the iteration process; the , in steps S2-S4 are taken as the physical model parameters to be updated, and the forward and residual interpretation for dynamic gravity correction strategy, where is the local lithology density parameter; is the regional background density parameter.

[0030] Build input data to be corrected: When passing the section for a shaft, build potential difference constraints for adjacent station pairs , and associate the constraints with the shaft geometric depth observations to generate vertical input data to be corrected; When the section is a flat tunnel extension section, calculate the physical elevation of the station relative to the reduced reference surface based on the potential number , and use the physical elevation to establish a reduced scale constraint for the flat tunnel geometric length observation to generate flat tunnel input data to be corrected. The physical elevation in this specification adopts the normal height or equivalent physical height caliber, which is calculated as: wherein, is the station the average normal gravity along the plumb line direction (from the reference ellipsoid to the station on the surface), which is the theoretical average of the gravity acceleration along the plumb line (from the reference ellipsoid to the station on the surface) path under the normal gravity field model, used to convert the gravity potential number into physical height.

[0031] The embodiment defines the physical model parameter vector: wherein, is the local lithology density parameter, is the regional background density parameter, is the vertical gradient correction parameter, with the dimension of , used to describe the residual vertical gradient correction term relative to the background model .

[0032] When the transmission section is a shaft, the depth points along the shaft axis are discretized , is the coordinate of the discretized depth point on the shaft axis, with the wellhead (the top of the shaft) as the zero point and the downward direction as the positive direction, obtained by discretizing the depth along the shaft axis; that is, the shaft axis is divided into M points, is the coordinate of the discretized depth point on the shaft axis, with the value of the depth value of the mth point; the gravity acceleration is calculated to form a gravity distribution sequence : wherein, is given by the global background gravity field model, is the local gravity correction term calculated by local mass volume integration (such as prism / polyhedral element gravity integration), affected by .

[0033] When the transmission section is a flat tunnel, the plumb deviation component is calculated for each station to form a plumb deviation distribution sequence . Wherein is taken as the north component, is taken as the east component; the survey line azimuth is taken as the azimuth (consistent with the coordinate convention of ) from north clockwise to the survey line direction In one implementation, the near-zone terrain mass gravity contribution and the far-zone Bouguer plate mass gravity contribution Then, the vertical deviation component is calculated from the horizontal derivative of the perturbation potential. During this process, a local tangent plane coordinate system is established, and the following definitions are made: The axis points eastward. If the axis points north, then: in, To disturb the potential. This is the normal gravity at the station. It can be obtained by converting the horizontal gravitational component from the mass volume integral model of the near and far regions.

[0034] The above The relationship with the derivative of the perturbation potential holds approximately within the local engineering area based on the local tangent plane (ENU); when the survey area spans a large area, a latitude-related spherical coordinate conversion factor (e.g., ...) can be further introduced. (etc.) to maintain strict physical consistency.

[0035] Vertical shaft transfer measurement section: The theoretical value of the potential difference between stations is obtained by performing gravity path integration on the vertical shaft axis based on the gravity distribution sequence. : in The number of discrete depth points along the shaft axis. , For the first Each depth point. In the vertical shaft survey section, according to the upper station... → Next stop Order definition and order and Take the positive value accumulated in the same direction; if the order of the endpoints does not meet the above convention, use... and It is used in the ratio calculation.

[0036] Under potential difference constraints, geometric depth observations The vertical physical distance measurement data was obtained by performing the reduction calculation. One feasible reduction method is scale-constrained reduction: when If the value is too small or abnormal, an alarm is triggered and the edge is assigned a low weight.

[0037] Extended measurement section of the horizontal tunnel: The tilt correction and reduction correction terms are calculated using the vertical deviation distribution sequence to correct the observed geometric length of the horizontal tunnel and reduce it to a preset reference surface. First, the projection angle of the vertical deviation on the side azimuth is calculated: wherein , .

[0038] Definition of equivalent zenith angle under normal base: Accordingly, the geometric horizontal length under the normal base is obtained: According to the average value of the physical elevation of the two end stations The length is reduced to the preset reference surface: wherein is the equivalent curvature radius corresponding to the preset reference surface or the constant of the average radius of the earth. Finally, the physical reduction surveying data of the flat tunnel is obtained The physical reduction surveying data of the vertical shaft and the flat tunnel are collectively referred to as the physical reduction surveying data.

[0039] Collecting the measured gravity observation value at the station When it is a scalar observation, is the gravity acceleration scalar; when it is a vector component observation, contains the component observation values under the same coordinate base. The theoretical gravity observation value is calculated from the physical model parameters The gravity residual is constructed: When the measured gravity observation value is a vector component observation, the residual samples are unfolded into a one-dimensional residual sequence according to (station index and component index ), that is, let , and the standard deviation and the robust weight are defined accordingly; or the vector residual is taken as a norm to form a scalar residual , and one of the two ways is selected. The sorted along the section is constructed into a gravity residual sequence. The gravity residual sequence is decomposed to obtain the high-frequency fluctuation component and the low-frequency trend component: In one implementation, for the flat tunnel extension section, the is subjected to multi-scale spectrum analysis or wavelet decomposition, the short spatial scale component is extracted as , and the long spatial scale component is extracted as ; for the vertical shaft transmission section, the is subjected to trend decomposition and mutation detection, the mutation component is extracted as , and the linear or polynomial trend component is extracted as wherein is the high-frequency fluctuation component, is the low-frequency trend component The residual is mapped to the misclosure residual error based on the sensitivity matrix. To ensure computability and clarity, we use: wherein is the mapping sensitivity matrix from residual samples to observation-side residual error, and is the parameter sensitivity matrix used in the inversion update

[0040] wherein, is the gravity residual vector organized by stations or sampling points ( is the number of residual samples), is the misclosure residual error vector organized by observation sides ( is the number of misclosure observation sides), is the sensitivity matrix whose elements . Wherein can be organized by station sequence or along the line sampling points, organized by observation sides, and the corresponding relationship between the two is determined by the association rule of the observation side endpoints and residual samples (for example, the weighted combination of the endpoint residuals approximates the error contribution of the side).

[0041] Substitute as the input vector into the traverse network adjustment function to calculate the transverse position deviation projection value and the longitudinal position deviation projection value at the through surface, and define the through error projection value as the norm or weighted norm of the two. One definition is: wherein is the transverse and longitudinal projection deviation at the through surface, is the preset non-negative weight and at least one is positive.

[0042] If ( is the preset threshold), it is determined that the accuracy is qualified and the physical reduction range data is output; otherwise, check whether the iteration termination condition is met. The iteration termination conditions include: the number of iterations reaches the maximum value , the parameter update amount of the adjacent two iterations is less than the convergence threshold , or the adjacent two through error projection values decrease by a preset proportion threshold .

[0043] If the termination condition is not met, perform inversion calculations to update the physical model parameters and return the policy generation and physics reduction. Inversion calculations can be performed using the sensitivity matrix method or gradient descent method. Taking weighted least squares as an example, the objective function is constructed as follows: in, Here, m is the parameter sensitivity matrix, and m is the physical model parameter vector. Dimensions, matrix elements ; This is the weight matrix. The weight matrix... The diagonal elements can be determined by the observation standard deviation and robust weights, for example, taking... ,in For the first The standard deviation of each residual sample For robust weights, the corresponding least squares solution is: And order .

[0044] Select update objects based on dominant component: When the projection value of the dominant penetration error exceeds the limit due to the residual error of the ranging correction corresponding to the high-frequency fluctuation component, the local lithological density parameter corresponding to the spatial location is updated first. ; When the residual error of the ranging correction corresponding to the low-frequency trend component of the extended section of the tunnel exceeds the limit, the regional background density parameter should be updated first. ; When the residual error of the ranging correction corresponding to the low-frequency trend component of the vertical shaft transfer section exceeds the limit, the vertical gradient correction parameters should be updated first. .

[0045] When the termination condition is met but still At that time, the physical distance measurement data with residual error alarm label is output, and the alarm label includes at least the over-limit type, the over-limit amplitude and the corresponding measurement segment label.

[0046] Example 2 A gravity anomaly correction system for high-altitude tunnel ranging, employing the aforementioned gravity anomaly correction method for high-altitude tunnel ranging, includes: a data initialization module, a strategy generation and calculation module, a residual and accuracy evaluation module, a closed-loop control and output module, and a data processing module. The data initialization module, strategy generation and calculation module, residual and accuracy evaluation module, and closed-loop control and output module are respectively connected to the data processing module. The data initialization module is configured to acquire original ranging observation data and station information, determine a section type, and construct input data to be corrected under potential difference or potential height constraints. The strategy generation and reduction module is configured to perform gravity field forward calculation based on physical model parameters to generate a shaft axis gravity distribution sequence or a tunnel axis deflection distribution sequence, and perform physical reduction on geometric observation data to obtain physical reduction ranging data. The residual and accuracy evaluation module is configured to construct a gravity residual sequence and calculate ranging correction residual error and through error projection values. The closed-loop control and output module is configured to determine whether the through error projection values are less than a preset threshold value; if yes, the physical reduction ranging data is output; if no and an iteration termination condition is not met, the physical model parameters are updated and a new round of calculation is triggered; if no and the iteration termination condition is met, the physical reduction ranging data containing an alarm identifier is output.

[0047] Embodiment Three For the calculation of near-zone terrain mass gravity contribution and far-zone Bouguer plate mass gravity contribution, the integration domain of the near zone and the far zone is mutually exclusive in this embodiment, and an adaptive boundary radius is introduced , which is driven by the spectral characteristics of the gravity residual , to suppress systematic bias caused by repeated counting or missing counting and accelerate closed-loop convergence.

[0048] For each station, a horizontal distance integration domain is established around the station, and is defined as the horizontal projection distance of the station to the voxel. The near-zone voxel set and the far-zone voxel set are defined, and it is ensured that the two sets are mutually exclusive and the union covers the calculation domain. The near zone adopts real DEM voxel integration (density taken ), and the far zone adopts Bouguer plate approximation or layered background model (density taken ), so as to obtain: and calculate the deflection component . In order to update adaptively, this embodiment defines a target function mainly based on residual high-frequency energy: wherein is a weighting coefficient, is the residual component calculated and decomposed under a given . In the closed-loop iteration, while updating , a small step update is made to by one-dimensional search or gradient approximation, so that ​The descending is followed by returning to S2-S4 of Embodiment One to continue iteration.

[0049] Embodiment Four In order to solve the problem of sensitivity to abnormal gravity observation or local abnormal density body in the inversion calculation in Embodiment One, this embodiment performs robust weighting on the gravity residual sequence and abnormal segment isolation on the basis of Embodiment One, reduces the destructive effect of outlier residual on parameter update in the inversion update, and thus improves the stability and availability of the closed-loop iteration. For each station residual The observation standard deviation is given (obtained by evaluating the accuracy of the gravity instrument, the observation time or the repeated observation variance), and the standardized residual is constructed: The Huber type weight function is used to define the weight: wherein is a preset inflection point constant. Thus, the robust weight diagonal matrix is formed. In an implementation mode, the maximum weight matrix is constructed in combination with the observation standard deviation: and the is directly used in the inversion objective function and solution of Embodiment One. At the same time, the abnormal segment is defined: if a plurality of continuous stations satisfy ( is an abnormal weight threshold), the continuous interval is marked as an abnormal segment. The local density parameter in the abnormal segment is allowed to be segmented and parameterized: the abnormal segment is numbered by segment number , the abnormal segment is used in , and the parameter vector can be expanded to: to avoid forcibly explaining the local anomaly with a single which leads to overall divergence. In the propagation stage of the systematic error, the observation edges corresponding to the abnormal segment are given a lower weight to suppress the amplification effect of the abnormal segment on the projection value of the through-going surface, so as to realize the robust closed loop of observation-inversion-adjustment consistency.

[0050] The specific embodiments described above further illustrate the purpose, technical solutions and beneficial effects of the present application. It should be understood that the above description is only a specific embodiment of the present application and is not used to limit the protection scope of the present application. Any modification, equivalent replacement, improvement, etc. made within the spirit and principles of the present application shall be included in the protection scope of the present application.

Claims

1. A gravity anomaly correction method for distance measurement in high-altitude tunnels, characterized in that, include: S1. Obtain the original distance measurement observation data and station spatial location information of the section to be processed, determine whether the section is a vertical shaft transfer section or a horizontal tunnel extension section, calculate the potential number of each station based on the global background gravity field model and the digital elevation model of the survey area, and construct the input data to be corrected for potential difference constraints or potential elevation constraints. S2, obtain the initial physical model parameters and perform gravity field forward modeling calculation on the survey area: generate the gravity distribution sequence of the vertical shaft axis for the vertical shaft transfer section, generate the vertical deviation distribution sequence along the tunnel axis for the horizontal tunnel extension section, and perform vertical depth reduction or length reduction on the input data to be corrected according to the sequence to obtain the physical reduction distance data; S3. Collect the measured gravity observation values ​​at the station, and calculate the difference between them and the theoretical gravity observation values ​​calculated from the physical model parameters to obtain the gravity residual sequence. Calculate the distance measurement correction residual error from the gravity residual sequence, and obtain the penetration error projection value through the traverse network error propagation. S4, determine whether the projection value of the through error is less than the preset threshold; If the output is physical distance measurement data, then if not and the iteration termination condition is not met, the physical model parameters are updated by inversion with the goal of minimizing the gravity residual and S2 is returned until the iteration termination condition is met and physical distance measurement data containing alarm flags is output.

2. The gravity anomaly correction method for high-altitude tunnel ranging according to claim 1, characterized in that, The input data to be corrected includes: Based on the global background gravity field model, an Earth gravity potential reference framework is constructed, and the potential values ​​of each station are calculated. When the measurement section is a vertical shaft transfer section, a potential difference constraint is formed based on the potential values ​​of the upper and lower stations, and the potential difference constraint is associated with the vertical shaft geometric depth observation to generate vertical input data to be corrected. When the measurement section is a horizontal tunnel extension section, the physical elevation of the station relative to the preset reduction reference surface is calculated based on the station potential values, and the reduction ratio constraint is established on the horizontal tunnel geometric length observation using the physical elevation to generate horizontal tunnel input data to be corrected.

3. The gravity anomaly correction method for high-altitude tunnel ranging according to claim 1, characterized in that, The process of obtaining initial physical model parameters and performing forward gravity field calculations on the survey area includes: generating a gravity distribution sequence along the shaft axis for the vertical shaft transfer survey section, and generating a vertical deviation distribution sequence along the tunnel axis for the horizontal tunnel extension survey section, including: The physical model parameters include local lithological density parameters, regional background density parameters, and vertical gradient correction parameters. When the measurement section is a vertical shaft transfer measurement section, the gravitational acceleration at each depth point along the shaft axis is calculated based on the vertical gradient correction parameters and combined with the local lithological density parameters to form the gravity distribution sequence along the shaft axis. When the measurement section is a horizontal tunnel extension measurement section, the gravity contribution of the near-field topographic mass is calculated based on the local lithological density parameters, the gravity contribution of the far-field Bougand plate mass is calculated based on the regional background density parameters, and integral modeling is performed on the gravity contribution to obtain the vertical deviation components at each measurement point along the tunnel axis to form the vertical deviation distribution sequence.

4. The gravity anomaly correction method for high-altitude tunnel ranging according to claim 1, characterized in that, The step of performing vertical depth or length reduction on the input data to be corrected based on the sequence to obtain physically reduced distance data includes: When the measurement section is a vertical shaft transfer measurement section, the gravity path integral is performed on the vertical shaft axis based on the gravity distribution sequence of the vertical shaft axis to obtain the theoretical value of the potential difference between the measurement stations, and the vertical shaft geometric depth observation is reduced under the potential difference constraint to obtain the vertical physical reduced distance measurement data. When the measurement section is an extension of a horizontal tunnel, the tilt correction term and the reduction correction term are calculated based on the vertical deviation distribution sequence. The geometric length observation of the horizontal tunnel is corrected and reduced to a preset reference surface to obtain the physical reduction distance data of the horizontal tunnel.

5. The gravity anomaly correction method for high-altitude tunnel ranging according to claim 1, characterized in that, The measured gravity observations include at least one of gravitational acceleration scalar values ​​and gravity vector component observations; the theoretical gravity observations are the corresponding scalar values ​​or corresponding component values ​​in the same coordinate base as the measured gravity observations; the gravity residual sequence is constructed from the difference between the measured gravity observations and the theoretical gravity observations.

6. The gravity anomaly correction method for high-altitude tunnel ranging according to claim 1, characterized in that, The step of calculating the distance correction residual error from the gravity residual sequence and obtaining the penetration error projection value through traverse network error propagation includes: The gravity residual sequence is decomposed to obtain high-frequency fluctuation components and low-frequency trend components; a mapping sensitivity matrix from residual samples to ranging observation edges is constructed to map the high-frequency fluctuation components as random disturbance terms and the low-frequency trend components as system disturbance terms, and the ranging correction residual error affecting the physical reduction ranging data is calculated; the ranging correction residual error is substituted as an input vector into the traverse network adjustment or error propagation function to calculate the lateral position deviation projection value and the longitudinal position deviation projection value at the connection surface, and the connection error projection value is determined accordingly.

7. A gravity anomaly correction method for high-altitude tunnel ranging according to claim 3, characterized in that, When the measurement section is a horizontal tunnel extension section, the computational domain of the integral modeling is divided into a near-field integral domain and a far-field integral domain, and an adaptive boundary radius is introduced. During the closed-loop iteration process, the adaptive boundary radius is updated according to the spectral characteristics of the gravity residual sequence to suppress the duplicate or missed calculation of the gravity contribution of the near-field terrain mass and the gravity contribution of the far-field Bouguer plate mass.

8. The gravity anomaly correction method for high-altitude tunnel ranging according to claim 1, characterized in that, Also includes: Robust weighting of gravity residual sequence and identification of outlier segments: Gravity residuals are standardized based on the observation standard deviation of gravity observations, and robust weights are obtained by using a piecewise weighting function; When the robust weight of a continuous station is lower than the anomaly threshold, the corresponding interval is marked as an anomaly segment; when updating the physical model parameters during inversion, the objective function is weighted based on the robust weight, and the local lithological density parameters corresponding to the anomaly segment are updated using piecewise parameterization.

9. A gravity anomaly correction method for high-altitude tunnel ranging according to claim 1, characterized in that, The iteration termination conditions include: The iteration count reaches the preset maximum value, the update amount of physical model parameters between two adjacent iterations is less than the convergence threshold, or the decrease in the projection value of the breakthrough error between two adjacent iterations is less than the preset proportion threshold.

10. A gravity anomaly correction method for high-altitude tunnel ranging according to claim 1, characterized in that, The process of acquiring the original distance measurement observation data and station spatial location information of the section to be processed, and determining whether the section is a vertical shaft transfer section or a horizontal tunnel extension section, includes: The average slope threshold of the edge is used for determination: when If the slope exceeds the preset threshold, it is determined to be a vertical shaft transfer measurement section; otherwise, it is determined to be a horizontal tunnel extension measurement section. in For the station Slope distance observation between , for The elevation difference between the two stations.

11. A gravity anomaly correction method for high-altitude tunnel ranging according to claim 4, characterized in that, When the measurement section is a vertical shaft transfer measurement section, the theoretical value of the inter-station potential difference is obtained by performing gravity path integration on the vertical shaft axis based on the gravity distribution sequence of the vertical shaft axis, and the vertical shaft geometric depth observation is reduced under the potential difference constraint to obtain vertical physical reduced distance data, including: The theoretical value of the inter-station potential difference is obtained by performing gravity path integration on the shaft axis based on the gravity distribution sequence. : ； in The number of discrete depth points along the shaft axis. , For the first Each depth point; in the vertical shaft survey section, according to the upper station → Next stop Order definition and order and Take the positive value accumulated in the same direction; if the order of the endpoints does not meet the above convention, use... and Involved in ratio calculation; Under potential difference constraints, geometric depth observations The vertical physical distance measurement data was obtained by performing the reduction calculation. : ； in For adjacent stations Potential difference constraint.

12. The gravity anomaly correction method for high-altitude tunnel ranging according to claim 4, characterized in that, When the measured section is an extension of a horizontal tunnel, the tilt correction term and the reduction correction term are calculated based on the vertical deviation distribution sequence. The observed geometric length of the horizontal tunnel is corrected and reduced to a preset reference surface to obtain the physical reduction distance data of the horizontal tunnel, including: The tilt correction term and reduction correction term are calculated using the vertical deviation distribution sequence. The observed geometric length of the horizontal tunnel is corrected and reduced to a preset reference surface. The projection angle of the vertical deviation on the side azimuth is calculated. ； in , ; Define the equivalent zenith angle under the normal basis: ； The geometric horizontal length under the normal basis is obtained as follows: ； Then, based on the average physical elevation of the two ends of the measuring station Reduce the length to a preset reference plane: ； in This is the equivalent radius of curvature or the Earth's average radius constant corresponding to the preset reference surface.

13. A gravity anomaly correction system for distance measurement in high-altitude tunnels, characterized in that, The gravity anomaly correction method for high-altitude tunnel ranging according to any one of claims 1-12 includes: a data initialization module, a strategy generation and calculation module, a residual and accuracy evaluation module, a closed-loop control and output module, and a data processing module. The data initialization module, strategy generation and calculation module, residual and accuracy evaluation module, and closed-loop control and output module are respectively connected to the data processing module. The data initialization module is used to acquire the original distance measurement observation data and station information, determine the measurement segment type, and construct the input data to be corrected under the constraints of potential difference or potential elevation. The strategy generation and reduction module is used to perform forward modeling of gravity field based on physical model parameters to generate a gravity distribution sequence of the shaft axis or a vertical deviation distribution sequence of the tunnel axis, and to perform physical reduction on the geometric observation data to obtain physical reduction distance data. The residual and accuracy evaluation module is used to construct the gravity residual sequence and calculate the distance measurement correction residual error and the projection value of the penetration error. The closed-loop control and output module is used to determine whether the projection value of the penetration error is less than a preset threshold; if so, it outputs the physical distance measurement data; if not and the iteration termination condition is not met, it updates the physical model parameters and triggers a new round of calculation; if not and the iteration termination condition is met, it outputs the physical distance measurement data containing alarm flags.

14. A gravity anomaly correction method for high-altitude tunnel ranging according to claim 13, characterized in that, The step of inverting and updating the physical model parameters if the iteration termination condition is not met includes: Select update objects based on dominant component: When the projection value of the dominant penetration error, which corresponds to the distance correction residual error of the high-frequency fluctuation component, exceeds the limit, the local lithological density parameter corresponding to the spatial location is updated. ; When the residual error of the ranging correction corresponding to the low-frequency trend component of the extended section of the tunnel exceeds the limit, the background density parameter of the updated area is updated. ; When the residual error of the ranging correction corresponding to the low-frequency trend component of the vertical shaft transfer section exceeds the limit, the vertical gradient correction parameters are updated. ; When the termination condition is met but still At that time, output the physical distance measurement data with residual error alarm flag. This is the projection threshold for the through-hole error.

15. A gravity anomaly correction method for high-altitude tunnel ranging according to claim 13, characterized in that, The alarm identifier includes at least the over-limit type, the over-limit amplitude, and the corresponding measurement segment identifier.

Images