Multi-array geomagnetic underground pipeline detection system and method

Abstract

The application discloses a multi-array geomagnetic underground pipeline detection system and method. The system comprises one or more geomagnetic sensor arrays, each array comprising a group of three-axis fluxgate sensors, an analog-to-digital converter (ADC), a geomagnetic background field cancellation signal generation circuit and an analog differential summation node. For a first time sequence segment, the ADC generates first digital geomagnetic time sequence data with a first quantization bit width; the cancellation signal generation circuit generates an analog cancellation signal according to prior background field data, and cancels a known strong magnetic interference component from subsequent analog signals at the analog differential summation node; then the ADC is dynamically adjusted to a second quantization bit width smaller than the first quantization bit width; and the equivalent magnetic moment vector data of a weak magnetic target body is determined by calculating a scaling factor corresponding to a least square error in combination with an amplitude reference benchmark in the first digital geomagnetic time sequence data.

Description

Technical Field

[0001] This invention belongs to the field of pipeline detection technology, and specifically relates to a multi-array geomagnetic underground pipeline detection system and method. Background Technology

[0002] As a key component of urban infrastructure, the accurate detection of underground pipelines is crucial for preventing construction accidents and ensuring the safe operation of cities. Magnetic detection utilizes the magnetic anomaly field generated by the Earth's magnetic field to locate underground ferrous pipelines. As a passive detection technology, it is unaffected by non-magnetic interfaces and can operate effectively even in harsh environments, thus gaining widespread application in the field of underground pipeline detection.

[0003] Existing geomagnetic underground pipeline detection systems typically consist of fluxgate sensors, signal conditioning circuits, analog-to-digital converters (ADCs), and digital signal processing units. A three-axis fluxgate sensor simultaneously measures the components of the geomagnetic field in three orthogonal directions. The output analog voltage signal is conditioned and amplified, then quantized by the ADC, and finally processed by the digital signal processor for magnetic anomaly signal extraction and inversion calculation. However, existing technologies still face the following pressing problems in the practical application of multi-array geomagnetic underground pipeline detection:

[0004] Strong magnetic background interference significantly obscures signals from weak magnetic targets. Large ferromagnetic structures widely distributed in the detection area (such as manhole covers, reinforced concrete structures, and existing large culverts) typically generate magnetic anomaly signals several orders of magnitude stronger than those from deeply buried small-diameter metal pipelines. When a strong magnetic interference object is spatially close to the target pipeline, the strong and weak magnetic signals superimpose at the sensor output. The strong magnetic component occupies the majority of the ADC input dynamic range, while the effective amplitude of the weak magnetic target signal is much smaller than the quantization interval corresponding to the least significant bit of the ADC, thus being submerged in quantization noise, preventing the system from effectively detecting the target pipeline.

[0005] The power consumption and data volume issues brought about by high-resolution ADCs limit the continuous operation capability of the system. To improve the sensitivity of weak magnetic signal detection, existing systems generally use 16-bit or 24-bit high-resolution ADCs. However, high-resolution ADCs mean higher sampling bit width, greater data throughput, and higher power consumption. In multi-array detection systems, the massive amount of data generated by multi-channel synchronous sampling places extremely high demands on storage and transmission bandwidth, and high power consumption severely limits the endurance of battery-powered mobile detection platforms. If a low-resolution ADC is simply used to reduce power consumption, the occupation of the dynamic range by strong magnetic interference will cause the weak magnetic signal to be completely lost, and the system will fall into a dilemma of "high resolution leads to excessive power consumption, while low resolution leads to detection failure". The dynamic changes of the geomagnetic background field interfere with the accurate extraction of magnetic anomaly signals. The geomagnetic diurnal variation can cause magnetic field fluctuations of several nanoteslas to tens of nanoteslas within a day, which is comparable in magnitude to the magnetic anomaly signals generated by deeply buried small-diameter pipelines. Existing systems usually perform differential processing on the background field baseline acquired at the beginning of the detection, but the geomagnetic diurnal variation is continuous, and the initial baseline is no longer accurate in subsequent detections. The background field estimation error introduced by this becomes an important source of noise. Furthermore, changes in the probe's attitude cause the sensor coordinate system to rotate relative to the geographic coordinate system, which in turn alters the background field components, further complicating the interference. Summary of the Invention

[0006] The purpose of this invention is to overcome the defects in the existing technology and provide a multi-array geomagnetic underground pipeline detection system and method.

[0007] One or more geomagnetic sensor arrays, wherein the geomagnetic sensor arrays are fixed on the detection carrier in a three-dimensional spatial arrangement and are configured to continuously acquire, passively, simulate voltage signals of magnetic anomaly fields formed by abnormal changes in the Earth's magnetic field caused by differences in magnetic susceptibility of one or more underground targets in the detection area of ​​the detection system. Each geomagnetic sensor array includes:

[0008] The three-axis fluxgate sensor group is used to synchronously sense the changes in magnetic field strength of the geomagnetic field in three orthogonal axes and output the simulated voltage signal of the magnetic anomaly field to the subsequent conditioning circuit.

[0009] An analog-to-digital converter (ADC) is configured to perform time-series quantization on the analog voltage signal at a preset sampling rate. For a first time-series segment in the analog voltage signal, the ADC is configured to quantize the first analog voltage signal segment using a first quantization bit width to generate first digital geomagnetic time-series data from the initial sampling period of the first analog voltage signal segment.

[0010] A geomagnetic background field cancellation signal generation circuit is configured to generate a digital cancellation waveform based on pre-stored geomagnetic field spatial distribution background field data or historical detection baseline data of the detection area. The digital cancellation waveform is then converted from digital to analog to form an analog cancellation signal. The amplitude and phase of the analog cancellation signal are set to cancel the strong magnetic interference components generated by one or more large fixed ferromagnetic objects or known dynamic magnetic moving bodies at known spatial locations in the detection area from the magnetic anomaly field analog voltage signal before analog-to-digital conversion.

[0011] An analog differential summing node is provided between the analog signal output terminal of the triaxial fluxgate sensor group and the analog signal input terminal of the analog-to-digital converter (ADC). The analog differential summing node is used to differentially superimpose the analog voltage signal and the analog cancellation signal before the analog voltage signal enters the ADC for quantization.

[0012] For subsequent time-series segments in the analog voltage signal, the geomagnetic sensor array is configured as follows:

[0013] After the differential superposition is completed at the analog differential summation node, and when the cumulative detection mileage reaches a preset detection mileage threshold or the cumulative detection time reaches a preset detection time threshold, the first quantization bit width of the analog-to-digital converter (ADC) is dynamically adjusted to a second quantization bit width that is smaller than the first quantization bit width; the detection mileage threshold is determined based on the spatial distribution density of known strong magnetic interference bodies in the detection area, and the detection time threshold is determined based on the time-varying rate of the diurnal variation of the geomagnetic field in the detection area;

[0014] Using the analog-to-digital converter (ADC) measurement output value at the second quantization bit width, and combining it with the amplitude reference in the first digital geomagnetic time series data obtained during the initial sampling period, the equivalent magnetic moment vector data of one or more weak magnetic targets within the detection area is determined by calculating the scaling factor corresponding to the least square error between the first digital geomagnetic time series data and the magnetic anomaly amplitude reconstructed by the ADC measurement output value.

[0015] Furthermore, each of the geomagnetic sensor arrays is configured to change the system power consumption state by dynamically switching the quantization bit width of the analog-to-digital converter (ADC) from the first quantization bit width to the second quantization bit width; the second quantization bit width is a one-bit zero-crossing quantization bit width, that is, only the voltage polarity sign of the differentially superimposed analog signal is output; each of the geomagnetic sensor arrays is also configured to use the signal amplitude envelope in the first digital geomagnetic time series data obtained during the initial sampling period to nonlinearly scale the magnetic anomaly amplitude corresponding to the polarity sequence output under the one-bit zero-crossing quantization bit width, the nonlinear scaling being achieved by constructing an amplitude mapping relationship between the polarity sequence and the first digital geomagnetic time series data.

[0016] Furthermore, the geomagnetic background field cancellation signal generation circuit and the analog differential summing node are used for:

[0017] After differential superposition of the strong magnetic interference components, a sparse signal reconstruction algorithm is used to eliminate the residual magnetic anomaly components caused by the one or more weak magnetic targets from the ADC measurement output value, so as to re-estimate the magnetic moment tensor data of the one or more large fixed ferromagnetic interference objects or known dynamic magnetic moving objects; for the polarity symbol sequence output by the ADC with the second quantization bit width, before the sparse signal reconstruction algorithm is executed, the polarity symbol sequence of multiple consecutive sampling frames is accumulated in the time domain, and the proportion of the number of sampling points with logic high level in the accumulation window to the total number of sampling points in the window is used as the amplitude observation value replacement amount in the time window; by sliding the accumulation window, an amplitude observation value sequence is generated as the input observation vector of the sparse signal reconstruction algorithm to compensate for the amplitude information loss caused by 1-bit quantization;

[0018] The equivalent magnetic moment vector data and the magnetic moment tensor data of the weak magnetic target are spatially superimposed to generate a comprehensive magnetic distribution dataset of the target in the detection area. The comprehensive magnetic distribution dataset contains the classification identifiers and spatial separation information of the weak magnetic target and the large fixed ferromagnetic interference object.

[0019] Furthermore, it also includes:

[0020] A central microcontroller and an embedded digital signal processing unit are coupled via a data bus to the digital output port of the analog-to-digital converter (ADC) of each of the one or more geomagnetic sensor arrays to receive one or more of the first digital geomagnetic time series data or the equivalent magnetic moment vector data.

[0021] And each of the one or more geomagnetic sensor arrays further includes:

[0022] A local solid-state storage unit is configured to cache the first digital geomagnetic time series data from the initial sampling period as reference data for subsequent amplitude recovery and background field update. Before being used as the reference frame for subsequent amplitude recovery, the first digital geomagnetic time series data stored in the local solid-state storage unit is preprocessed by the central microcontroller and embedded digital signal processing unit. The reference frame preprocessing includes: using the digital cancellation waveform corresponding to the analog cancellation signal provided by the geomagnetic background field cancellation signal generation circuit, processing the first digital geomagnetic time series data in the digital domain... A digital geomagnetic time series data is subjected to synchronous cancellation processing to remove the strong magnetic interference component corresponding to the known strong magnetic interference body from the first digital geomagnetic time series data, obtaining a pure reference signal containing only the geomagnetic background field and the magnetic anomaly component of the weak magnetic target body; the pure reference signal serves as the amplitude reference and the observed value in the least squares error calculation, used to determine the scaling factor; the digital cancellation waveform used in the digital domain synchronous cancellation processing corresponds precisely to the analog cancellation signal used at the analog differential summation node in terms of time axis and amplitude, and both originate from the same digital cancellation waveform data source; and

[0023] Geomagnetic interference waveform synthesis circuit, wherein the geomagnetic interference waveform synthesis circuit is configured to:

[0024] The system receives pre-calculated background field data of the geomagnetic field spatial distribution of the detection area from one or more of the local solid-state storage unit or the central microcontroller and embedded digital signal processing unit. The data includes the spatial coordinates of known large ferromagnetic structures and the spatial distribution function of the theoretical magnetic anomaly field they generate.

[0025] Based on the theoretical magnetic anomaly field spatial distribution function, the cached first digital geomagnetic time series data is subjected to time-domain waveform reconstruction and amplitude mapping processing to generate the digital cancellation waveform. The digital cancellation waveform is then converted from digital to analog to form the analog cancellation signal.

[0026] Furthermore, the geomagnetic interference waveform synthesis circuit includes:

[0027] A digital waveform generation and analog-to-digital conversion driving module is configured to generate multi-channel digital cancellation waveform data in real time based on the theoretical magnetic anomaly field spatial distribution function and the current spatial pose information of the sensor array.

[0028] A digital-to-analog converter (DAC) includes a digital input port coupled to the digital waveform generation and digital-to-analog conversion drive module and an analog output port configured to convert the multi-channel digital cancellation waveform data into a continuous analog cancellation voltage signal.

[0029] The analog differential summing node includes a first differential input port for receiving the original analog voltage signal output by the triaxial fluxgate sensor group, a second differential input port coupled to the analog output port of the digital-to-analog converter (DAC), and a differential signal output port for outputting the difference between the original analog voltage signal and the continuous analog cancelling voltage signal.

[0030] Furthermore, the analog-to-digital converter (ADC) includes an analog sampling input coupled to the differential signal output port of the analog differential summing node and a digital data output coupled to the central microcontroller and the embedded digital signal processing unit.

[0031] The present invention also provides a method for identifying multi-array geomagnetic underground pipelines, the method comprising:

[0032] At the triaxial fluxgate sensor of each geomagnetic sensor array in the multi-array geomagnetic underground pipeline detection system, the simulated voltage signal of the magnetic anomaly field generated by one or more underground targets in the detection area of ​​the geomagnetic sensor array due to the difference in magnetic susceptibility and superimposed on the Earth's background magnetic field is continuously received in a passive manner.

[0033] The analog-to-digital converter (ADC) of the geomagnetic sensor array is used to quantize the simulated voltage signal of the magnetic anomaly field within the first time window with a first quantization bit width to generate first digital geomagnetic time series data, and the data is stored in a local reference data memory.

[0034] The first digital geomagnetic time series data stored in the local reference data memory is subjected to digital domain cancellation processing. The magnetic field component contributed by the known strong magnetic interference body is subtracted from the first digital geomagnetic time series data using the digital cancellation waveform generated by the geomagnetic background field cancellation signal generation circuit, to obtain a pure weak magnetic target body plus background field reference signal. The pure reference signal replaces the original first digital geomagnetic time series data as the observation value for solving the amplitude scaling factor in the subsequent steps.

[0035] Using the geomagnetic background field cancellation signal generation circuit of the geomagnetic sensor array, based on the pre-loaded geomagnetic background grid model of the detection area and the list of known strong magnetic interference bodies, combined with the low-frequency trend term in the first digital geomagnetic time series data, a digital cancellation waveform for canceling strong magnetic interference components is determined, and an analog cancellation signal is generated via digital-to-analog conversion.

[0036] Based on the simulated cancellation signal, the simulated differential summation node of the geomagnetic sensor array is used to remove the strong magnetic interference voltage component caused by the one or more large fixed ferromagnetic interference objects or known dynamic magnetic moving bodies from the real-time magnetic anomaly field simulated voltage signal in the subsequent time window to obtain the differential mode simulated signal.

[0037] After the removal is completed at the analog differential summation node, the quantization bit width of the analog-to-digital converter (ADC) is dynamically adjusted from the first quantization bit width to a second quantization bit width that is smaller than the first quantization bit width.

[0038] The differential analog signal is quantized using the analog-to-digital converter (ADC) with the second quantization bit width to generate a time-series differential digital output. Combined with the clean reference signal in the local reference data memory, an objective optimization function is constructed with the clean reference signal as the observation and the symbol sequence of the time-series differential digital output as the prior constraint. The amplitude scaling coefficient in the objective optimization function is solved to determine the equivalent magnetic moment data and spatial location data corresponding to one or more weak magnetic targets in the detection area.

[0039] Further, a signal intensity scaling factor is calculated between the first digital geomagnetic time series data and the inverted magnetic moment data corresponding to the one or more weak magnetic targets. The scaling factor is obtained by solving the least squares error between the first digital geomagnetic time series data and the signal reconstructed by scaling the time series differential digital output using the scaling factor.

[0040] Furthermore, after estimating the equivalent magnetic moment data of the one or more weak magnetic targets, based on the estimated spatial distribution data of the weak magnetic targets, the contribution of the magnetic anomaly field theoretically generated by the weak magnetic targets is subtracted from the first digital geomagnetic time series data to separate the residual geomagnetic signal component that mainly includes the contribution of large fixed ferromagnetic interference objects (the residual geomagnetic signal component may still contain geomagnetic background field components that have not been completely eliminated). In subsequent processing, the geomagnetic background field components are removed as known interference terms through least squares fitting, thereby estimating the magnetic moment tensor data of the large fixed ferromagnetic interference objects.

[0041] The magnetic moment tensor data is spatially registered and numerically superimposed with the equivalent magnetic moment data corresponding to the one or more weakly magnetic targets to generate a complete magnetic source distribution map of all magnetic targets in the detection area.

[0042] Further, estimating the magnetic moment tensor data of the large fixed ferromagnetic interference object includes: after the simulated cancellation signal has eliminated the first-order major contribution of the large fixed ferromagnetic interference object, using the residual higher-order components in the first digital geomagnetic time series data, and using the least squares fitting method to back-calculate the higher-order magnetic multipole moment parameters of the large fixed ferromagnetic interference object from the time series difference digital output of the subsequent time window.

[0043] Furthermore, before determining the digital cancellation waveform, the method includes determining the geomagnetic background grid model of the detection area and the spatial location list of the known strong magnetic interference bodies based on one or more of the following: historical geomagnetic diurnal variation observation data, current array spatial pose data provided by the inertial navigation system, or a geographic information system map of urban underground infrastructure corresponding to the detection area. The geographic information system map of urban underground infrastructure includes the design coordinates of existing large ferromagnetic culverts, manhole covers, reinforced concrete structures, and subway tunnels, as well as their theoretical magnetic anomaly feature fingerprint database.

[0044] Further, determining the digital cancellation waveform includes:

[0045] The first digital geomagnetic time series data is stored in the local solid-state storage unit of the geomagnetic sensor array as a background field calibration baseline;

[0046] At the geomagnetic interference waveform synthesis circuit of the geomagnetic sensor array, the geomagnetic field background grid model of the detection area and the spatial location list of the known strong magnetic interference bodies are received from one or more of the local solid-state storage unit or the central microcontroller and the embedded digital signal processing unit.

[0047] Based on the geomagnetic field background grid model of the detection area and the list of spatial locations of the known strong magnetic interference objects, the background field calibration baseline is subjected to spatial coordinate transformation and magnetic field attenuation function weighting to reconstruct a theoretical interference waveform that matches the current sensor array attitude in the digital domain. The analog cancellation signal is then generated via digital-to-analog conversion to eliminate the strong magnetic interference components corresponding to one or more large fixed ferromagnetic interference objects in the detection area in the analog domain, and simultaneously eliminate the time-varying magnetic interference components corresponding to magnetic moving objects with known motion trajectories in the detection area.

[0048] Furthermore, the second quantization bit width is a one-bit zero-crossing quantization bit width, and the first quantization bit width is a sixteen-bit or higher quantization bit width; the quantization bit width of the analog-to-digital converter (ADC) is dynamically adjusted on the time axis so that the first quantization bit width is used to collect mixed data of background field and target field in the early stage of detection, and the second quantization bit width is used to collect differential polarity data after interference cancellation in the subsequent continuous scanning process.

[0049] The beneficial effects are as follows:

[0050] The geomagnetic background field cancellation signal generation circuit generates a digital cancellation waveform based on pre-stored geomagnetic field spatial distribution background field data or historical detection baseline data of the detection area. After digital-to-analog conversion, an analog cancellation signal is formed. At the analog differential summation node, the magnetic anomaly field analog voltage signal and the analog cancellation signal are differentially superimposed. This ensures that before the analog-to-digital converter (ADC) quantizes, the strong magnetic interference component generated by a large fixed ferromagnetic interference object or a known dynamic magnetic moving body at a known spatial location in the detection area is canceled in the analog domain. This avoids the strong magnetic interference signal occupying most of the quantization dynamic range of the ADC during the quantization process, overcomes the technical defect that the weak magnetic target signal is submerged by quantization noise due to strong magnetic interference, and effectively improves the detection sensitivity of weak magnetic underground pipelines.

[0051] After differential superposition is completed at the analog differential summing node, the first quantization bit width of the analog-to-digital converter (ADC) is dynamically adjusted to a second quantization bit width smaller than the first quantization bit width. The ADC output value at the second quantization bit width is then used to measure the output value. Combined with the amplitude reference in the first digital geomagnetic time series data obtained during the initial sampling period, the scaling factor corresponding to the least square error between the first digital geomagnetic time series data and the magnetic anomaly amplitude reconstructed by the ADC measurement output value is calculated. This determines the equivalent magnetic moment vector data of the weak magnetic target in the detection area. After the strong magnetic interference is canceled by the analog domain, the amplitude information of the weak magnetic target signal can be acquired by the ADC with a low quantization bit width. The accurate magnetic anomaly amplitude can be recovered using the first digital geomagnetic time series data as the amplitude calibration reference. This reduces the quantization bit width of the ADC while maintaining the quantitative inversion accuracy of the magnetic parameters of the weak magnetic target.

[0052] By dynamically switching the quantization bit width of the analog-to-digital converter (ADC) from the first quantization bit width to the second quantization bit width, and the second quantization bit width is a one-bit zero-crossing quantization bit width, the ADC only outputs the voltage polarity sign of the differentially superimposed analog signal during subsequent sampling after the analog domain differential superposition is completed. This significantly reduces the amount of data at each sampling moment, thereby reducing the data storage and transmission bandwidth requirements.

[0053] After the geomagnetic background field cancellation signal generation circuit and the analog differential summing node complete the differential superposition of the strong magnetic interference components, the residual magnetic anomaly components caused by the one or more weak magnetic targets are eliminated from the ADC measurement output value using a sparse signal reconstruction algorithm. This allows for the re-estimation of the magnetic moment tensor data of the one or more large fixed ferromagnetic interference objects or known dynamic magnetic moving objects. The equivalent magnetic moment vector data of the weak magnetic target is then spatially superimposed with the magnetic moment tensor data. This enables the magnetic contributions of the strong magnetic interference object and the weak magnetic target in the same detection area to be separated and reconstructed in stages and at different levels. This achieves the classification, identification, and spatial separation of the strong magnetic interference object and the weak magnetic target, solving the problem of the strong magnetic interference object doubly obscuring the adjacent weak magnetic target in both the spatial and amplitude domains, and improving the spatial resolution capability of underground pipelines.

[0054] The first digital geomagnetic time series data from the initial sampling period is cached in the local solid-state storage unit as a reference frame for subsequent amplitude recovery and background field update. The geomagnetic interference waveform synthesis circuit receives pre-calculated geomagnetic spatial distribution background field data of the detection area from the local solid-state storage unit or the central microcontroller and embedded digital signal processing unit. The cached first digital geomagnetic time series data is reconstructed in the time domain and the amplitude is mapped according to the theoretical magnetic anomaly field spatial distribution function to generate the digital cancellation waveform. This makes the generation of the analog cancellation signal not only dependent on the preset background field model, but also integrates the first digital geomagnetic time series data measured in the current detection period, thereby realizing the dynamic tracking and adaptive update of the analog cancellation signal to the real-time environmental magnetic field changes. Attached Figure Description

[0055] The following figures are for illustrative purposes only and are not intended to limit the scope of the invention, wherein:

[0056] Figure 1 : A flowchart of the system in this invention;

[0057] Figure 2 : A schematic diagram of the working process of the geomagnetic sensor array in this invention;

[0058] Figure 3 : A flowchart illustrating the method in this invention. Detailed Implementation

[0059] To make the objectives, technical solutions, design methods, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are for illustrative purposes only and are not intended to limit the invention.

[0060] like Figure 1 As shown, the multi-array geomagnetic underground pipeline detection system provided by the present invention generally includes one or more geomagnetic sensor arrays and a central microcontroller and embedded digital signal processing unit coupled to the geomagnetic sensor arrays; Figure 2 The internal structure of a single geomagnetic sensor array and its signal flow are shown; Figure 3 The overall flow sequence of the detection method of the present invention is shown.

[0061] This invention provides a multi-array geomagnetic underground pipeline detection system, comprising: one or more geomagnetic sensor arrays, which are fixed on a detection carrier in a three-dimensional spatial arrangement and configured to continuously acquire, passively, simulated voltage signals of magnetic anomaly fields formed by abnormal changes in the Earth's magnetic field caused by differences in magnetic susceptibility of one or more underground targets in the detection area of ​​the system. Each geomagnetic sensor array includes: a triaxial fluxgate sensor group, which is used to synchronously sense changes in the magnetic field strength of the Earth's magnetic field along three orthogonal axes and output the simulated voltage signals of the magnetic anomaly fields to a downstream device. The circuit includes a conditioning circuit; an analog-to-digital converter (ADC) configured to perform time-series quantization on the analog voltage signal at a preset sampling rate; for a first time-series segment of the analog voltage signal, the ADC is configured to quantize the first analog voltage signal segment using a first quantization bit width to generate first digital geomagnetic time-series data from the initial sampling period of the first analog voltage signal segment; and a geomagnetic background field cancellation signal generation circuit configured to generate digital data based on pre-stored geomagnetic field spatial distribution background field data of the detection area or historical detection baseline data. The digital cancellation waveform is converted into an analog cancellation signal after digital-to-analog conversion. The amplitude and phase of the analog cancellation signal are set to cancel the strong magnetic interference components generated by one or more large fixed ferromagnetic interference objects or known dynamic magnetic moving bodies at known spatial locations in the detection area from the magnetic anomaly field analog voltage signal before analog-to-digital conversion. An analog differential summing node is also included between the analog signal output of the triaxial fluxgate sensor group and the analog signal input of the analog-to-digital converter (ADC). This analog differential summing node is used to convert the analog voltage signal into a signal before it enters the ADC for quantization. The signal is differentially superimposed with the analog cancellation signal; wherein, for subsequent time series segments in the analog voltage signal, the geomagnetic sensor array is configured to: after the differential superposition is completed at the analog differential summation node and when the cumulative detection mileage reaches a preset detection mileage threshold or the cumulative detection time reaches a preset detection time threshold, dynamically adjust the first quantization bit width of the analog-to-digital converter (ADC) to a second quantization bit width smaller than the first quantization bit width; the detection mileage threshold is determined based on the spatial distribution density of known strong magnetic interference bodies in the detection area, and the detection time threshold is determined based on the time-varying rate of the diurnal variation of the geomagnetic field in the detection area;Using the analog-to-digital converter (ADC) measurement output value at the second quantization bit width, and combining it with the amplitude reference in the first digital geomagnetic time series data obtained during the initial sampling period, the equivalent magnetic moment vector data of one or more weakly magnetic targets within the detection area is determined by calculating the scaling factor corresponding to the least square error between the first digital geomagnetic time series data and the magnetic anomaly amplitude reconstructed from the ADC measurement output value.

[0062] The triaxial fluxgate sensor assembly includes three orthogonally arranged fluxgate probes mounted on the detection carrier along the X-axis, Y-axis, and Z-axis directions in a Cartesian coordinate system, respectively. The X-axis is parallel to the direction of travel of the detection carrier, the Y-axis is parallel to the lateral direction of the detection carrier, and the Z-axis is parallel to the direction of gravity. Each fluxgate probe independently outputs an analog voltage signal corresponding to its axis, and the voltage amplitude of the analog voltage signal is linearly related to the magnetic field strength along the corresponding axis. The triaxial fluxgate sensor assembly also includes a signal conditioning circuit that performs low-pass filtering and amplification on the analog voltage signals along the three axes, respectively. The cutoff frequency of the low-pass filter is... The frequency is set to be higher than the highest frequency component of the magnetic anomaly signal generated by the target in the detection area, so as to completely retain the frequency band components of the magnetic anomaly signal while filtering out high-frequency noise. The gain coefficient of the amplification process is dynamically adjusted according to the geomagnetic field background level of the detection area. When the geomagnetic field background level is higher than a preset background threshold, the gain coefficient is set to a first gain value. When the geomagnetic field background level is lower than or equal to the preset background threshold, the gain coefficient is set to a second gain value greater than the first gain value. The three axial analog voltage signals after the low-pass filtering and amplification process are simultaneously output to the analog differential summing node in parallel. The parallel mode ensures that the signals of the three axes are strictly synchronized in time and there is no time delay difference introduced by serial transmission.

[0063] The preset sampling rate of the analog-to-digital converter (ADC) is dynamically determined based on the travel speed of the probe, which is provided by one or more of a GPS receiver or odometer. When the travel speed of the probe is higher than a preset speed threshold, the area of ​​the probe swept by the probe per unit time increases, and the preset sampling rate of the ADC is set to a first sampling rate to ensure that the geomagnetic data sampling density corresponding to each spatial sampling point is not lower than a preset density threshold under high-speed travel conditions. When the travel speed of the probe is lower than or equal to the preset speed threshold, the area of ​​the probe swept by the probe per unit time decreases, and the preset sampling rate of the ADC is set to a second sampling rate lower than the first sampling rate to reduce the amount of data collected under low-speed travel conditions, thereby reducing the storage occupation of the local solid-state storage unit and the transmission load of the data bus. The preset speed threshold is jointly determined based on the minimum spatial scale of the target in the probe area and the preset sampling rate of the ADC, such that when sampling at the first sampling rate, the spatial interval between two adjacent sampling points is less than half of the minimum spatial scale, thus satisfying the requirements of the Nyquist spatial sampling theorem.

[0064] The geomagnetic background field cancellation signal generation circuit includes a waveform lookup table memory and a real-time interpolation unit. The waveform lookup table memory pre-stores theoretical magnetic anomaly waveform sequences corresponding to different spatial grid points in the detection area. The different spatial grid points cover the entire three-dimensional spatial range of the detection area. The theoretical magnetic anomaly waveform sequence stored at each spatial grid point contains the complete waveform of the magnetic field vector generated by all known strong magnetic interference bodies at that grid point as a function of time. The real-time interpolation unit reads the theoretical magnetic anomaly waveform sequences corresponding to the eight spatial grid points adjacent to the real-time spatial coordinates from the waveform lookup table memory based on the real-time spatial coordinates of the current sensor array in the detection area. The eight spatial grid points form a spatial cube surrounding the real-time spatial coordinates. The real-time interpolation unit performs trilinear spatial interpolation on the eight theoretical magnetic anomaly waveform sequences to generate the digital cancellation waveform that matches the real-time spatial coordinates. The trilinear spatial interpolation operation performs linear weighting in three spatial dimensions, and the weighting coefficient of each dimension is determined by the relative position of the real-time spatial coordinates with respect to the adjacent grid points in the corresponding dimension.

[0065] The equivalent magnetic moment vector data is a six-tuple data containing three orthogonal directional components and spatial gradient information of the three orthogonal directional components. The spatial gradient information is obtained by performing a differential operation on the measurement values ​​of multiple triaxial fluxgate sensor groups at different spatial locations in the geomagnetic sensor array. The differential operation includes calculating the difference between the coaxial signals output by the multiple triaxial fluxgate sensor groups at the same time and dividing the difference by the corresponding sensor spatial spacing to obtain the spatial rate of change of the axial magnetic field component in that direction. The six-tuple data is used to uniquely characterize the magnetic dipole source parameters of the weak magnetic target. The magnetic dipole source parameters include the spatial centroid coordinates and equivalent magnetic dipole moment vector of the weak magnetic target. The magnetic dipole source parameters are obtained by substituting the six-tuple data into the spatial distribution equation of the magnetic dipole field and solving the inverse problem of the equation. The spatial distribution equation of the magnetic dipole field is established with the spatial centroid coordinates and equivalent magnetic dipole moment vector as unknowns and the magnetic field vector and its spatial gradient as observed values.

[0066] The amplitude reference benchmark in the first digital geomagnetic time series data is obtained by statistically averaging the amplitudes of all sampling points within the first time series segment. The statistical average includes calculating the magnitude of the triaxial magnetic field strength vector amplitude at each sampling time in the first time series segment, and the arithmetic mean of all magnitude calculation results. The amplitude reference benchmark is used as a normalization factor for amplitude recovery during the reconstruction of the analog-to-digital converter (ADC) measurement output value in subsequent time series segments. The calculation of the least squares error uses the normalization factor as the initial estimate of the scaling factor. An iterative optimization algorithm searches for the scaling factor near the initial estimate that makes the least squares error reach a global minimum. The iterative optimization algorithm can be either gradient descent or Newton's method.

[0067] The phase of the simulated cancellation signal is set to differ from the phase of the strong magnetic interference component in the simulated voltage signal of the magnetic anomaly field by 180 degrees. This ensures that when the simulated cancellation signal and the simulated voltage signal of the magnetic anomaly field are differentially superimposed at the simulated differential summation node, the instantaneous voltage value of the simulated cancellation signal and the instantaneous voltage value of the strong magnetic interference component are opposite in sign and equal in amplitude at every moment, thereby achieving complete cancellation of the strong magnetic interference component. The amplitude of the simulated cancellation signal is normalized according to the midpoint between the upper and lower limits of the input dynamic range of the analog-to-digital converter (ADC), so that the dynamic range of the differential analog signal is symmetrically distributed around the midpoint, thereby fully utilizing the effective quantization range of the ADC at the second quantization bit width.

[0068] Furthermore, each of the geomagnetic sensor arrays is configured to change the system power consumption state by dynamically switching the quantization bit width of the analog-to-digital converter (ADC) from the first quantization bit width to the second quantization bit width; the second quantization bit width is a one-bit zero-crossing quantization bit width, that is, only the voltage polarity sign of the differentially superimposed analog signal is output; each of the geomagnetic sensor arrays is also configured to use the signal amplitude envelope in the first digital geomagnetic time series data obtained during the initial sampling period to nonlinearly scale the magnetic anomaly amplitude corresponding to the polarity sequence output under the one-bit zero-crossing quantization bit width, the nonlinear scaling being achieved by constructing an amplitude mapping relationship between the polarity sequence and the first digital geomagnetic time series data.

[0069] The nonlinear scaling is achieved by constructing an amplitude mapping relationship between the polarity sequence and the first digital geomagnetic time series data. This amplitude mapping relationship is characterized by a piecewise linear mapping function, which maps the run length of consecutive identical signs in the polarity sequence to the corresponding magnetic anomaly amplitude value. The run length is the number of consecutive identical sign bits in the polarity sequence. Each linear segment of the piecewise linear mapping function corresponds to a preset run length interval. Each linear segment is defined by its starting run length, ending run length, and corresponding starting and ending amplitude values. When the run length of a certain run in the polarity sequence falls within the preset run length interval, the piecewise linear mapping function performs linear interpolation between the starting and ending amplitude values ​​based on the relative position of the run length between the starting and ending run lengths to obtain the corresponding magnetic anomaly amplitude value. The starting and ending amplitude values ​​are determined by least-squares fitting of the relationship between the run length of the polarity sequence and the amplitude value at the corresponding time position in the first digital geomagnetic time series data.

[0070] The geomagnetic sensor array also includes a power management unit coupled to the quantization bit width control port of the analog-to-digital converter (ADC). When the ADC operates at the second quantization bit width, the power management unit cuts off the power supply to the high-order comparator array inside the ADC. The high-order comparator array consists of all comparator circuits in the ADC used to achieve a quantization resolution higher than the second quantization bit width. While cutting off the power supply to the high-order comparator array, the power management unit keeps the bias current path of the least significant bit comparator inside the ADC in a conducting state. The least significant bit comparator is a single comparator circuit in the ADC used to implement the quantization operation corresponding to the second quantization bit width. When the ADC switches from the second quantization bit width back to the first quantization bit width, the power management unit restores the power supply to the high-order comparator array. The power restoration operation is performed according to a predetermined power-on sequence, which includes applying a bias voltage first and then applying a bias current to avoid the high-order comparator array generating instantaneous spike current during power-on.

[0071] The polarity sequence with a one-bit zero-crossing quantization width is stored in the local solid-state storage unit in the form of a binary symbol stream. Each symbol bit in the binary symbol stream corresponds to the voltage polarity at a sampling time. When the voltage of the differential analog signal is higher than zero, the symbol bit is encoded as a logic high level; when the voltage of the differential analog signal is lower than zero, the symbol bit is encoded as a logic low level. A predetermined number of adjacent symbol bits form a data frame. The predetermined number is eight, sixteen, or thirty-two. The frame header of each data frame contains timestamp information. The timestamp information is obtained by counting the sampling clock cycles of the analog-to-digital converter (ADC). The counting starts from the beginning of the initial sampling period and is triggered by the rising edge of each sampling clock cycle of the ADC. The timestamp information is used to precisely align the polarity sequence with the first digital geomagnetic time series data on the time axis. The precise alignment is achieved by comparing the timestamp information with the time reference recorded in the first digital geomagnetic time series data, calculating the time deviation between the two, and performing time axis translation correction on the polarity sequence based on the time deviation.

[0072] The amplitude mapping relationship is implemented using a lookup table. The table is stored in the local solid-state storage unit. The index of the table is a symbol combination pattern consisting of a predetermined number of consecutive symbol bits in the polarity sequence. Each entry in the table stores a magnetic anomaly amplitude estimate corresponding to the symbol combination pattern. The magnetic anomaly amplitude estimate is obtained by statistically averaging the amplitude values ​​of all time positions in the first digital geomagnetic time series data that match the symbol combination pattern. When a symbol combination pattern appears in the polarity sequence, the geomagnetic sensor array directly reads the corresponding entry value from the table as the magnetic anomaly amplitude recovery value corresponding to the symbol combination pattern.

[0073] The signal amplitude envelope is obtained by performing a Hilbert transform on the first digital geomagnetic time series data and taking the modulus of the Hilbert transform result. The Hilbert transform transforms the first digital geomagnetic time series data from the real domain to the complex domain. The real part of the transformed complex signal is the original sequence, and the imaginary part is the orthogonal component of the original sequence. The orthogonal component is calculated by convolving the original sequence with the Hilbert transform kernel function. The expression of the Hilbert transform kernel function in the time domain is the negative two times the reciprocal of the time variable divided by pi. The modulus is the arithmetic square root of the sum of the squares of the real and imaginary parts of the complex signal. The instantaneous values ​​of the signal amplitude envelope at each sampling time constitute a slowly changing envelope waveform. The envelope waveform is used to determine the local amplitude reference corresponding to each sampling time in the nonlinear scaling.

[0074] The second quantization bit width also supports switching between the one-bit zero-crossing quantization bit width and the two-bit quantization bit width. When the background magnetic field gradient of the detection area exceeds a preset gradient threshold, the geomagnetic sensor array switches the second quantization bit width from the one-bit zero-crossing quantization bit width to the two-bit quantization bit width. The two-bit quantization bit width outputs four quantization levels, which correspond to the differential-mode analog signal voltage being located in four intervals: less than the negative reference voltage, between the negative reference voltage and zero voltage, between zero voltage and the positive reference voltage, and greater than the positive reference voltage.

[0075] Furthermore, the geomagnetic background field cancellation signal generation circuit and the analog differential summing node are used to: after completing the differential superposition of the strong magnetic interference components, eliminate the residual magnetic anomaly components caused by the one or more weak magnetic targets from the analog-to-digital converter (ADC) measurement output value using a sparse signal reconstruction algorithm, so as to re-estimate the magnetic moment tensor data of the one or more large fixed ferromagnetic interference objects or known dynamic magnetic moving bodies; for the polarity symbol sequence output by the ADC with the second quantization bit width, the sparse signal reconstruction algorithm first performs time-domain analysis on the polarity symbol sequence of multiple consecutive sampling frames before execution. Accumulation: The proportion of sampling points with a logic high level within the accumulation window to the total number of sampling points in the window is used as the amplitude observation value replacement quantity within that time window; by sliding the accumulation window, an amplitude observation value sequence is generated as the input observation vector of the sparse signal reconstruction algorithm to compensate for the loss of amplitude information caused by 1-bit quantization; the equivalent magnetic moment vector data of the weak magnetic target and the magnetic moment tensor data are spatially superimposed to generate a comprehensive magnetic distribution dataset of the target in the detection area, which includes the classification identifier and spatial separation information of the weak magnetic target and the large fixed ferromagnetic interference object.

[0076] The sparse signal reconstruction algorithm is a basis pursuit denoising algorithm. The basis pursuit denoising algorithm uses the measured output value of the analog-to-digital converter (ADC) as the observation vector and a pre-constructed dictionary matrix as the sparse representation basis. Each column of the dictionary matrix corresponds to the theoretical magnetic field response vector generated by a magnetic dipole source with a unit magnetic moment at a spatial grid point in the detection area at the location of the geomagnetic sensor array. The dimension of the theoretical magnetic field response vector is equal to the product of the total number of sensor channels in the geomagnetic sensor array and the number of sampling times. The basis pursuit denoising algorithm obtains the sparse coefficient vector by solving the following optimization problem: minimizing the sum of squares of the residuals between the norm of the sparse coefficient vector and the product of the observation vector and the dictionary matrix multiplied by the sparse coefficient vector, multiplied by the sum of regularization parameters. The regularization parameters are used to balance the weight between sparsity and data fit.

[0077] The comprehensive magnetic distribution dataset of the target objects is stored in the form of a three-dimensional grid map. Each grid node in the three-dimensional grid map records the total magnetic moment vector, classification label, and confidence score at that grid node. The total magnetic moment vector is the sum of the magnetic moment vectors of all magnetic targets at that grid node. The classification label is determined by comparing the amplitude of the total magnetic moment vector with a preset amplitude threshold and the direction of the total magnetic moment vector with a preset direction threshold. When the amplitude of the total magnetic moment vector is greater than the preset amplitude threshold and the angle between the direction of the total magnetic moment vector and the direction of the geomagnetic field of the detection area is greater than a preset angle threshold, the classification label is marked as a large fixed ferromagnetic interference object. When the amplitude of the total magnetic moment vector is less than or equal to the preset amplitude threshold and the angle between the direction of the total magnetic moment vector and the direction of the geomagnetic field of the detection area is less than or equal to the preset angle threshold, the classification label is marked as a weak magnetic target object. The confidence score is obtained by normalizing the signal-to-noise ratio of the signal at the grid node.

[0078] The spatial separation information includes the spatial centroid coordinates of each magnetic target, the spatial boundary of the magnetic target, and the minimum distance between adjacent magnetic targets. The spatial boundary is determined by performing connected component analysis on continuous high-value regions in the comprehensive magnetic distribution dataset. The connected component analysis marks all grid nodes in the 3D grid map with confidence scores higher than a preset confidence threshold as candidate nodes, and groups adjacent nodes among the candidate nodes into the same connected component. Each connected component corresponds to one magnetic target. The spatial centroid coordinates are the weighted centroid coordinates obtained by weighting the spatial coordinates of all grid nodes in the same connected component according to the confidence scores of each grid node. The spatial boundary of the magnetic target is determined by solving the minimum bounding box of the connected component in three orthogonal directions. The minimum distance between adjacent magnetic targets is obtained by calculating the Euclidean distance between any two grid nodes in different connected components and taking the minimum value.

[0079] The sparse signal reconstruction algorithm further includes an iterative refinement step for the magnetic moment tensor data of the large fixed ferromagnetic interference object. After the sparse signal reconstruction algorithm first estimates the magnetic moment tensor data of the large fixed ferromagnetic interference object, the iterative refinement step substitutes the magnetic moment tensor data as a known parameter into the forward model to calculate the theoretical magnetic field response generated by the large fixed ferromagnetic interference object at the geomagnetic sensor array. The theoretical magnetic field response is then subtracted from the measurement output value of the analog-to-digital converter (ADC). The sparse signal reconstruction algorithm is then re-executed using the residual data after subtraction to update the equivalent magnetic moment vector data of the weak magnetic target. The above process is repeated until the relative change in the equivalent magnetic moment vector data of the weak magnetic target between two adjacent iterations is less than a preset convergence threshold.

[0080] The magnetic moment tensor data is a second-order tensor. The nine components of the second-order tensor correspond to the magnetic field components generated by the magnetic moment components in the X-axis direction, Y-axis direction, Z-axis direction, Y-axis direction, Y-axis direction, Z-axis direction, X-axis direction, Y-axis direction, Z-axis direction, X-axis direction, Y-axis direction, Y-axis direction, Z-axis direction, X-axis direction, Y-axis direction, Z-axis direction, Y-axis direction, and Z-axis direction.

[0081] After the analog differential summing node completes the differential superposition of the strong magnetic interference components, the dynamic range of the differential analog signal output to the analog-to-digital converter (ADC) is controlled to not exceed the full-scale input range of the ADC under the second quantization bit width, and preferably less than half of the full-scale input range, so as to ensure that the signal components of the weak magnetic target are not truncated or lost under the second quantization bit width, and at the same time reserve dynamic margin for diurnal variations of the geomagnetic field and residual interference that has not been completely canceled.

[0082] Furthermore, it also includes: a central microcontroller and an embedded digital signal processing unit, the central microcontroller and the embedded digital signal processing unit being coupled to the digital output port of the analog-to-digital converter (ADC) of each of the one or more geomagnetic sensor arrays via a data bus to receive one or more of the first digital geomagnetic time series data or the equivalent magnetic moment vector data; and each of the one or more geomagnetic sensor arrays further includes: a local solid-state storage unit, the local solid-state storage unit being configured to cache the first digital geomagnetic time series data from the initial sampling period as reference data for subsequent amplitude recovery and background field update; before the first digital geomagnetic time series data stored in the local solid-state storage unit is used as the reference frame for subsequent amplitude recovery, the central microcontroller and the embedded digital signal processing unit first perform reference frame preprocessing, the reference frame preprocessing including: using the digital cancellation waveform corresponding to the analog cancellation signal provided by the geomagnetic background field cancellation signal generation circuit, performing synchronous cancellation processing on the first digital geomagnetic time series data in the digital domain, so as to obtain the first digital geomagnetic time series data from the first digital geomagnetic time series data. The strong magnetic interference component corresponding to the known strong magnetic interference body is removed from the inter-series data to obtain a pure reference signal containing only the geomagnetic background field and the magnetic anomaly component of the weak magnetic target body; the pure reference signal is used as the amplitude reference and the observation value in the least square error calculation to determine the scaling factor; the digital cancellation waveform used in the digital domain synchronous cancellation processing corresponds precisely to the analog cancellation signal used at the analog differential summation node in terms of time axis and amplitude, and both originate from the same digital cancellation waveform data source; and a geomagnetic interference waveform synthesis circuit is configured to: receive pre-calculated geomagnetic field spatial distribution background field data of the detection area from one or more of the local solid-state storage unit or the central microcontroller and embedded digital signal processing unit, the data containing the spatial coordinates of the known large ferromagnetic structure and the spatial distribution function of the theoretical magnetic anomaly field it generates; perform time-domain waveform reconstruction and amplitude mapping processing on the cached first digital geomagnetic time series data according to the theoretical magnetic anomaly field spatial distribution function to generate the digital cancellation waveform, which is converted from digital to analog to form the analog cancellation signal.

[0083] The data bus is a high-speed serial peripheral interface bus. The high-speed serial peripheral interface bus connects the one or more geomagnetic sensor arrays in series in a daisy-chain topology to a single bus master port of the central microcontroller and embedded digital signal processing unit. In the daisy-chain topology, each geomagnetic sensor array acts as a slave device on the bus. Each slave device has a unique bus address, which is determined by the cascading position of the slave device in the daisy chain. The cascading position is set by either a hardware coding switch or an automatic address allocation mechanism during power-on initialization. The central microcontroller and embedded digital signal processing unit, as bus masters, sequentially select each slave device for data transmission via a chip select signal line. At any given time, only the selected slave device occupies the data bus for data transmission, while the remaining unselected slave devices place their respective data output ports in a high-impedance state to avoid bus conflicts.

[0084] The local solid-state storage unit includes a dual-port static random access memory (SRAM). The first port of the SRAM is coupled to the digital output port of the analog-to-digital converter (ADC) to directly write the first digital geomagnetic time series data. The second port of the SRAM is coupled to the geomagnetic interference waveform synthesis circuit to provide read access to the reference frame. The first port and the second port each have independent address buses and data buses, so that write operations performed on the first port and read operations performed on the second port can be performed simultaneously without interference. The SRAM also includes an arbitration logic circuit. When the first port and the second port access the same memory address simultaneously, the arbitration logic circuit determines which port's access request is responded to first according to a preset access priority rule. The preset access priority rule is either write operation priority over read operation or read operation priority over write operation.

[0085] The theoretical magnetic anomaly field spatial distribution function is a magnetic dipole field distribution function. The parameters of the magnetic dipole field distribution function include the spatial coordinates of the known large ferromagnetic structure, the equivalent magnetic dipole moment vector, and the magnetic dipole orientation angle. The spatial coordinates of the known large ferromagnetic structure are extracted from the geographic information system map of the urban underground infrastructure in the detection area. The equivalent magnetic dipole moment vector is obtained by multiplying the volume and magnetic susceptibility of the known large ferromagnetic structure by the geomagnetic field strength vector of the detection area. The magnetic dipole orientation angle is determined by the angle between the major axis direction of the known large ferromagnetic structure and the geomagnetic field direction of the detection area.

[0086] The multi-array data fusion operation performed by the central microcontroller and the embedded digital signal processing unit includes spatial coordinate registration, time synchronization, and weighted superposition processing of measurement data from different geomagnetic sensor arrays. The spatial coordinate registration is achieved by transforming the local coordinate system of each sensor array to a unified global geographic coordinate system. The time synchronization is achieved by aligning the timestamps of the sampled data of each sensor array using a common time base reference. The weighted superposition determines the fusion weight of the data from each sensor array based on the signal-to-noise ratio or confidence level of each sensor array.

[0087] After the reference frame stored in the local solid-state storage unit is read by the geomagnetic interference waveform synthesis circuit, it continues to be stored in the local solid-state storage unit until the next initial sampling period arrives. When the next initial sampling period arrives, the reference frame stored in the local solid-state storage unit is overwritten by the newly acquired first digital geomagnetic time series data. The overwritten old reference frame is transmitted to the central microcontroller and embedded digital signal processing unit as historical reference data for permanent storage or archiving before the overwriting operation is performed.

[0088] When the geomagnetic background field cancellation signal generation circuit receives the pre-calculated geomagnetic field spatial distribution background field data, it also receives error correction parameters provided by the central microcontroller and the embedded digital signal processing unit. The error correction parameters are calculated based on the root mean square value of the residual of the differential signal output by the analog differential summing node in the previous sampling period. The root mean square value of the residual is used to adjust the amplitude and phase of the digital cancellation waveform so that the root mean square value of the residual of the differential signal in the current sampling period is less than that in the previous sampling period.

[0089] Further, the geomagnetic interference waveform synthesis circuit includes: a digital waveform generation and analog-to-digital conversion driving module, configured to generate multi-channel digital cancellation waveform data in real time based on the theoretical magnetic anomaly field spatial distribution function and the current spatial pose information of the sensor array; a digital-to-analog converter (DAC), including a digital input port coupled to the digital waveform generation and DAC driving module and an analog output port configured to convert the multi-channel digital cancellation waveform data into a continuous analog cancellation voltage signal; an analog differential summing node, including a first differential input port for receiving the original analog voltage signal output by the triaxial fluxgate sensor group, a second differential input port coupled to the analog output port of the DAC, and a differential mode signal output port for outputting the difference between the original analog voltage signal and the continuous analog cancellation voltage signal; and wherein the analog-to-digital converter (ADC) includes an analog sampling input terminal coupled to the differential mode signal output port of the analog differential summing node and a digital data output terminal coupled to the central microcontroller and the embedded digital signal processing unit.

[0090] The digital waveform generation and analog-to-digital conversion driving module includes a direct digital frequency synthesizer. The direct digital frequency synthesizer takes the frequency control word and phase control word provided by the theoretical magnetic anomaly field spatial distribution function as input, and generates the multi-channel digital cancellation waveform data in real time through phase accumulation and amplitude lookup table operations. The phase accumulation is achieved by accumulating the frequency control word to the current value of the phase accumulator in each analog-to-digital conversion sampling clock cycle. The current value of the phase accumulator is updated once in each sampling clock cycle. The amplitude lookup table operation uses the output value of the phase accumulator as an address index to read the corresponding waveform amplitude value from the waveform lookup table. The waveform lookup table stores a complete cycle of reference waveform sequence in advance.

[0091] The update rate of the digital-to-analog converter (DAC) is synchronized with the preset sampling rate of the analog-to-digital converter (ADC). The DAC outputs an updated analog cancellation voltage value in each analog-to-digital conversion sampling clock cycle, so that the differential signal output by the analog differential summing node corresponds to the instantaneous difference between the original analog voltage signal and the analog cancellation signal at the current sampling time at each sampling moment. The synchronization is achieved by driving the DAC and the ADC with the same clock signal output by the central microcontroller and the embedded digital signal processing unit. The same clock signal is fed to the clock input of the DAC and the clock input of the ADC through a clock distribution network.

[0092] The analog differential summing node is an analog subtractor composed of a fully differential operational amplifier. The fully differential operational amplifier includes a non-inverting input terminal, an inverting input terminal, and a differential output terminal. The non-inverting input terminal serves as the first differential input port to receive the original analog voltage signal, the inverting input terminal serves as the second differential input port to receive the continuous analog cancellation voltage signal, and the differential output terminal serves as the differential signal output port to output the difference value. The fully differential operational amplifier also includes a common-mode feedback circuit. The common-mode feedback circuit detects the common-mode voltage at the differential output terminal and compares the common-mode voltage with a preset common-mode reference voltage. Based on the comparison result, it adjusts the bias current inside the fully differential operational amplifier to stabilize the common-mode voltage, so that the common-mode voltage at the differential output terminal remains near the preset common-mode reference voltage.

[0093] A low-pass reconstruction filter is provided between the output terminal of the digital-to-analog converter (DAC) and the second differential input port of the analog differential summing node. The low-pass reconstruction filter is used to filter out the high-frequency quantization noise component in the continuous analog cancellation voltage signal output by the DAC. The frequency of the high-frequency quantization noise component is higher than half of the Nyquist sampling frequency of the analog-to-digital converter (ADC). The cutoff frequency of the low-pass reconstruction filter is set to be equal to half of the Nyquist sampling frequency of the analog-to-digital converter (ADC).

[0094] The digital waveform generation and analog-to-digital conversion driving module includes multiple output channels. The number of output channels is equal to the number of three-axis fluxgate sensor groups in the geomagnetic sensor array multiplied by three. Each output channel independently generates digital cancellation waveform data corresponding to one sensor axis. The digital cancellation waveform data corresponding to each of the multiple output channels are respectively sent to multiple digital-to-analog converters (DACs) for digital-to-analog conversion. The converted multi-channel analog cancellation voltage signals are sent in parallel to the corresponding multiple analog differential summing nodes.

[0095] The quantization bit width of the digital-to-analog converter (DAC) is greater than or equal to the first quantization bit width of the analog-to-digital converter (ADC) to ensure that the amplitude resolution of the analog cancellation signal is not lower than the amplitude resolution of the original analog voltage signal, so that the differential signal output after differential superposition does not contain additional quantization noise introduced by insufficient quantization accuracy of the DAC.

[0096] This invention also provides a method for identifying underground pipelines in a multi-array geomagnetic system. The method includes: continuously receiving, in a passive manner, simulated voltage signals of magnetic anomalies generated by differences in magnetic susceptibility of one or more underground targets within the detection area of ​​the geomagnetic sensor array and superimposed on the Earth's background magnetic field at the triaxial fluxgate sensor of each geomagnetic sensor array; quantizing the simulated voltage signals of the magnetic anomalies within a first time window using the analog-to-digital converter (ADC) of the geomagnetic sensor array with a first quantization bit width to generate first digital geomagnetic time series data, and storing this data in a local reference data store. The first digital geomagnetic time series data stored in the local reference data memory is subjected to digital domain cancellation processing. Using the digital cancellation waveform generated by the geomagnetic background field cancellation signal generation circuit, the magnetic field component contributed by the known strong magnetic interference body is subtracted from the first digital geomagnetic time series data to obtain a pure weak magnetic target body plus background field reference signal. This pure reference signal replaces the original first digital geomagnetic time series data as the observed value for solving the amplitude scaling factor in subsequent steps. The geomagnetic background field cancellation signal generation circuit of the geomagnetic sensor array is used to calculate the background magnetic field signal based on the pre-loaded geomagnetic background network of the detection area. Based on the lattice model and a list of known strong magnetic interference objects' spatial locations, combined with the low-frequency trend term in the first digital geomagnetic time series data, a digital cancellation waveform for canceling the strong magnetic interference component is determined, and an analog cancellation signal is generated via digital-to-analog conversion. Based on the analog cancellation signal, the strong magnetic interference voltage component caused by the one or more large fixed ferromagnetic interference objects or known dynamic magnetic moving bodies is removed from the real-time magnetic anomaly field analog voltage signal within subsequent time windows using the analog differential summing node of the geomagnetic sensor array, to obtain a differential-mode analog signal. After the analog differential summing node completes the removal, the signal from the analog-to-digital converter (ADC) is... The quantization bit width is dynamically adjusted from the first quantization bit width to a second quantization bit width smaller than the first quantization bit width; the differential analog signal is quantized using the analog-to-digital converter (ADC) with the second quantization bit width to generate a time-series differential digital output; and combined with the clean reference signal in the local reference data memory, an objective optimization function is constructed with the clean reference signal as the observation value and the symbol sequence of the time-series differential digital output as the prior constraint. The amplitude scaling coefficient in the objective optimization function is solved to determine the equivalent magnetic moment data and spatial location data corresponding to one or more weak magnetic targets in the detection area.

[0097] Determining the digital cancellation waveform by combining the low-frequency trend term in the first digital geomagnetic time series data includes: performing low-pass filtering on the first digital geomagnetic time series data to extract the low-frequency trend term, wherein the cutoff frequency of the low-pass filtering is twice the highest frequency component of the diurnal variation of the geomagnetic field in the detection area, the highest frequency component of the diurnal variation of the geomagnetic field corresponds to the periodic variation of the geomagnetic field with the Earth's rotation throughout the day; performing a weighted average of the low-frequency trend term and the theoretical value of the background field at the corresponding spatial location in the geomagnetic field background grid model of the detection area, wherein the weight coefficient of the weighted average is determined based on the signal-to-noise ratio of the first digital geomagnetic time series data, and when the signal-to-noise ratio is higher than a preset signal-to-noise ratio threshold, the weight of the low-frequency trend term is set to be higher than the weight of the theoretical value of the background field; and using the result of the weighted average as the amplitude control parameter of the digital cancellation waveform.

[0098] The objective optimization function includes a data fitting term and a sparse regularization term. The data fitting term measures the difference between the first digital geomagnetic time series data and the signal reconstructed by scaling the time series differential digital output using the amplitude scaling factor. The specific form of the data fitting term is half of the sum of squares of the differences between the first digital geomagnetic time series data and the reconstructed signal. The sparse regularization term is the norm of the amplitude scaling factor, which is the number of non-zero elements in the amplitude scaling factor. The objective optimization function is solved iteratively using the alternating direction multiplier method. The alternating direction multiplier method decomposes the objective optimization function into a subproblem concerning the amplitude scaling factor and a subproblem concerning the auxiliary variable. The two subproblems are solved iteratively and alternately until convergence.

[0099] The spatial location data is obtained by substituting the equivalent magnetic moment data into the magnetic dipole forward model for inversion calculation. The magnetic dipole forward model uses spatial coordinates and magnetic moment vector as unknown parameters and the equivalent magnetic moment data as observed values. It solves for the unknown parameters through nonlinear least squares iteration. The nonlinear least squares iteration adopts the Levenberg-Marquardt algorithm. In each iteration, the algorithm calculates the Jacobian matrix of the objective function with respect to the unknown parameters and constructs an approximate Hessian matrix based on the Jacobian matrix. The update step size of the current iteration is determined by solving a system of linear equations in which the product of the approximate Hessian matrix and the update step size equals the negative gradient direction.

[0100] The length of the first time window is determined based on the spatial distribution density of known strong magnetic interference bodies in the detection area. When the spatial distribution density of the known strong magnetic interference bodies is higher than a preset density threshold, the length of the first time window is set as a first window length. When the spatial distribution density of the known strong magnetic interference bodies is lower than or equal to the preset density threshold, the length of the first time window is set as a second window length that is greater than the first window length.

[0101] The local reference frame memory is a dedicated storage area in the local solid-state storage unit of the geomagnetic sensor array. Before writing the first digital geomagnetic time series data, the dedicated storage area is first erased. The erase operation resets all data bits stored in the dedicated storage area to a logic low level to ensure that the first digital geomagnetic time series data is not mixed with historical data when it is written.

[0102] The time-varying magnetic interference component of the known dynamic magnetic moving body is eliminated by superimposing the time-varying magnetic field component corresponding to the trajectory of the dynamic magnetic moving body into the analog cancellation signal. The time-varying magnetic field component is calculated in real time based on the current spatial position, speed and direction of the dynamic magnetic moving body. The current spatial position, speed and direction of the moving body are provided by one or more of the external traffic monitoring system or the radar speed measuring device on the detection carrier.

[0103] Further, a signal intensity scaling factor is calculated between the first digital geomagnetic time series data and the inverted magnetic moment data corresponding to the one or more weak magnetic targets. The scaling factor is obtained by solving the least squares error between the first digital geomagnetic time series data and the signal reconstructed by scaling the time series differential digital output using the scaling factor.

[0104] The least squares error is calculated using a recursive least squares algorithm. This algorithm iteratively updates the gain vector and covariance matrix at each sampling time to output the current scale factor estimate in real time without needing to store all historical sampling data. The gain vector is obtained by left-multiplying the current covariance matrix by the current input vector and then dividing by one, plus the result of left-multiplying the transpose of the current input vector by the covariance matrix and then left-multiplying by the current input vector. The covariance matrix is ​​updated by subtracting the product of the gain vector and the transpose of the current input vector multiplied by the previous covariance matrix from the covariance matrix of the previous time step, and then dividing by a forgetting factor. The forgetting factor is a constant between zero and one.

[0105] The signal strength scaling factor is calculated independently for each of the three orthogonal axes of the three-axis fluxgate sensor group. Each axis corresponds to an independent signal strength scaling factor. The independent scaling factors of the three axes are used to reconstruct the magnetic moment components of the weak magnetic target in the three orthogonal axes. The independent scaling factors of the three axes each adopt the same least squares error calculation method, but each uses the first digital geomagnetic time series data of its own axis and the time series difference digital output of its own axis as input. The independent scaling factors of the three axes do not affect each other during the calculation process and are updated independently.

[0106] When the change in the scaling factor is less than a preset stability threshold within a predetermined number of consecutive update cycles, the method locks the current scaling factor to a fixed value and pauses the iterative update of the scaling factor to reduce computational overhead. The predetermined number is five, ten, or twenty, and the stability threshold is any one-hundredth or one-thousandth of the current scaling factor value. After the scaling factor is locked, the method continues to monitor the statistical characteristics of the time series differential digital output. When the change in the statistical characteristics exceeds a preset statistical change threshold, the method unlocks the scaling factor and resumes the iterative update of the scaling factor.

[0107] The sliding time window used in the calculation of the least squares error is dynamically adjusted according to the traveling speed of the probe. When the traveling speed of the probe is higher than a preset speed threshold, the window length is set to a first length. When the traveling speed of the probe is lower than or equal to the preset speed threshold, the window length is set to a second length greater than the first length.

[0108] The scaling factor is a positive real number. When the scaling factor is less than a preset lower threshold or greater than a preset upper threshold, the method triggers the acquisition of new digital geomagnetic time series data with the first quantization bit width to update the reference frame. The lower threshold is 50% of the estimated scaling factor value in the previous update cycle, and the upper threshold is 200% of the estimated scaling factor value in the previous update cycle.

[0109] The signal reconstructed after scaling the time series differential digital output by the scaling factor is obtained by multiplying the sign value of the differential digital output with the current estimated value of the scaling factor at each sampling time. When the differential digital output is at a logic high level, the reconstructed signal value is a positive scaling factor value, and when the differential digital output is at a logic low level, the reconstructed signal value is a negative scaling factor value.

[0110] Further, after estimating the equivalent magnetic moment data of the one or more weak magnetic targets, based on the estimated spatial distribution data of the weak magnetic targets, the contribution of the magnetic anomaly field theoretically generated by the weak magnetic targets is subtracted from the first digital geomagnetic time series data to separate the residual geomagnetic signal component mainly containing the contribution of large fixed ferromagnetic interference objects (the residual geomagnetic signal component may still contain geomagnetic background field components that have not been completely eliminated). In subsequent processing, the geomagnetic background field component is removed as a known interference term through least squares fitting, thereby estimating the magnetic moment tensor data of the large fixed ferromagnetic interference objects. The magnetic moment tensor data is then spatially registered and numerically superimposed with the equivalent magnetic moment data corresponding to the one or more weak magnetic targets to generate a complete magnetic source distribution map of all magnetic targets in the detection area.

[0111] The magnetic moment tensor data is a second-order tensor. The nine components of the second-order tensor correspond to the contribution coefficients of the magnetic moment components of the large fixed ferromagnetic interference object to the magnetic field components of the three orthogonal axes. Among the nine components of the second-order tensor, the diagonal components correspond to the magnetic field components generated by the magnetic moment components of the same axis in that axis, and the off-diagonal components correspond to the magnetic field components generated by the magnetic moment components of different axes in different axes. The second-order tensor is a symmetric tensor, that is, the magnetic field component generated by the magnetic moment component corresponding to axis i in axis j is equal to the magnetic field component generated by the magnetic moment component corresponding to axis j in axis i, where i and j represent any two different axes among the three orthogonal axes.

[0112] The contribution of the magnetic anomaly field generated by the weak magnetic target is obtained by calculation through a forward model. The forward model takes the spatial coordinates and equivalent magnetic moment data of the weak magnetic target as input, the spatial coordinates of the geomagnetic sensor array as the field point position, and outputs the theoretical magnetic anomaly field vector at the field point position. The calculation formula of the forward model is the magnetic dipole field formula. The magnetic dipole field formula includes the displacement vector from the field point position to the spatial coordinates of the weak magnetic target, the cube of the magnitude of the displacement vector, the dot product of the equivalent magnetic moment vector and the displacement vector, and the vector product term of the displacement vector and the equivalent magnetic moment vector.

[0113] The complete magnetic source distribution map is output in the form of a two-dimensional horizontal slice map. The two-dimensional horizontal slice map corresponds to the distribution of magnetic target bodies at a preset depth layer in the detection area. Each preset depth layer generates an independent two-dimensional horizontal slice map. The preset depth layers are set at equal depth intervals, which are half of the minimum vertical resolution of the target body in the detection area. The grayscale value or color value of each pixel in the two-dimensional horizontal slice map corresponds to the total magnetic moment amplitude at the horizontal position of the pixel and the preset depth layer.

[0114] The operation of subtracting the contribution of the magnetic anomaly field generated by the theory of the weak magnetic target from the first digital geomagnetic time series data is performed in the frequency domain. First, the first digital geomagnetic time series data is transformed to the frequency domain to obtain a first spectrum, and the theoretical magnetic anomaly field contribution is transformed to the frequency domain to obtain a second spectrum. The second spectrum is subtracted from the first spectrum point by point, and then the subtracted spectrum is restored to the time domain through inverse Fourier transform to obtain the residual geomagnetic signal component.

[0115] The spatial coordinate registration is achieved by aligning the spatial coordinate system of the magnetic moment tensor data with the spatial coordinate system of the equivalent magnetic moment data through rotation and translation transformation. The parameters of the rotation and translation transformation are determined by minimizing the weighted distance between the spatial centroid of the strong magnetic disturbance body in the magnetic moment tensor data and the spatial centroid of all weak magnetic target bodies in the equivalent magnetic moment data.

[0116] The complete magnetic source distribution map also includes a confidence ellipsoid for each magnetic target. The confidence ellipsoid is obtained by eigenvalue decomposition of the covariance matrix of the spatial position estimate of the magnetic target. The diagonal elements of the covariance matrix correspond to the estimated variance of each component of the spatial position, and the off-diagonal elements of the covariance matrix correspond to the estimated covariance between each component of the spatial position. The three principal axes of the confidence ellipsoid are determined by the three eigenvectors of the covariance matrix, and the half-lengths of the three principal axes of the confidence ellipsoid are determined by multiplying the arithmetic square root of the three eigenvalues ​​by a preset confidence scaling factor.

[0117] Further, estimating the magnetic moment tensor data of the large fixed ferromagnetic interference object includes: after the simulated cancellation signal has eliminated the first-order major contribution of the large fixed ferromagnetic interference object, using the residual higher-order components in the first digital geomagnetic time series data, and using the least squares fitting method to back-calculate the higher-order magnetic multipole moment parameters of the large fixed ferromagnetic interference object from the time series difference digital output of the subsequent time window.

[0118] The spatial decay rate of the quadrupole moment component in the higher-order magnetic multipole moment parameters is inversely proportional to the fourth power of the distance, and the spatial decay rate of the octrupole moment component is inversely proportional to the sixth power of the distance. The least squares fitting method utilizes the difference in spatial decay rate to separate the contribution of each order of magnetic multipole moment parameters from the measurements at different distances. By introducing the reciprocal of the fourth power and the reciprocal of the sixth power of the distance as coefficients for different orders of magnetic multipole moment parameters into the observation equation, the measurement values ​​at closer distances contribute to all orders of magnetic multipole moment parameters when solving the observation equation, while the measurement values ​​at farther distances mainly contribute to the lower-order magnetic multipole moment parameters. Thus, the separation of each order of magnetic multipole moment parameters is achieved by utilizing the sensitivity difference of the measurement values ​​at different distances.

[0119] The least squares fitting method employs a weighted least squares algorithm. In this algorithm, the weight coefficient of each sampling point is inversely proportional to the residual amplitude of the simulated cancellation signal at that sampling point. Sampling points with larger residual amplitudes are assigned smaller weights. The residual amplitude is obtained by power estimation of the differential-mode signal output by the simulated differential summing node. The weight coefficients are recalculated in each iteration based on the estimated value of the magnetic multipole moment parameter obtained in the current iteration, thus forming an adaptive iterative weighted relationship between the weight coefficients and the residual amplitude.

[0120] The higher-order magnetic multipole moment parameters obtained by reverse calculation are used to update the theoretical magnetic anomaly field spatial distribution function in the geomagnetic background field cancellation signal generation circuit, so that the simulated cancellation signal in subsequent sampling periods includes the contribution of the higher-order magnetic multipole moment parameters, thereby achieving the cancellation of higher-order magnetic field components of the large fixed ferromagnetic interference object. The update is achieved by substituting the estimated value of the higher-order magnetic multipole moment parameters into the magnetic multipole expansion to calculate the higher-order correction term of the magnetic field spatial distribution of the large fixed ferromagnetic interference object, and superimposing the higher-order correction term onto the original theoretical magnetic anomaly field spatial distribution function to form the updated theoretical magnetic anomaly field spatial distribution function.

[0121] The residual higher-order components are obtained by performing multi-scale wavelet decomposition on the first digital geomagnetic time series data. The multi-scale wavelet decomposition decomposes the first digital geomagnetic time series data into multiple detail components of different scales and an approximate component. The approximate component corresponds to the magnetic dipole term of the large fixed ferromagnetic interference object, and the detail component corresponds to the higher-order magnetic multipole moment parameter. The scale range of the retained detail components is controlled by selecting the number of decomposition layers. The selected number of decomposition layers makes the spatial frequency corresponding to the retained detail components higher than the spatial frequency of the magnetic dipole term.

[0122] Each column of the observation matrix in the linear observation equation set corresponds to the theoretical magnetic field response generated at the location of the geomagnetic sensor array by an unknown high-order magnetic multipole moment parameter. The theoretical magnetic field response is obtained by evaluating the field expressions of each order magnetic multipole moment source at the spatial coordinates of the geomagnetic sensor array. The field expressions of each order magnetic multipole moment source are obtained by performing a Taylor expansion of the magnetic dipole field formula and retaining it to the corresponding order.

[0123] When the root mean square value of the fitting residual of the least square fitting method exceeds a preset residual threshold, the method increases the order of the higher-order magnetic multipole moment parameter from containing only the magnetic quadrupole moment component to containing both the magnetic quadrupole moment component and the magnetic octupole moment component, and re-executes the least square fitting method. The preset residual threshold is five percent of the root mean square value of the first digital geomagnetic time series data.

[0124] Furthermore, before determining the digital cancellation waveform, the method includes determining the geomagnetic background grid model of the detection area and the spatial location list of the known strong magnetic interference bodies based on one or more of the following: historical geomagnetic diurnal variation observation data, current array spatial pose data provided by the inertial navigation system, or a geographic information system map of urban underground infrastructure corresponding to the detection area. The geographic information system map of urban underground infrastructure includes the design coordinates of existing large ferromagnetic culverts, manhole covers, reinforced concrete structures, and subway tunnels, as well as their theoretical magnetic anomaly feature fingerprint database.

[0125] The geomagnetic field background grid model of the detection area is obtained by Kriging spatial interpolation of the historical geomagnetic diurnal variation observation data. The Kriging spatial interpolation takes the latitude and longitude grid nodes of the detection area as the interpolation target points, and takes the observation point locations and observation values ​​of the historical geomagnetic diurnal variation observation data as inputs, and outputs the background magnetic field estimate at each grid node. The Kriging spatial interpolation first calculates the variogram between the observation point locations, which describes the spatial correlation of the observation values ​​between observation points at different distances. Then, it determines the weight coefficient of each observation point to each interpolation target point based on the variogram, and uses the weight coefficients to perform a weighted summation of the observation values ​​to obtain the background magnetic field estimate.

[0126] The theoretical magnetic anomaly feature fingerprint database consists of type labels, spatial size parameters, material magnetic susceptibility parameters, and standard magnetic anomaly spatial distribution patterns for each strong magnetic interference body. The standard magnetic anomaly spatial distribution patterns are pre-calculated through finite element numerical simulation. The finite element numerical simulation determines the corresponding geometric model based on the type labels of each strong magnetic interference body, determines the size of the geometric model based on the spatial size parameters, and determines the magnetic properties of the geometric model based on the material magnetic susceptibility parameters. A uniform background magnetic field boundary condition is applied to the space surrounding the geometric model, and Maxwell's equations are solved to obtain the magnetic anomaly field distribution at each point in the space surrounding the geometric model. The solution results are used as the standard magnetic anomaly spatial distribution patterns.

[0127] The geographic information system map of the urban underground infrastructure is stored in vector layer format. Each strong magnetic interference body in the vector layer is represented by a polygonal surface feature. The attribute table of the polygonal surface feature contains information on the burial depth, orientation, length, and material type of the strong magnetic interference body. The geometric contour of the polygonal surface feature corresponds to the projection range of the strong magnetic interference body on the ground surface. The burial depth is the vertical distance from the top of the strong magnetic interference body to the ground surface. The orientation is the azimuth angle of the major axis direction of the strong magnetic interference body relative to the geographic north direction. The material type is used to index the corresponding standard magnetic anomaly spatial distribution pattern in the theoretical magnetic anomaly feature fingerprint database.

[0128] The historical geomagnetic diurnal variation observation data includes geomagnetic field intensity data continuously recorded at fixed time intervals within at least one complete Earth rotation cycle in the detection area. The fixed time interval is less than half of the period corresponding to the highest frequency component of the geomagnetic diurnal variation, so as to satisfy the Nyquist time sampling theorem. The complete Earth rotation cycle is 24 hours, and the highest frequency component of the geomagnetic diurnal variation corresponds to the second harmonic of the Earth rotation frequency.

[0129] The current array spatial pose data provided by the inertial navigation system is obtained through strapdown inertial navigation calculation. The strapdown inertial navigation calculation takes the angular velocity data output by the three-axis gyroscope and the specific force data output by the three-axis accelerometer in the inertial navigation system as input. The attitude angle is obtained by integrating the angular velocity data. The velocity and position are obtained by integrating the specific force data after coordinate transformation and gravitational acceleration compensation. The attitude angle is used to transform the geomagnetic field background grid model from the geographic coordinate system to the sensor coordinate system.

[0130] When the historical geomagnetic diurnal variation observation data, the current array spatial pose data provided by the inertial navigation system, and the urban underground infrastructure geographic information system map are all available, the geomagnetic background grid model of the detection area is obtained by Bayesian fusion of the three. The Bayesian fusion uses the historical geomagnetic diurnal variation observation data as the observation likelihood function, the prior distribution provided by the urban underground infrastructure geographic information system map as the prior probability, and the current array spatial pose data provided by the inertial navigation system as the coordinate transformation constraint. The fused background magnetic field estimate is obtained by maximizing the posterior probability.

[0131] Further, determining the digital cancellation waveform includes: storing the first digital geomagnetic time series data in the local solid-state storage unit of the geomagnetic sensor array as a background field calibration baseline; at the geomagnetic interference waveform synthesis circuit of the geomagnetic sensor array, receiving the geomagnetic field background grid model of the detection area and the spatial location list of the known strong magnetic interference bodies from the local solid-state storage unit or one or more of the central microcontroller and embedded digital signal processing unit; according to the geomagnetic field background grid model of the detection area and the spatial location list of the known strong magnetic interference bodies, performing spatial coordinate transformation and magnetic field attenuation function weighting processing on the background field calibration baseline to reconstruct a theoretical interference waveform matching the current sensor array attitude in the digital domain, and generating the analog cancellation signal via digital-to-analog conversion, thereby eliminating the strong magnetic interference component corresponding to the one or more large fixed ferromagnetic interference bodies in the detection area in the analog domain, and simultaneously eliminating the time-varying magnetic interference component corresponding to the magnetic moving bodies with known motion trajectories in the detection area.

[0132] The spatial coordinate transformation includes converting the spatial position of the known strong magnetic interference body from the geographic coordinate system to a local Cartesian coordinate system with the centroid of the geomagnetic sensor array as the origin. The transformation is achieved through a combination operation of rotation matrix and translation vector. The rotation matrix is ​​determined by the attitude angles of the geomagnetic sensor array, including pitch angle, roll angle and heading angle. The translation vector is the negative of the coordinate value of the centroid of the geomagnetic sensor array in the geographic coordinate system. In the combination operation, the translation operation is performed first and then the rotation operation, or the rotation operation is performed first and then the translation operation.

[0133] The magnetic field attenuation function is the reciprocal of the cube of the distance between the centroid of the geomagnetic sensor array and the spatial location of the known strong magnetic interference body. The amplitude of each sampling point of the theoretical interference waveform is obtained by multiplying the amplitude of the corresponding sampling point in the background field calibration baseline by the reciprocal of the cube of the distance. When the distance between the centroid of the geomagnetic sensor array and the spatial location of the known strong magnetic interference body is less than a preset distance threshold, the magnetic field attenuation function is truncated to the reciprocal of the cube of the distance threshold to avoid the amplitude of the theoretical interference waveform being too large at close range due to singularity.

[0134] The time-varying magnetic interference component of a magnetic moving body with a known trajectory is obtained by discretizing the known trajectory of the magnetic moving body into time-series position points, calculating the corresponding theoretical magnetic anomaly field value for each time-series position point, and then splicing them in time order. The time interval between the time-series position points is equal to the sampling period of the analog-to-digital converter (ADC). The theoretical magnetic anomaly field value is calculated using the magnetic dipole field formula, and the magnetic dipole moment vector in the magnetic dipole field formula is determined by the product of the mass, material magnetic susceptibility, and volume of the magnetic moving body.

[0135] Before performing the spatial coordinate transformation, the background field calibration baseline is first subjected to high-pass filtering to remove the DC component and ultra-low frequency component in the background field calibration baseline. The cutoff frequency of the high-pass filtering is half of the lowest frequency component of the diurnal variation of the geomagnetic field in the detection area, so as to remove the influence of the long-term slow change of the geomagnetic field while retaining the frequency band component corresponding to the strong magnetic interference component.

[0136] Before generating the analog cancellation signal via the digital-to-analog converter, the theoretical interference waveform reconstructed in the digital domain undergoes digital domain delay compensation processing. This processing calculates the total delay time based on the sampling delay of the analog-to-digital converter (ADC), the setup time of the digital-to-analog converter (DAC), and the transmission delay of the analog differential summing node. The theoretical interference waveform is then shifted forward by this total delay time on the time axis, aligning the arrival time of the analog cancellation signal at the analog differential summing node with the arrival time of the original analog voltage signal at the analog differential summing node.

[0137] The weighted processing of the magnetic field attenuation function also includes the correction of the shape factor of the known strong magnetic interference body. The shape factor is determined according to the aspect ratio of the known strong magnetic interference body. When the aspect ratio of the known strong magnetic interference body is greater than a preset aspect ratio threshold, a shape factor correction term is introduced into the magnetic field attenuation function. The shape factor correction term is the square of the cosine of the angle between the long axis direction of the known strong magnetic interference body and the direction of the line connecting the geomagnetic sensor array to the known strong magnetic interference body.

[0138] Furthermore, the second quantization bit width is a one-bit zero-crossing quantization bit width, and the first quantization bit width is a sixteen-bit or higher quantization bit width; the quantization bit width of the analog-to-digital converter (ADC) is dynamically adjusted on the time axis so that the first quantization bit width is used to collect mixed data of background field and target field in the early stage of detection, and the second quantization bit width is used to collect differential polarity data after interference cancellation in the subsequent continuous scanning process.

[0139] The detection range threshold is determined based on the spatial distribution density of known strong magnetic interference bodies in the detection area. When the distribution density of the known strong magnetic interference bodies is higher than a preset density threshold, the detection range threshold is set to a smaller value so that the second quantization bit width is used before entering the high-density strong magnetic interference area, thereby reducing the data acquisition power consumption in the high-density interference area. When the distribution density of the known strong magnetic interference bodies is lower than or equal to the preset density threshold, the detection range threshold is set to a larger value so that the time of working with the first quantization bit width is extended in the low-density interference area, thereby acquiring more mixed background field and target field data as a reference frame for subsequent amplitude recovery.

[0140] The frequency of the recalibration sampling is positively correlated with the rate of change of the scaling factor. When the change of the scaling factor between adjacent update cycles increases, the frequency of the recalibration sampling increases accordingly. The recalibration sampling is achieved by inserting a calibration frame operating at the first quantization bit width between two consecutive sampling frames operating at the second quantization bit width. The number of sampling points in the calibration frame is equal to the number of sampling points in each sampling frame under the first quantization bit width. During the calibration frame, the geomagnetic sensor array suspends the differential superposition operation of the analog differential summing node, or maintains the differential superposition operation but samples the differential signal at the first quantization bit width.

[0141] During operation with the second quantization bit width, the sampling rate of the analog-to-digital converter (ADC) is higher than that during operation with the first quantization bit width, in order to compensate for the amplitude resolution loss caused by the reduction in quantization bit width through time oversampling. The oversampling rate of the time oversampling is the ratio of the sampling rate during operation with the first quantization bit width to the sampling rate during operation with the second quantization bit width. The oversampling rate is set to be equal to the square of the difference between the first quantization bit width and the second quantization bit width, such that the decibel reduction in quantization noise power due to oversampling cancels out the decibel resolution loss caused by the reduction in quantization bit width.

[0142] The operation of switching from the first quantization bit width to the second quantization bit width is triggered by a mode switching instruction issued by the central microcontroller and the embedded digital signal processing unit. The mode switching instruction is automatically generated by the central microcontroller and the embedded digital signal processing unit and sent to the quantization bit width control port of the analog-to-digital converter (ADC) when the cumulative detection mileage reaches the detection mileage threshold or the cumulative detection time reaches the detection time threshold. The mode switching instruction includes a target quantization bit width parameter and a switching effective time. The switching effective time is the rising edge of the next sampling clock after the mode switching instruction is issued.

[0143] The detection time threshold is determined based on the time-varying rate of the diurnal variation of the geomagnetic field in the detection area. When the time-varying rate of the diurnal variation of the geomagnetic field is higher than a preset rate threshold, the detection time threshold is set to be less than half of the time-varying period of the diurnal variation of the geomagnetic field, so as to ensure that the change in the geomagnetic background field within the detection time threshold does not exceed the preset background field change tolerance.

[0144] During the operation of the analog-to-digital converter (ADC) with the second quantization bit width, the method also monitors the peak amplitude of the differential analog signal. When the peak amplitude continuously exceeds the full-scale input range of the ADC with the second quantization bit width for multiple consecutive sampling periods, the method temporarily switches the quantization bit width of the ADC back to the first quantization bit width and increases the amplitude of the analog cancellation signal. After the peak amplitude of the differential analog signal returns to the full-scale input range, the method switches the quantization bit width of the ADC back to the second quantization bit width.

Claims

1. A multi-array geomagnetic underground pipeline detection system, characterized in that, include: One or more geomagnetic sensor arrays are fixed on a detection carrier in a three-dimensional spatial arrangement and configured to continuously acquire, passively, simulated voltage signals of magnetic anomaly fields formed by abnormal changes in the Earth's magnetic field caused by differences in magnetic susceptibility of one or more underground targets in the detection area of ​​the detection system. Each geomagnetic sensor array includes: The three-axis fluxgate sensor group is used to synchronously sense the changes in magnetic field strength of the geomagnetic field in three orthogonal axes and output the simulated voltage signal of the magnetic anomaly field to the subsequent conditioning circuit. An analog-to-digital converter (ADC) is configured to perform time-series quantization on the analog voltage signal at a preset sampling rate. For a first time-series segment in the analog voltage signal, the ADC is configured to quantize the first analog voltage signal segment using a first quantization bit width to generate first digital geomagnetic time-series data from the initial sampling period of the first analog voltage signal segment. A geomagnetic background field cancellation signal generation circuit is configured to generate a digital cancellation waveform based on pre-stored geomagnetic field spatial distribution background field data or historical detection baseline data of the detection area. The digital cancellation waveform is then converted from digital to analog to form an analog cancellation signal. The amplitude and phase of the analog cancellation signal are set to cancel the strong magnetic interference components generated by one or more large fixed ferromagnetic objects or known dynamic magnetic moving bodies at known spatial locations in the detection area from the magnetic anomaly field analog voltage signal before analog-to-digital conversion. An analog differential summing node is provided between the analog signal output terminal of the triaxial fluxgate sensor group and the analog signal input terminal of the analog-to-digital converter (ADC). The analog differential summing node is used to differentially superimpose the analog voltage signal and the analog cancellation signal before the analog voltage signal enters the ADC for quantization. For subsequent time-series segments in the analog voltage signal, the geomagnetic sensor array is configured as follows: After the differential superposition is completed at the analog differential summation node, and when the cumulative detection mileage reaches a preset detection mileage threshold or the cumulative detection time reaches a preset detection time threshold, the first quantization bit width of the analog-to-digital converter (ADC) is dynamically adjusted to a second quantization bit width that is smaller than the first quantization bit width; the detection mileage threshold is determined based on the spatial distribution density of known strong magnetic interference bodies in the detection area, and the detection time threshold is determined based on the time-varying rate of the diurnal variation of the geomagnetic field in the detection area; Using the analog-to-digital converter (ADC) measurement output value at the second quantization bit width, and combining it with the amplitude reference in the first digital geomagnetic time series data obtained during the initial sampling period, the equivalent magnetic moment vector data of one or more weak magnetic targets within the detection area is determined by calculating the scaling factor corresponding to the least square error between the first digital geomagnetic time series data and the magnetic anomaly amplitude reconstructed by the ADC measurement output value.

2. The multi-array geomagnetic underground pipeline detection system according to claim 1, characterized in that, Each of the geomagnetic sensor arrays is configured to change the system power consumption state by dynamically switching the quantization bit width of the analog-to-digital converter (ADC) from the first quantization bit width to the second quantization bit width; the second quantization bit width is a one-bit zero-crossing quantization bit width, that is, only the voltage polarity sign of the differentially superimposed analog signal is output; each of the geomagnetic sensor arrays is also configured to use the signal amplitude envelope in the first digital geomagnetic time series data obtained during the initial sampling period to nonlinearly scale the magnetic anomaly amplitude corresponding to the polarity sequence output under the one-bit zero-crossing quantization bit width, the nonlinear scaling being achieved by constructing an amplitude mapping relationship between the polarity sequence and the first digital geomagnetic time series data.

3. The multi-array geomagnetic underground pipeline detection system according to claim 1, characterized in that, The geomagnetic background field cancellation signal generation circuit and the analog differential summation node are used for: After the differential superposition of the strong magnetic interference components is completed, the residual magnetic anomaly components caused by the one or more weak magnetic targets are eliminated from the sparse signal reconstruction algorithm from the analog-to-digital converter (ADC) measurement output value in order to re-estimate the magnetic moment tensor data of the one or more large fixed ferromagnetic interference objects or known dynamic magnetic moving objects. The equivalent magnetic moment vector data and the magnetic moment tensor data of the weak magnetic target are spatially superimposed to generate a comprehensive magnetic distribution dataset of the target in the detection area. The comprehensive magnetic distribution dataset contains the classification identifiers and spatial separation information of the weak magnetic target and the large fixed ferromagnetic interference object.

4. The multi-array geomagnetic underground pipeline detection system according to claim 1, characterized in that, Also includes: A central microcontroller and an embedded digital signal processing unit are coupled via a data bus to the digital output port of the analog-to-digital converter (ADC) of each of the one or more geomagnetic sensor arrays to receive one or more of the first digital geomagnetic time series data or the equivalent magnetic moment vector data. And each of the one or more geomagnetic sensor arrays further includes: A local solid-state storage unit is configured to cache the first digital geomagnetic time series data from the initial sampling period as reference data for subsequent amplitude recovery and background field update; as well as Geomagnetic interference waveform synthesis circuit, wherein the geomagnetic interference waveform synthesis circuit is configured to: The system receives pre-calculated background field data of the geomagnetic field spatial distribution of the detection area from one or more of the local solid-state storage unit or the central microcontroller and embedded digital signal processing unit. The data includes the spatial coordinates of known large ferromagnetic structures and the spatial distribution function of the theoretical magnetic anomaly field they generate. Based on the theoretical magnetic anomaly field spatial distribution function, the cached first digital geomagnetic time series data is subjected to time-domain waveform reconstruction and amplitude mapping processing to generate the digital cancellation waveform. The digital cancellation waveform is then converted from digital to analog to form the analog cancellation signal.

5. The multi-array geomagnetic underground pipeline detection system according to claim 4, characterized in that, The geomagnetic interference waveform synthesis circuit includes: A digital waveform generation and analog-to-digital conversion driving module is configured to generate multi-channel digital cancellation waveform data in real time based on the theoretical magnetic anomaly field spatial distribution function and the current spatial pose information of the sensor array. A digital-to-analog converter (DAC) includes a digital input port coupled to the digital waveform generation and digital-to-analog conversion drive module and an analog output port configured to convert the multi-channel digital cancellation waveform data into a continuous analog cancellation voltage signal. The analog differential summing node includes a first differential input port for receiving the original analog voltage signal output by the triaxial fluxgate sensor group, a second differential input port coupled to the analog output port of the digital-to-analog converter (DAC), and a differential signal output port for outputting the difference between the original analog voltage signal and the continuous analog cancelling voltage signal. Furthermore, the analog-to-digital converter (ADC) includes an analog sampling input coupled to the differential signal output port of the analog differential summing node and a digital data output coupled to the central microcontroller and the embedded digital signal processing unit.

6. A method for identifying underground pipelines using multi-array geomagnetic arrays, characterized in that, The method includes: At the triaxial fluxgate sensor of each geomagnetic sensor array in the multi-array geomagnetic underground pipeline detection system, the simulated voltage signal of the magnetic anomaly field generated by one or more underground targets in the detection area of ​​the geomagnetic sensor array due to the difference in magnetic susceptibility and superimposed on the Earth's background magnetic field is continuously received in a passive manner. The analog-to-digital converter (ADC) of the geomagnetic sensor array is used to quantize the simulated voltage signal of the magnetic anomaly field within the first time window with a first quantization bit width to generate first digital geomagnetic time series data, and the data is stored in a local reference data memory. The first digital geomagnetic time series data stored in the local reference data memory is subjected to digital domain cancellation processing. The magnetic field component contributed by the known strong magnetic interference body is subtracted from the first digital geomagnetic time series data using the digital cancellation waveform generated by the geomagnetic background field cancellation signal generation circuit, to obtain a pure weak magnetic target body plus background field reference signal. The pure reference signal replaces the original first digital geomagnetic time series data as the observation value for solving the amplitude scaling factor in the subsequent steps. Using the geomagnetic background field cancellation signal generation circuit of the geomagnetic sensor array, based on the pre-loaded geomagnetic background grid model of the detection area and the list of known strong magnetic interference bodies, combined with the low-frequency trend term in the first digital geomagnetic time series data, a digital cancellation waveform for canceling strong magnetic interference components is determined, and an analog cancellation signal is generated via digital-to-analog conversion. Based on the simulated cancellation signal, the simulated differential summation node of the geomagnetic sensor array is used to remove the strong magnetic interference voltage component caused by the one or more large fixed ferromagnetic interference objects or known dynamic magnetic moving bodies from the real-time magnetic anomaly field simulated voltage signal in the subsequent time window to obtain the differential mode simulated signal. After the removal is completed at the analog differential summation node, the quantization bit width of the analog-to-digital converter (ADC) is dynamically adjusted from the first quantization bit width to a second quantization bit width that is smaller than the first quantization bit width. The differential analog signal is quantized using the analog-to-digital converter (ADC) with the second quantization bit width to generate a time-series differential digital output. Combined with the clean reference signal in the local reference data memory, an objective optimization function is constructed with the clean reference signal as the observation and the symbol sequence of the time-series differential digital output as the prior constraint. The amplitude scaling coefficient in the objective optimization function is solved to determine the equivalent magnetic moment data and spatial location data corresponding to one or more weak magnetic targets in the detection area.

7. The method according to claim 6, characterized in that, Also includes: Calculate the signal intensity scaling factor between the first digital geomagnetic time series data and the inverted magnetic moment data corresponding to the one or more weak magnetic targets. The scaling factor is obtained by solving the least squares error between the first digital geomagnetic time series data and the signal reconstructed by scaling the time series differential digital output using the scaling factor.

8. The method according to claim 6, characterized in that, Also includes: After estimating the equivalent magnetic moment data of the one or more weak magnetic targets, based on the estimated spatial distribution data of the weak magnetic targets, the contribution of the magnetic anomaly field theoretically generated by the weak magnetic targets is subtracted from the first digital geomagnetic time series data to separate the residual geomagnetic signal component that mainly includes the contribution of large fixed ferromagnetic interference objects (the residual geomagnetic signal component may still contain geomagnetic background field components that have not been completely eliminated). In subsequent processing, the geomagnetic background field components are removed as known interference terms through least squares fitting, thereby estimating the magnetic moment tensor data of the large fixed ferromagnetic interference objects. The magnetic moment tensor data is spatially registered and numerically superimposed with the equivalent magnetic moment data corresponding to the one or more weakly magnetic targets to generate a complete magnetic source distribution map of all magnetic targets in the detection area.

9. The method according to claim 8, characterized in that, Estimating the magnetic moment tensor data of the large fixed ferromagnetic interference object includes: after the simulated cancellation signal has eliminated the first-order major contribution of the large fixed ferromagnetic interference object, using the residual higher-order components in the first digital geomagnetic time series data, and using the least squares fitting method to back-calculate the higher-order magnetic multipole moment parameters of the large fixed ferromagnetic interference object from the time series difference digital output of the subsequent time window.

10. The method according to claim 6, characterized in that, Before determining the digital cancellation waveform, the method includes determining a geomagnetic background grid model of the detection area and a list of spatial locations of the known strong magnetic interference bodies based on one or more of the following: historical geomagnetic diurnal variation observation data, current array spatial pose data provided by an inertial navigation system, or a geographic information system map of urban underground infrastructure corresponding to the detection area. The geographic information system map of urban underground infrastructure includes the design coordinates of existing large ferromagnetic culverts, manhole covers, reinforced concrete structures, and subway tunnels, as well as their theoretical magnetic anomaly feature fingerprint database.

11. The method according to claim 6, characterized in that, Determining the digital cancellation waveform includes: The first digital geomagnetic time series data is stored in the local solid-state storage unit of the geomagnetic sensor array as a background field calibration baseline; At the geomagnetic interference waveform synthesis circuit of the geomagnetic sensor array, the geomagnetic field background grid model of the detection area and the list of spatial locations of the known strong magnetic interference bodies are received from one or more of the local solid-state storage unit or the central microcontroller and the embedded digital signal processing unit. Based on the geomagnetic field background grid model of the detection area and the list of spatial locations of the known strong magnetic interference objects, the background field calibration baseline is subjected to spatial coordinate transformation and magnetic field attenuation function weighting to reconstruct a theoretical interference waveform that matches the current sensor array attitude in the digital domain. The analog cancellation signal is then generated via digital-to-analog conversion to eliminate the strong magnetic interference components corresponding to one or more large fixed ferromagnetic interference objects in the detection area in the analog domain, and simultaneously eliminate the time-varying magnetic interference components corresponding to magnetic moving objects with known motion trajectories in the detection area.

12. The method according to claim 6, characterized in that, The second quantization bit width is a one-bit zero-crossing quantization bit width, and the first quantization bit width is a sixteen-bit or higher quantization bit width; the quantization bit width of the analog-to-digital converter (ADC) is dynamically adjusted on the time axis so that the first quantization bit width is used to collect mixed data of background field and target field in the early stage of detection, and the second quantization bit width is used to collect differential polarity data after interference cancellation in the subsequent continuous scanning process.

Images