Electromagnetic wave CT underground cavity identification method based on phase information

Through the joint inversion model constrained by phase gradient terms and the dual-frequency differential phase method, combined with the morphological enhancement algorithm, the problems of high missed detection rate and poor stability of electromagnetic wave CT in identifying small-scale underground voids are solved, and high-precision void identification and quantification are achieved.

CN120802365APending Publication Date: 2025-10-17INST OF GEOPHYSICAL & GEOCHEMICAL EXPLORATION CHINESE ACAD OF GEOLOGICAL SCI
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202511030392.0
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2025-07-25
Publication Date
2025-10-17

AI Technical Summary

Technical Problem

Existing electromagnetic wave CT methods have problems such as high missed detection rate, insufficient resolution and poor stability when identifying small-scale underground cavities. Traditional phase inversion models are easily affected by noise in complex environments, making it difficult to accurately identify and quantify cavities.

Method used

An electromagnetic wave CT method based on phase information is adopted. Through a joint inversion model with the phase gradient term as the core constraint, combined with the dual-frequency differential phase method and morphological enhancement algorithm, phase ambiguity and multipath effect are eliminated, the weight is dynamically adjusted, and the dielectric constant distribution is accurately extracted.

Benefits of technology

The ability to identify small-scale voids is significantly improved, the robustness and adaptability of the method are enhanced, and it can accurately identify the location and volume of voids in complex noise environments and provide confidence probability.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN120802365A_ABST
    Figure CN120802365A_ABST
Patent Text Reader

Abstract

The invention discloses an electromagnetic wave CT underground cavity identification method based on phase information, and the method comprises the steps: A1, data collection: arranging a transmitting antenna and receiving antenna array in an underground space, transmitting an electromagnetic wave signal, and synchronously recording the amplitude attenuation and phase offset of the electromagnetic wave of each receiving point; a2, phase decoupling processing: performing full-period phase unwrapping operation on the original phase data, eliminating 2pi period fuzziness of the phase data, and calculating a normalized phase gradient of each propagation path; and A3, joint inversion modeling: constructing a target optimization function taking a phase gradient term as a core constraint and an amplitude attenuation term as a secondary constraint. The invention relates to the technical field of geophysical exploration. According to the electromagnetic wave CT underground cavity identification method based on the phase information, a joint inversion model is constructed by taking a phase gradient term as a core constraint, and the characteristic that the phase information is highly sensitive to tiny change of a dielectric constant of a medium is fully utilized.
Need to check novelty before this filing date? Find Prior Art

Description

TECHNICAL FIELD

[0001] The application relates to the technical field of geophysical prospecting, in particular to an electromagnetic wave CT underground cavity identification method based on phase information. BACKGROUND

[0002] Underground cavity detection is of great significance to tunnel engineering, mine safety and urban underground space development. Traditional methods mainly rely on seismic wave CT or electromagnetic wave amplitude attenuation CT technology. Although seismic wave CT has a large penetration depth, its resolution is limited (usually > 1 meter), and it is difficult to identify small-scale cavities. Existing electromagnetic wave CT methods generally rely on amplitude attenuation inversion, but when electromagnetic waves propagate in a non-uniform medium, the amplitude signal is easily affected by environmental noise and scattering loss, and is not sensitive enough to small changes in dielectric constant (such as the sudden drop in dielectric constant caused by cavities), resulting in a high detection rate of more than 40% for small cavities (diameter < 0.5 meters) or weak anomaly areas. In addition, although phase information contains higher precision medium characteristics, it is limited by two technical bottlenecks: one is the 2π cycle ambiguity of the original phase, and the unwinding error will accumulate with the propagation distance; the other is the random jitter of phase caused by multipath effect, especially in the underground complex structure environment, the stability of traditional single-frequency phase inversion is poor. Existing joint inversion models mostly use fixed weight distribution (such as phase / amplitude weight 6:4), which cannot adapt to noise fluctuations and is prone to false anomalies in low signal-to-noise ratio areas. The extraction of cavity boundaries generally relies on simple threshold segmentation, ignoring the geometric continuity, resulting in a cavity volume estimation error of more than 20%. Therefore, it is urgent to develop a high-sensitivity, anti-interference phase-driven electromagnetic wave CT method to realize the accurate identification and quantitative analysis of underground cavities. SUMMARY

[0003] To overcome the above-mentioned defects, the application is implemented by the following technical scheme: an electromagnetic wave CT underground cavity identification method based on phase information, comprising the following steps:

[0004] Step A1, data acquisition:

[0005] The transmitting antenna and receiving antenna array are arranged in the underground space, electromagnetic wave signals are transmitted, and the amplitude attenuation and phase shift of the electromagnetic wave at each receiving point are recorded synchronously;

[0006] Step A2, phase decoupling processing:

[0007] The original phase data is subjected to full-cycle phase unwinding operation to eliminate the 2π cycle ambiguity of the phase data, and the normalized phase gradient of each propagation path is calculated;

[0008] Step A3, joint inversion modeling:

[0009] An objective optimization function is constructed with a phase gradient term as a core constraint and an amplitude attenuation term as a secondary constraint, and a specific calculation formula of the objective optimization function is:

[0010] Objective optimization function = minimum phase gradient to dielectric constant mapping error weighted sum + amplitude measured value and calculation value error weighted sum;

[0011] Wherein, the phase gradient term weight coefficient a is greater than the amplitude attenuation term weight coefficient b;

[0012] By constructing a joint inversion model with the phase gradient term as the core constraint (weight a> b), the characteristics of the phase information being highly sensitive to small changes in the dielectric constant of the medium are fully utilized. Compared with the traditional method mainly relying on amplitude attenuation, this method can detect the area of abnormal reduction of dielectric constant (cavity) earlier and more accurately, and significantly improves the identification ability of small-scale or weak abnormal cavities.

[0013] Step A4, cavity identification:

[0014] The dielectric constant spatial distribution is solved by using a regularization iterative algorithm, and the abnormal area with dielectric constant less than or equal to 3.0 is extracted to determine the underground cavity.

[0015] Preferably, the frequency band of the electromagnetic wave signal in step A1 is wide frequency of 10MHz to 1GHz.

[0016] Preferably, the specific calculation formula of the normalized phase gradient in step A2 is:

[0017] Normalized phase gradient = (observed phase - reference phase) / (angular frequency x propagation path length).

[0018] Preferably, in the phase decoupling processing of step A2, a dual-frequency differential phase method is used to eliminate system error, and the specific operation steps of the dual-frequency differential phase method are:

[0019] Step S1, synchronously transmitting low frequency f1 and high frequency f2 signals, wherein the low frequency f1 is greater than or equal to 5 times the high frequency f2;

[0020] Step S2, calculating the differential phase:

[0021] Delta Phi_diff = Phi(f1) - k*Phi(f2);

[0022] Wherein, k is a frequency ratio coefficient. The phase jitter caused by multipath effect is suppressed by Delta Phi_diff.

[0023] The phase data is processed by using a dual-frequency differential phase method, and through combination operation of high-frequency signals and low-frequency signals, phase jitter and instability caused by multipath propagation due to complex underground environment can be effectively offset or significantly weakened, and reliability of the phase data and stability of the inversion result are improved.

[0024] Preferably, in the joint inversion model in step A3, a phase-amplitude confidence weight factor is introduced, and the specific operation steps are as follows:

[0025] In step D1, a confidence coefficient is defined, and the specific calculation formula of the confidence coefficient is as follows:

[0026] γ=e^(–σ 2 / ΔΦ_norm 2 );

[0027] Wherein, σ is the environmental noise variance.

[0028] In step D2, a dynamic adjustment objective function is set, and the specific calculation formula of the dynamic adjustment objective function is as follows:

[0029] α=γ·α0,β=(1–γ)β0;

[0030] In step D3, the weight is dynamically distributed.

[0031] The initial weight α0 of the phase term and the initial weight β0 of the amplitude term are set, and α0 is greater than β0.

[0032] The phase-amplitude confidence weight factor is introduced, and the weight of the phase term and the amplitude term in the objective function is dynamically adjusted according to the real-time calculated confidence γ. This makes the inversion process adapt to the environmental noise level and the quality of the current phase data. When the signal-to-noise ratio is high, more accurate phase information is relied on. When the signal-to-noise ratio is low, the weight of the relatively robust amplitude information is appropriately increased, which significantly improves the robustness and adaptability of the method in a complex noise environment.

[0033] Preferably, in step A4, the dielectric constant spatial distribution is solved by using a morphological enhancement algorithm, and the specific steps include the following steps:

[0034] In step F1, anisotropic diffusion filtering is performed on the dielectric constant spatial distribution to retain the edge of the cavity.

[0035] In step F2, a curvature-driven level set method is used to extract the three-dimensional boundary of the cavity.

[0036] In step F3, the position, volume and confidence probability of the cavity are output.

[0037] Preferably, in step F2, the three-dimensional boundary of the cavity is extracted by using the following specific operation steps:

[0038] In step F21, curvature calculation is performed.

[0039] Calculate the average curvature of the current boundary surface, reflecting the concave and convex degree of the surface;

[0040] Gradient operation is performed on the boundary function in three-dimensional space to obtain the directional change rate of each point;

[0041] Based on the gradient result, further divergence calculation is performed to obtain the accurate curvature value, and the curvature value of the convex region is positive and the curvature value of the concave region is negative;

[0042] Step F21, evolution rule execution;

[0043] According to the curvature value, the boundary moving direction is dynamically adjusted:

[0044] The positive curvature area is convex: the boundary shrinks inward;

[0045] The negative curvature area is concave: the boundary expands outward;

[0046] The zero curvature area is flat: the boundary remains stationary.

[0047] After solving the dielectric constant distribution, an advanced morphological enhancement algorithm is used for cavity extraction. This algorithm can effectively filter out noise interference, retain key edge information, and intelligently evolve according to the curvature of the boundary surface to accurately capture the three-dimensional boundary profile of complex shape cavities. The final output not only includes the cavity position, but also includes its accurate volume and confidence probability, providing more abundant and reliable information for engineering decision-making. Combined with wideband detection, it can adapt to the identification needs of cavities of different scales.

[0048] The application provides an electromagnetic wave CT underground cavity identification method based on phase information.

[0049] (1) The electromagnetic wave CT underground cavity identification method based on phase information constructs a joint inversion model by taking the phase gradient term as the core constraint, and fully utilizes the characteristics that phase information is highly sensitive to small changes in dielectric constant. Compared with the traditional method mainly relying on amplitude attenuation, this method can detect areas with abnormal reduction of dielectric constant earlier and more accurately, and significantly improves the identification ability of small-scale or weak abnormal cavities.

[0050] (2) The electromagnetic wave CT underground cavity identification method based on phase information processes phase data by using a dual-frequency differential phase method, and through the combination operation of high-frequency and low-frequency signals, it can effectively offset or significantly weaken the phase jitter and instability caused by multipath propagation caused by complex underground environment, thereby improving the reliability of the phase data and the stability of the inversion result.

[0051] (Three), the electromagnetic wave CT underground cavity recognition method based on phase information, the phase-amplitude confidence weight factor is introduced, and the weight of the phase item and the amplitude item in the objective function is dynamically adjusted according to the real-time calculated confidence γ. This makes the inversion process adapt to the noise level of the environment and the quality of the current phase data, and more depends on the accurate phase information when the signal-to-noise ratio is high, and appropriately increases the weight of the relatively robust amplitude information when the signal-to-noise ratio is low, which significantly improves the robustness and adaptability of the method in complex noise environment.

[0052] (Four), the electromagnetic wave CT underground cavity recognition method based on phase information, after solving the dielectric constant distribution, an advanced morphological enhancement algorithm is used for cavity extraction. The algorithm can effectively filter out noise interference, retain key edge information, and intelligently evolve according to the boundary surface curvature to accurately capture the three-dimensional boundary profile of complex shape cavities. The final output is not only the position of the cavity, but also its accurate volume and confidence probability, providing more abundant and reliable information for engineering decision-making. Combined with wideband detection, it can adapt to the recognition needs of cavities of different scales. BRIEF DESCRIPTION OF DRAWINGS

[0053] Figure 1 It is a structural schematic diagram of the whole application. DETAILED DESCRIPTION

[0054] The technical solutions in the embodiments of the application will be described clearly and completely below with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are only part of the embodiments of the application, not all. Based on the embodiments in the application, all other embodiments obtained by those skilled in the art without creative labor are within the scope of protection of the application.

[0055] Embodiment one, please refer to Figure 1 The application provides a technical solution:

[0056] An electromagnetic wave CT underground cavity recognition method based on phase information, comprising the following steps:

[0057] Step A1, data acquisition:

[0058] The transmitting antenna and receiving antenna array are arranged in the underground space, the electromagnetic wave signal is transmitted, and the amplitude attenuation and phase shift of the electromagnetic wave at each receiving point are recorded synchronously;

[0059] The frequency band of the electromagnetic wave signal is 1GHz wideband.

[0060] Step A2, phase decoupling processing:

[0061] The original phase data is subjected to full-cycle phase unwrapping operation to eliminate the 2π cycle ambiguity of the phase data, and the normalized phase gradient of each propagation path is calculated.

[0062] The specific calculation formula of the normalized phase gradient is:

[0063] The normalized phase gradient = (observed phase - reference phase) / (angular frequency x propagation path length).

[0064] In the phase decoupling processing, a dual-frequency differential phase method is used to eliminate system errors, and the specific operation steps of the dual-frequency differential phase method are:

[0065] Step S1, synchronously transmitting low-frequency f1 and high-frequency f2 signals, wherein the low-frequency f1 is equal to 7.5 times the high-frequency f2;

[0066] Step S2, calculating the differential phase:

[0067] ΔΦ_diff = Φ(f1) - k·Φ(f2);

[0068] Wherein, k is the frequency ratio coefficient. The phase jitter caused by the multipath effect is suppressed through ΔΦ_diff.

[0069] Step A3, joint inversion modeling:

[0070] An objective optimization function is constructed, which takes the phase gradient term as the core constraint and the amplitude attenuation term as the secondary constraint, and the specific calculation formula of the objective optimization function is:

[0071] The objective optimization function = minimum mapping error weighted sum of phase gradient to dielectric constant + error weighted sum of amplitude measured value and calculated value;

[0072] A phase-amplitude confidence weight factor is introduced, and the specific operation steps are:

[0073] Step D1, defining the confidence coefficient, and the specific calculation formula of the confidence coefficient is:

[0074] γ = e^(–σ 2 / ΔΦ_norm 2 );

[0075] Wherein, σ is the environmental noise variance;

[0076] Step D2, setting a dynamic adjustment objective function, and the specific calculation formula of the dynamic adjustment objective function is:

[0077] α = γ·α0, β = (1-γ)β0;

[0078] Step D3, dynamic weight distribution;

[0079] Set the initial weight α0 phase term and β0 amplitude term;

[0080] Wherein, the phase gradient term weight coefficient a is greater than the amplitude attenuation term weight coefficient b;

[0081] Step A4, cavity recognition:

[0082] The dielectric constant spatial distribution is solved by using a regularization iterative algorithm, and an abnormal area with a dielectric constant less than or equal to 3.0 is extracted to determine the underground cavity.

[0083] The dielectric constant spatial distribution is solved by using a morphological enhancement algorithm, specifically including the following steps:

[0084] Step F1, anisotropic diffusion filtering is performed on the dielectric constant spatial distribution to retain the cavity edge;

[0085] Step F2, a three-dimensional boundary of the cavity is extracted based on a curvature-driven level set method;

[0086] The specific operation steps of extracting the three-dimensional boundary of the cavity are:

[0087] Step F21, curvature calculation;

[0088] The average curvature of the current boundary surface is calculated to reflect the concave-convex degree of the surface;

[0089] Gradient operation is performed on the boundary function in three-dimensional space to obtain the directional change rate of each point;

[0090] Based on the gradient result, the divergence is further calculated to obtain the accurate curvature value, and the curvature value in the convex region is positive and the curvature value in the concave region is negative;

[0091] Step F21, evolution rule execution;

[0092] The boundary moving direction is dynamically adjusted according to the curvature value:

[0093] The convex region of the positive curvature protrudes: the boundary shrinks inward;

[0094] The concave region of the negative curvature is recessed: the boundary expands outward;

[0095] The flat region of the zero curvature is flat: the boundary remains stationary;

[0096] Step F3, output the cavity position, volume and confidence probability.

[0097] Example two, please refer to Figure 1 The application provides a technical solution: an electromagnetic wave CT underground cavity recognition method based on phase information, including the following steps:

[0098] Step A1, data acquisition:

[0099] The transmitting antenna and receiving antenna array are arranged in the underground space, the electromagnetic wave signal is transmitted, and the amplitude attenuation and phase shift of the electromagnetic wave at each receiving point are recorded synchronously;

[0100] The frequency band of the electromagnetic wave signal is wideband of 10MHz.

[0101] Step A2, phase decoupling processing:

[0102] Full-cycle phase unwrapping is performed on the original phase data to eliminate the 2π cycle ambiguity of the phase data, and the normalized phase gradient of each propagation path is calculated.

[0103] The specific calculation formula of the normalized phase gradient is:

[0104] Normalized phase gradient=(observed phase-reference phase) / (angular frequency*propagation path length).

[0105] In the phase decoupling processing, a dual-frequency differential phase method is used to eliminate system errors, and the specific operation steps of the dual-frequency differential phase method are:

[0106] Step S1, synchronously transmitting low-frequency f1 and high-frequency f2 signals, wherein the low-frequency f1 is equal to 5 times the high-frequency f2;

[0107] Step S2, calculating the differential phase:

[0108] ΔΦ_diff=Φ(f1)-k·Φ(f2);

[0109] Wherein, k is the frequency ratio coefficient. Through ΔΦ_diff, the phase jitter caused by multipath effect is suppressed.

[0110] Step A3, joint inversion modeling:

[0111] An objective optimization function is constructed with the phase gradient term as the core constraint and the amplitude attenuation term as the secondary constraint, and the specific calculation formula of the objective optimization function is:

[0112] Objective optimization function=minimum mapping error weighted sum of phase gradient to dielectric constant+error weighted sum of amplitude measured value and calculated value;

[0113] A phase-amplitude confidence weight factor is introduced, and the specific operation steps are:

[0114] Step D1, defining the confidence coefficient, and the specific calculation formula of the confidence coefficient is:

[0115] γ=e^(–σ 2 / ΔΦ_norm 2 );

[0116] Wherein, σ is the environmental noise variance;

[0117] Step D2, setting a dynamic adjustment objective function, and the specific calculation formula of the dynamic adjustment objective function is:

[0118] a = y a0, b = (1 - y) b0;

[0119] Step D3, dynamic weight distribution;

[0120] Set initial weight a0phase term, b0amplitude term;

[0121] Wherein, the phase gradient term weight coefficient a is greater than the amplitude attenuation term weight coefficient b;

[0122] Step A4, cavity recognition:

[0123] The dielectric constant space distribution is solved by using a regularization iterative algorithm, and the abnormal area with dielectric constant less than or equal to 3.0 is extracted to determine the underground cavity.

[0124] The dielectric constant space distribution is solved by using a morphological enhancement algorithm, which includes the following steps:

[0125] Step F1, anisotropic diffusion filtering is performed on the dielectric constant space distribution to retain the cavity edge;

[0126] Step F2, extract the three-dimensional boundary of the cavity based on the curvature driven level set method;

[0127] The specific operation steps of extracting the three-dimensional boundary of the cavity are:

[0128] Step F21, curvature calculation;

[0129] Calculate the average curvature of the current boundary surface to reflect the concave-convex degree of the surface;

[0130] Gradient operation is performed on the boundary function in three-dimensional space to obtain the directional change rate of each point;

[0131] Based on the gradient result, further calculate the divergence to obtain the accurate curvature value, the curvature value in the convex region is positive, and the curvature value in the concave region is negative;

[0132] Step F21, evolution rule execution;

[0133] According to the curvature value, dynamically adjust the moving direction of the boundary:

[0134] Positive curvature area protrusion: the boundary shrinks inward;

[0135] Negative curvature area depression: the boundary expands outward;

[0136] Zero curvature area flat: the boundary remains stationary;

[0137] Step F3, output the cavity position, volume and confidence probability.

[0138] In use, first, the transmitting antenna array and the receiving antenna array are laid out in the underground detection area, the array spacing is set to 1-10 meters according to the target depth; the dual-frequency electromagnetic wave signals are synchronously transmitted, and the original amplitude attenuation and the original phase offset of all propagation paths are recorded in real time. Then, the phase decoupling processing is performed: the original phase is fully periodically disentangled to eliminate the 2π ambiguity, the normalized phase gradient is calculated, and the dual-frequency differential phase method is used to suppress the multipath effect interference. Subsequently, the joint inversion model is constructed: the phase gradient mapping error is taken as the core constraint (the initial weight α0=0.7), and the amplitude attenuation error is taken as the secondary constraint (the initial weight β0=0.3), and the confidence coefficient is dynamically calculated according to the environmental noise variance σ.

[0139] The regularization iterative algorithm is used to solve the three-dimensional distribution of the dielectric constant, and morphological enhancement is performed on the result: anisotropic diffusion filtering is first performed to preserve the hollow edge, and then the curvature-driven level set method (convex boundary contraction and concave boundary expansion) is used to optimize the three-dimensional profile, and finally the region with a dielectric constant ≤3.0 is extracted as a hollow, and its spatial position, volume and confidence probability are output.

[0140] In a typical implementation (such as tunnel detection), the system can identify a hollow with a diameter ≥0.2m in 10-30 minutes, with a positioning error <5% and a volume error <10%.

[0141] Implementation example;

[0142] Scenario: hollow detection behind tunnel lining.

[0143] Parameter setting:

[0144] Frequency band: f L =100MHz, f H =20MHz;

[0145] Array: 12 channels for transmitting / receiving each, spacing 0.5m;

[0146] Weight initial value: α0=0.7, β0=0.3;

[0147] Hollow criterion less than or equal to 3.0 (air / filler dielectric constant).

[0148] Result:

[0149] Successfully identified a hollow with a diameter ≥0.2m, with a positioning error <5% and a volume estimation error <10%.

[0150] It is to be understood that the terminology used herein is for the purpose of describing particular embodiments only and is not intended to be limiting; it is not intended to exclude myriad other embodiments of the present application that other inventors can develop based on the same general inventive concepts embodied by the described embodiments. That is, although the present application is described in terms of particular embodiments and implementations, it is to be understood that the terminology used is for the purpose of descriptive clarity and that it should be taken in its broadest possible sense. For example, the terms "a", "an", and "the" include both singular and plural referents unless the context clearly dictates otherwise. The terms "comprises", "comprising", "includes", "including" and the like can be used in conjunction with the term "consisting of to include the elements or steps listed after such conjunctive language, but not to the exclusion of other elements or steps. The singular forms "a", "an" and "the" include plural referents unless the context clearly dictates otherwise. The term "plurality" means two or more. The term "consisting essentially of to provide that the composition or process include additional steps, elements, compounds, compositions of matter, or materials not specifically recited. The use of the terms "first", "second", and the like does not imply any particular ordering, but rather are used to denote distinct and separate steps. Where the context requires, the singular forms "a", "an" and "the" include their corresponding plural referents unless the context clearly dictates otherwise. Terms such as "above" and "below" refer to positions relative to the orientation of the figures and are used for purposes of illustration and description only. Terms such as "first" and "second" are used to identify various elements, but the elements should not be limited by these terms. The use of the terms "first" and "second" are not necessarily intended to connote chronological order, but rather serve as labels to distinguish one element from another. The use of the terms "a", "an", and "the" and / or the use of articles in the singular and plural sense are intended to cover a

[0151] While the embodiments of the application have been shown and described herein, it is understood that modifications, substitutions, combinations, and variations of the embodiments can be undertaken by those skilled in the art without departing from the spirit and scope of the present application, which is defined by the following claims and their equivalents.

Claims

1. A method for identifying underground cavities using electromagnetic wave CT based on phase information, characterized by: The method comprises the following steps: Step A1, data collection; An array of transmitting and receiving antennas is deployed in the underground space to transmit electromagnetic wave signals and simultaneously record the amplitude attenuation and phase offset of the electromagnetic waves at each receiving point; Step A2: phase decoupling processing; Perform full-cycle phase unwrapping on the original phase data to eliminate the 2π periodic ambiguity of the phase data and calculate the normalized phase gradient of each propagation path; Step A3, joint inversion modeling; Construct a target optimization function with the phase gradient term as the core constraint and the amplitude attenuation term as the secondary constraint. The specific calculation formula of the target optimization function is: Objective optimization function = minimize the weighted sum of squares of the mapping error from phase gradient to dielectric constant + the weighted sum of squares of the error between the measured and calculated amplitudes; Among them, the phase gradient term weight coefficient α is greater than the amplitude attenuation term weight coefficient β; Step A4: cavity identification; A regularized iterative algorithm is used to solve the spatial distribution of dielectric constant, and abnormal areas with dielectric constants less than or equal to 3.0 are extracted and determined to be underground cavities.

2. The method for identifying underground cavities using electromagnetic wave CT based on phase information according to claim 1, characterized in that: The frequency band of the electromagnetic wave signal in step A1 is a wideband frequency band from 10 MHz to 1 GHz.

3. The method for identifying underground cavities using electromagnetic wave CT based on phase information according to claim 1, characterized in that: The specific calculation formula of the normalized phase gradient in step A2 is: Normalized phase gradient = (observed phase - reference phase) / (angular frequency x propagation path length).

4. The method for identifying underground cavities using electromagnetic wave CT based on phase information according to claim 3, characterized in that: In the phase decoupling process of step A2, a dual-frequency differential phase method is used to eliminate system errors. The specific operation steps of the dual-frequency differential phase method are as follows: Step S1, synchronously transmitting low frequency f1 and high frequency f2 signals, wherein the low frequency f1 is greater than or equal to 5 times the high frequency f2; Step S2: Calculate the differential phase: ΔΦ_diff=Φ(f1)-k·Φ(f2); Where k is the frequency ratio coefficient.

5. The method for identifying underground cavities using electromagnetic wave CT based on phase information according to claim 1, characterized in that: In the joint inversion model in step A3, a phase-amplitude confidence weight factor is introduced. The specific operation steps are as follows: Step D1: define the confidence coefficient. The specific calculation formula of the confidence coefficient is: γ=e^(–σ 2 / ΔΦ_norm 2 ); Where, σ is the variance of the ambient noise; Step D2: Set a dynamic adjustment objective function. The specific calculation formula of the dynamic adjustment objective function is: α=γ·α0,β=(1–γ)β0; Step D3: Dynamic weight allocation; Set the initial weights α0 for the phase term and β0 for the amplitude term to ensure that α0 is greater than β0.

6. The method for identifying underground cavities using electromagnetic wave CT based on phase information according to claim 1, characterized in that: In step A4, the spatial distribution of the dielectric constant is solved by using a morphological enhancement algorithm, which specifically includes the following steps: Step F1, performing anisotropic diffusion filtering on the spatial distribution of dielectric constant to retain the edges of the voids; Step F2, extracting the three-dimensional boundary of the cavity based on the curvature-driven level set method; Step F3: Output the cavity position, volume and confidence probability.

7. The method for identifying underground cavities using electromagnetic wave CT based on phase information according to claim 6, characterized in that: In step F2, the specific steps for extracting the three-dimensional boundary of the hole are as follows: Step F21, curvature calculation; Calculate the average curvature of the current boundary surface; Perform gradient calculation on the boundary function in three-dimensional space to obtain the directional change rate of each point; The divergence is further calculated based on the gradient result to obtain the exact curvature value. The curvature of the convex area is positive and that of the concave area is negative. Step F21, execution of evolution rules; Dynamically adjust the boundary movement direction according to the curvature value: The area of ​​positive curvature is convex: the boundary shrinks inward; The negative curvature region is concave: the boundary expands outward; The region of zero curvature is flat: the boundary remains stationary.