A method for inverting the dielectric constant of the lunar subsurface based on diffraction wave velocity analysis

Abstract

The present invention belongs to the technical field of lunar subsurface material property inversion, and specifically relates to a lunar subsurface dielectric constant inversion method based on diffraction wave velocity analysis, the method comprising the following steps: reading, screening and merging lunar radar data to obtain effective detection data; preprocessing the effective detection data to obtain high signal-to-noise ratio data; performing diffraction wave velocity analysis on the preprocessed effective detection data to obtain the best offset velocity; performing Dix conversion on the obtained best offset velocity to obtain the layer velocity; and estimating the dielectric constant based on the layer velocity. The method solves the problems of accuracy and large amount of calculation in the traditional hyperbola fitting method when processing lunar radar data with low signal-to-noise ratio, complex morphology and large data volume, and realizes fast, automatic and high-resolution dielectric constant calculation, which has good applicability to deep space detection data.

Description

Technical Field

[0001] The present invention belongs to the technical field of lunar subsurface material property inversion, and specifically relates to a lunar subsurface dielectric constant inversion method based on diffraction wave velocity analysis. Background Art

[0002] As the closest celestial body to the Earth, the study of the evolutionary history of the Moon's morphological structure and material composition plays a key role in understanding the Earth and terrestrial planets. During its evolution, the Moon has been continuously acted upon by both internal and external forces, such as volcanic eruptions, weathering erosion, and meteorite impacts. These forces have continuously transformed the lunar surface, but the history of transformation can be traced back through the subsurface materials and structures. Especially in the Antarctic Aiken Basin, as the largest, deepest, and oldest impact mechanism on the Moon, it is most likely to excavate deep materials such as the lunar crust and even the lunar mantle, making it an important research area for lunar evolution analysis.

[0003] On January 3, 2019, the Chang'e-4 satellite landed on the back of the moon, marking the first soft landing on the back of the moon by humans. The lander then separated from the rover, and the rover began to patrol and detect the lunar surface. The lunar radar carried by the rover is responsible for detecting the structure and material of the subsurface layer along the way. The lunar radar includes dual-channel antennas with center frequencies of 60MHz and 500MHz, respectively. The 500MHz detection depth is shallow but the resolution is high, which is more suitable for the inversion of material properties. Compared with the structure, the properties contain more detailed information and are more suitable for understanding the properties of the material. In particular, the dielectric constant is controlled by the material composition and density, which directly affects the wave velocity and signal attenuation. Therefore, it is an important channel to understand the material on the lunar surface by inferring the dielectric constant through the wave velocity and inferring the material composition.

[0004] At present, the dielectric constant inversion methods of radar data mainly include metal plate method, impedance inversion method, dual antenna delay method and hyperbola fitting method. Among them, the hyperbola fitting method can best reflect the spatial anisotropy of the medium and is widely used. The principle of the hyperbola method is that the point target will diffract after receiving the radar signal, and an obvious hyperbola shape will appear on the radar profile. The hyperbola shape is determined by the depth of the target and the wave velocity. Therefore, the wave velocity can be determined by fitting the hyperbola, and then the dielectric constant can be inverted.

[0005] The traditional hyperbola fitting method is to manually pick up hyperbolas, which is subjective in judging the number and shape, and has a large workload. There are many sputtering objects in the lunar subsurface, and the number of hyperbolas is large and overlaps each other. In addition, the environmental and instrument noises in the lunar radar detection process can easily cause hyperbola deformation and distortion. The traditional hyperbola fitting method has many obstacles such as accuracy and large amount of calculation when processing lunar radar data with low signal-to-noise ratio, complex shape and large data volume. Summary of the invention

[0006] The purpose of the present invention is to provide a method for inverting the dielectric constant of the lunar subsurface layer based on diffraction wave velocity analysis, so as to solve the problems of accuracy and large amount of calculation in the traditional hyperbola fitting method when processing lunar radar data with low signal-to-noise ratio, complex morphology and large data volume.

[0007] The present invention is achieved in this way.

[0008] A method for inverting the dielectric constant of the lunar subsurface layer based on diffraction wave velocity analysis, the method comprising the following steps:

[0009] S1 reads, filters and merges the lunar radar data to obtain effective detection data;

[0010] S2 pre-processes the effective detection data to obtain high signal-to-noise ratio data;

[0011] S3 performs diffraction wave velocity analysis on the effective detection data after preprocessing to obtain the optimal migration velocity;

[0012] S4 performs Dix transformation on the obtained optimal migration velocity to obtain the interval velocity;

[0013] S5 estimates the dielectric constant based on the layer velocity.

[0014] Furthermore, the diffraction wave velocity analysis in S3 includes:

[0015] S3-1: Use the plane wave destruction method to separate the diffraction waves;

[0016] S3-2: Use the velocity continuation method to establish the three-dimensional velocity spectrum of the diffraction wave;

[0017] S 3-3: Using negative entropy as a measure of the focusing degree of the diffraction wave, calculate the focusing velocity spectrum corresponding to the three-dimensional velocity spectrum;

[0018] S 3-4: From the focused velocity spectrum, the offset velocity is obtained using an automatic picking algorithm.

[0019] Furthermore, the plane wave destruction method described in S3-1 is implemented based on the following formula:

[0020] ,

[0021] in, is the spatial position of the radar antenna, For walking time, For TRON, For and The local inclination angle of the in-phase axis is the local inclination angle of the reflected wave, which is a constant, and the local inclination angle of the diffracted wave is nonlinear. By calculating the local inclination angle of the in-phase axis, the reflected wave and the diffracted wave can be separated.

[0022] Furthermore, the velocity extension method described in S3-2 is implemented in the frequency-wave number domain based on the following formula:

[0023] ,

[0024] in, is the initial velocity, is the migration velocity to be selected, is the Fourier transform result of the initial wave field, is the wave number, is the frequency, is the wave field corresponding to the migration velocity to be selected. During the calculation, the fast Fourier forward and inverse transform and the phase shift operator are used to solve the migration velocity to be selected. The set of two-dimensional diffraction wave migration images corresponding to all the migration velocities to be selected constitutes a three-dimensional velocity spectrum.

[0025] Furthermore, the negative entropy described in S3-3 is realized based on the following formula:

[0026] ,

[0027] ,

[0028] in, is the pixel number of the time dimension in the sliding window, is the pixel number of the spatial dimension in the sliding window, and is the size of the sliding window, That is, the total number of pixels in the sliding window, is the Hilbert transform result of the original wave field amplitude, It is a sliding window The Hilbert transform result of the original wave field amplitude at the position, yes The normalized amplitude within the sliding window, Represents negative entropy, and the calculation formula of negative entropy is simplified to:

[0029] ,

[0030] ,

[0031] in, Represents AGC automatic gain, is the sampling point number of each data channel, yes At sampling point The amplitude of the velocity spectrum is calculated, and the negative entropy of each point in the three-dimensional velocity spectrum is obtained to obtain the focused velocity spectrum.

[0032] Furthermore, the calculation method of the automatic picking algorithm in S3-4 includes: extracting the negative entropy corresponding to each track from the focusing velocity spectrum to obtain the optimal offset velocity that can maximize the following formula: :

[0033] ,

[0034] in, is the optimal offset speed The corresponding trace in the corresponding two-dimensional diffraction wave migration image is The negative entropy of the moment, i.e. the degree of focus, Indicates the start time of radar recording, Indicates the end time of radar recording. The optimal displacement speed The integral of the corresponding focusing degree in the time dimension, that is, the sum of the focusing degrees of each channel, finds the optimal migration velocity corresponding to the highest focusing value at all times for each channel data, and combines the optimal migration velocities of each channel to form the initial velocity profile;

[0035] The initial velocity profile is smoothed in time and space domains using a shaping regularization window to remove outliers with excessively large gradients, thereby obtaining an optimal migration velocity profile that is heterogeneous in time and space but changes smoothly. The degree of smoothing is adjusted by the smoothing radius.

[0036] Furthermore, the Dix conversion in step 4 converts the optimal migration velocity into the layer velocity according to the following formula:

[0037] ,

[0038] in, For the The layer velocity of the layer, To the The two-way travel time of the layer is Indicates The optimal migration speed of the layer, For the round-trip travel time to the n-1th floor, use Indicates Optimal offset speed for the layer.

[0039] Furthermore, the dielectric constant in step 5 is estimated using the following formula:

[0040] ,

[0041] in, is the dielectric constant, is the propagation speed of radar waves in vacuum, When calculating the layer, correspond .

[0042] Compared with the prior art, the present invention has the following beneficial effects:

[0043] The method of the present invention separates the diffraction waves that record velocity information, eliminating the interference of plane waves and noise; utilizes velocity continuation, focusing analysis and automatic picking methods to quickly, automatically and objectively perform velocity, reduces the subjective uncertainty of visual identification and manual picking, and can meet the needs of high-precision velocity analysis, and realizes fast, automatic and high-resolution dielectric constant calculation, has good applicability to deep space exploration data, and is beneficial to subsequent subsurface layer division and attribute calculation research. BRIEF DESCRIPTION OF THE DRAWINGS

[0044] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the drawings required for use in the description of the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present application. For ordinary technicians in this field, other drawings can be obtained based on these drawings without creative work.

[0045] Figure 1 is a flow chart of a method for inverting the dielectric constant of the lunar subsurface layer based on diffraction wave velocity analysis provided by an embodiment of the present invention;

[0046] Figure 2 1 is a comparison diagram of Chang'e-4 lunar radar data before and after preprocessing provided by an embodiment of the present invention, wherein (a) is before processing, and (b) is after processing;

[0047] Figure 3 The reflected wave and the diffracted wave separated from the Chang'e-4 lunar radar data provided by the embodiment of the present invention include (a) the reflected wave and (b) the diffracted wave;

[0048] Figure 4 This is a stratification result diagram of the Chang'e-4 lunar radar data provided by an embodiment of the present invention, wherein (a) the stratification is based on the preprocessing result, and (b) the stratification is based on the inverted dielectric constant result. DETAILED DESCRIPTION

[0049] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the embodiments. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0050] The basic principle of a lunar subsurface dielectric constant inversion method based on diffraction wave velocity analysis is: separate the diffraction wave containing velocity information and remove the interference of plane waves and noise that are irrelevant to velocity; perform velocity analysis on the diffraction wave to obtain the optimal offset velocity at all positions of the radar profile; use Dix transformation to convert the offset velocity into layer velocity; and use the layer velocity to estimate the dielectric constant.

[0051] See also Figure 1 FIG. 1 is a flow chart of a method for inverting the dielectric constant of the lunar subsurface layer based on diffraction wave velocity analysis of the present invention, and the method comprises the following steps:

[0052] Step 1: read, filter and merge the lunar radar data to obtain effective detection data;

[0053] In this embodiment of the present invention, the selected lunar radar data is the data of Chang'e-4, Channel 2-500MHz antenna, 2B level, and the first 26 lunar days, which is released by the National Astronomical Observatory of the Chinese Academy of Sciences. Its format is the PDS4 format developed by NASA for storing planetary exploration missions. According to the data label format, radar echo data and rover position and other information can be read.

[0054] The rover equipped with lunar radar has a variety of detection methods, such as mobile detection, on-the-spot turning detection, fixed-point detection, etc. Therefore, the original radar data will be repeatedly detected at a fixed point when it stops or turns, resulting in data redundancy. At the same time, due to the influence of the machine or the environment, there are also abnormal road conditions, so data screening and splicing are required.

[0055] Step 2, preprocessing the effective detection data to obtain high signal-to-noise ratio data; including time zero point correction of the effective detection data to eliminate the influence of the start-up time difference between the receiving antenna and the transmitting antenna in the radar antenna; performing event axis correction on the data after time zero point correction, and adjusting the channel gather according to the continuity of the event axis; performing data gain on the data after event axis correction, and then removing horizontal background noise, bandpass filtering, reducing high-frequency and low-frequency noise outside the transmission signal frequency band, and random noise removal, reducing random noise similar to the transmission signal frequency; after the noise is removed, the radar profile is resampled according to the position information of part of the channel gather to supplement the position information.

[0056] in:

[0057] Time zero correction: Since the radar's receiving antenna is turned on 28.203ns earlier than the transmitting antenna, it is necessary to delete the invalid data recorded within the first 28.203ns.

[0058] Phase axis correction. During data collection, the lunar radar has multiple sleep and restart times, and there are slight deviations in the start time of collection of different data segments, which will cause longitudinal misalignment of the image. Therefore, it is necessary to adjust the gather for the misalignment phenomenon to reduce the error.

[0059] Data gain. Radar waves are lost and attenuated during propagation, and the energy in deep areas becomes lower and lower, which is not conducive to later waveform recognition and analysis. Therefore, automatic gain control AGC is used to compensate the amplitude gain of deep echoes to control the amplitude energy to be relatively balanced across the entire radar profile.

[0060] Background removal: Direct waves, ground reflections, equipment influences, etc. will generate horizontal background noise. By removing the average channel, the background can be suppressed and the abnormal body information can be highlighted.

[0061] Bandpass filtering. During radar data collection, various types of noise will also be recorded. By setting a filter for filtering, information near the center frequency can be retained and noise outside the transmission signal band can be reduced;

[0062] Random noise removal. For some noises similar to the frequency band of the transmitted signal, bandpass filtering cannot eliminate them. Most of them are random noises. Methods such as streaming orthogonal adaptive predictive filtering can be used to remove random noise and reconstruct signals.

[0063] Position information supplement. Although the patroller does not move at a uniform speed during detection, the corresponding waypoint information during data collection is recorded in the 2B level data. Therefore, the radar profile can be interpolated at equal distances based on the waypoint position information to obtain a radar profile with equal spacing and position information.

[0064] In an embodiment of the present invention, Figure 2 is the comparison of radar profiles before and after preprocessing. Figure 2 In the middle (a), after preprocessing, Figure 2 Figure (b) in the figure. The horizontal axis is the spatial distance from the landing point of the rover when the radar antenna detects, and the vertical axis is the two-way travel time of the detection data. The signal-to-noise ratio of the radar profile before preprocessing is low, making it difficult to identify effective signals. After preprocessing, the noise is suppressed, the effective signal is enhanced, and the position information is supplemented, so that the main stratigraphic structure can be clearly identified.

[0065] Step 3, performing diffraction wave velocity analysis on the pre-processed effective detection data to obtain the optimal migration velocity;

[0066] The basic principle of velocity inversion using diffraction waves is that when the radar signal propagates in the medium and encounters anomalies, reflection and diffraction phenomena will occur. The reflected wave is formed by the medium layer and records the layer information of the medium. The diffraction wave is formed by anomalies such as rocks and faults. Location As it moves forward, the radar waves encounter When there is an abnormal body, reflection will occur. is the distance from the radar antenna to the anomaly, is the double-layer travel time of the radar wave, assuming that the radar wave speed is , then the radar echo satisfies:

[0067] ,

[0068] The equation can be rewritten as:

[0069] ,

[0070] because and Fixed, then and The hyperbolic relationship is satisfied, and the hyperbolic shape is controlled by the radar wave velocity. Therefore, the radar wave velocity can be inverted by analyzing the hyperbolic shape.

[0071] Step 3-1: Use the plane wave destruction method to separate the diffraction waves.

[0072] The present invention separates the diffraction wave from the reflection wave in this step. The diffraction wave has a weaker amplitude than the reflection wave, so the velocity analysis using the separated diffraction wave has higher accuracy and credibility. The plane wave destruction method is a signal extraction method based on linear prediction. The reflection wave has a strong correlation in space, and the local inclination is a constant. The diffraction wave has a poor correlation and the local inclination is nonlinear. Therefore, the reflection wave can be predicted and eliminated by calculating the local inclination of the phase axis, and the diffraction wave field is retained. The plane wave destruction method is implemented based on the following formula:

[0073] ,

[0074] in, is the spatial position of the radar antenna, For walking time, For TRON, For and The local inclination angle of the in-phase axis is constant, and the local inclination angle of the diffracted wave is nonlinear. By calculating the local inclination angle of the in-phase axis, the reflected wave and the diffracted wave are separated. The wave field separation result is shown in Figure 3 As shown, Figure 3 (a) is the reflected wave, Figure 3(b) in the figure is the diffraction wave. The horizontal axis is the spatial distance from the landing point of the rover when the radar antenna detects, and the vertical axis is the two-way travel time of the detection data. The separated reflection waves show the continuous layered structure more prominently, and the separated diffraction waves show the details of the hyperbolic shape more clearly, providing better analysis data for velocity analysis.

[0075] Step 3-2: Use the velocity continuation method to establish the three-dimensional velocity spectrum of the diffraction wave:

[0076] In this step, the present invention calculates the radar profile in the time-space domain. Expand the velocity dimension, that is, use a series of velocities to migrate the diffraction wave field and obtain the migration image corresponding to each velocity, that is, the three-dimensional velocity spectrum , which is the basic data for subsequent search for the optimal migration speed.

[0077] The traditional method of calculating the three-dimensional velocity spectrum is to use the test speed set at a certain step size to offset the radar profile one by one, which is computationally intensive and has low computational efficiency. The present invention uses a simple algorithm based on fast Fourier transform to perform velocity extension, uses a single offset to determine the initial offset image, and then uses a velocity-related phase shift filter to calculate the offset results of other velocities. The calculation formula in the frequency-wavenumber domain is:

[0078] ,

[0079] in, is the initial velocity, is the migration velocity to be selected, is the Fourier transform result of the initial wave field, is the wave number, is the frequency, is the wave field corresponding to the migration velocity to be selected. During the calculation, the fast Fourier forward and inverse transform and the phase shift operator are used to solve the migration velocity to be selected. The set of two-dimensional diffraction wave migration images corresponding to all the migration velocities to be selected constitutes a three-dimensional velocity spectrum.

[0080] The specific calculation steps are:

[0081] (1) At initial speed Performing migration processing on the radar profile image;

[0082] (2) The offset image is Regular grid to Regular grid;

[0083] (3) Performing fast Fourier transform on the result of step (2);

[0084] (4) Multiply the result of step (3) by the all-pass phase shift filter ;

[0085] (5) Performing an inverse Fourier transform on the result of step (4);

[0086] (6) Convert the result of step (4) from Regular grid to Regular grid.

[0087] The set of two-dimensional diffraction wave migration images corresponding to all the migration velocities to be selected constitutes a three-dimensional velocity spectrum.

[0088] Step 3-3: Using negative entropy as a measure of the focusing degree of the diffraction wave, calculate the focusing velocity spectrum corresponding to the three-dimensional velocity spectrum:

[0089] The present invention analyzes the degree of concentration of the three-dimensional velocity spectrum in this step. Only the correct offset velocity can make the hyperbola converge and return to its position, and the focusing degree is the highest at this time. Therefore, the accuracy of the focusing degree directly affects the accuracy of the velocity analysis.

[0090] The present invention uses the negative entropy method to perform focus analysis, that is, the negative entropy of the image is used as a quantitative indicator of the focus degree. Negative entropy is the negative change of entropy and is a measure of the orderliness of the system. For a 2D image, the higher the focus degree, the greater the negative entropy, and the lower the focus degree, the smaller the negative entropy. In the specific operation, a sliding window is used for calculation, and the formula is as follows:

[0091] in, is the pixel number of the time dimension in the sliding window, is the pixel number of the spatial dimension in the sliding window, and is the size of the sliding window, That is, the total number of pixels in the sliding window, is the Hilbert transform result of the original wave field amplitude, It is a sliding window The Hilbert transform result of the original wave field amplitude at the position, yes The normalized amplitude within the sliding window, Represents negative entropy.

[0092] However, setting a sliding window and calculating pixel by pixel requires a large amount of calculation, so the calculation method is simplified as follows:

[0093] ,

[0094] ,

[0095] in, Represents AGC automatic gain, is the sampling point number of each data channel, yes At sampling point The amplitude of the focused velocity spectrum is obtained by calculating the negative entropy for each point in the three-dimensional velocity spectrum.

[0096] This operation can greatly improve the computational efficiency and realize the pixel-by-pixel focusing measurement of all velocity corresponding offset results, namely the focusing velocity spectrum.

[0097] Step 3-4: From the focused velocity spectrum, use the automatic picking algorithm to obtain the best offset velocity:

[0098] The present invention performs automatic speed picking in this step. After obtaining the focusing results of the offset images corresponding to each offset speed, i.e., the focusing speed spectrum, the optimal offset speed can be found accordingly. Theoretically, in the focusing result slices of each channel, the optimal offset speed It should satisfy the following equation to be maximized:

[0099] ,

[0100] in, is the optimal offset speed The corresponding trace in the corresponding two-dimensional diffraction wave migration image is The negative entropy of the moment is the degree of focus. Indicates the start time of radar recording, Indicates the end time of radar recording. The optimal displacement speed The integral of the corresponding focusing degree in the time dimension, that is, the sum of the focusing degrees of each channel, finds the migration speed corresponding to the highest focusing value at all times for each channel data as the optimal migration speed, and merges each channel to form an initial velocity profile; at this time, the initial velocity profile is mostly rough and changes dramatically in space.

[0101] In order to make the velocity results closer to reality, the present invention uses the shaping regularization method proposed by Fomel to perform regularization constraints on the velocity results in the time domain and space domain, remove outliers with excessive gradients, and obtain the optimal migration velocity profile that is heterogeneous in time and space but changes smoothly, where the degree of smoothing is adjusted by the smoothing radius.

[0102] Step 4, converting the obtained optimal migration velocity to obtain the layer velocity;

[0103] The best migration velocity estimated based on the diffraction wave velocity analysis is the root mean square velocity, that is, the average velocity of the medium between the ground and the diffraction point, not the true velocity of each pixel. According to lunar physiognomy, the velocity of the lunar subsurface is related to density, etc., and usually changes layer by layer. Therefore, it is necessary to convert the root mean square velocity into the true velocity of each layer before the subsequent attribute inversion can be carried out.

[0104] The present invention introduces Dix conversion, and the Dix conversion calculation formula is as follows:

[0105] ,

[0106] ,

[0107] in, For the The layer velocity of the layer, To the The two-way travel time of the layer is Indicates The displacement speed of the layer, To the The two-way travel time of the layer is Indicates The offset velocity of the layer. The result is the two-dimensional layer velocity result corresponding to the radar profile.

[0108] Step 5, it is helpful to estimate the layer velocity and obtain the dielectric constant;

[0109] According to lunar physiognomy, the typical mass magnetic susceptibility of lunar soil is in the range of 100-10000*10 -9 m 3 / kg, so the relative magnetic permeability is very close to 1. At the same time, the electrical conductivity is also very small, about 10 at 450MHz. -3 -10 -4 S / m, therefore, the lunar soil can be regarded as an ideal dielectric, and the radar wave speed is mainly determined by the relative dielectric constant of the medium, that is:

[0110] ,

[0111] in, is the relative dielectric constant, is the propagation speed of radar waves in vacuum (3*10 8 m / s), for When calculating the layer, The corresponding Therefore, the calculation formula for the relative dielectric constant is:

[0112] ,

[0113] The calculation result is the dielectric constant distribution under the moving track of the patroller. Figure 4 As shown, Figure 4 (a) The layering with the preprocessing result as the background, Figure 4The dielectric constant result inverted in (b) is the background layering. In the figure, the horizontal axis is the spatial distance from the rover landing point when the radar antenna detects, the vertical axis is the two-way travel time of the detection data, the red solid line is the layer interface division result, the white dashed line is the sub-layer interface division result, and the white solid line is the interface of the hidden pit. It can be seen that the interface division result based on the radar profile reflection intensity can match the interface division result based on the interface constant distribution well. In the vertical direction, the first layer of 0-150ns is the lunar soil layer formed by surface weathering, which contains a small amount of stones on average and has a dielectric constant between 2.3 and 4; the second layer is the sputtering layer, which contains multiple sublayers and is destroyed by impact craters at a spatial distance of 120-280. Each sublayer includes an upper layer of sputtering materials with a low degree of weathering and a lower layer of lunar soil with a high degree of weathering. The dielectric constant of the sputtering materials is between 4 and 5.5, and the lunar soil is slightly larger than that of the lunar surface. The dielectric constant of the upper part of the impact crater is large, and the value of the lower layer is small, that is, the sputtering materials filled in the surface layer have not been completely weathered and are covered by new sputtering materials in the later period; the third layer is the basalt layer, and the dielectric constant is significantly higher than that of the first two layers. Among them, the upper sublayer of 400-500ns has a large dielectric constant at the top and contains multiple low dielectric constant areas, which are speculated to be the weathered basalt layer. The lower sublayer of 500-612ns has a dielectric constant mostly higher than 6, which is speculated to be a more complete basalt layer. The above structural stratification and dielectric constant distribution are highly consistent with the vast majority of current research results, so the method has a high degree of credibility.

[0114] Compared with the traditional lunar dielectric constant calculation method, the lunar subsurface dielectric constant inversion method based on diffraction wave velocity analysis proposed in the present invention separates the diffraction wave that records the velocity information and eliminates the interference of plane waves and noise; the velocity is quickly, automatically and objectively calculated by using velocity continuation, focusing analysis and automatic picking methods, which reduces the subjective uncertainty of visual identification and manual picking, meets the needs of high-precision velocity analysis, and realizes fast, automatic and high-resolution dielectric constant calculation. It has good applicability to deep space exploration data and is beneficial to the subsequent subsurface layer division and attribute calculation research.

[0115] Those skilled in the art will appreciate that the present invention can be used for the inversion of the subsurface dielectric constant of the Chang'e-4 lunar radar, and can also be applied to other deep space detection radars.

[0116] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions and improvements made within the spirit and principles of the present invention should be included in the protection scope of the present invention.

Claims

1. A method for inverting the dielectric constant of the lunar subsurface layer based on diffraction wave velocity analysis, characterized in that: The method comprises the following steps: S1 reads, filters and merges the lunar radar data to obtain effective detection data; S2 preprocesses the effective detection data to obtain high signal-to-noise ratio data; including time zero point correction of the effective detection data to eliminate the influence of the start-up time difference between the receiving antenna and the transmitting antenna in the radar antenna; the data after time zero point correction is corrected for the event axis, and the channel gather is adjusted according to the continuity of the event axis; the data after event axis correction is performed for data gain, and then the horizontal background noise is removed in turn, and bandpass filtering is performed to reduce the high-frequency and low-frequency noise outside the transmission signal frequency band, and random noise is removed to reduce the random noise similar to the transmission signal frequency; the radar profile is resampled according to the position information of part of the channel gather after noise removal to supplement the position information, and the radar profile is interpolated at equal distances according to the waypoint position information to obtain the radar profile with equal spacing and position information; S3 performs diffraction wave velocity analysis on the effective detection data after preprocessing to obtain the optimal migration velocity; S4 performs Dix transformation on the obtained optimal migration velocity to obtain the interval velocity; S5 estimates the dielectric constant based on the layer velocity; The diffraction wave velocity analysis described in S3 includes: S3-1: Use the plane wave destruction method to separate the diffraction waves; S3-2: Use the velocity continuation method to establish the three-dimensional velocity spectrum of the diffraction wave; S 3-3: Using negative entropy as a measure of the focusing degree of the diffraction wave, calculate the focusing velocity spectrum corresponding to the three-dimensional velocity spectrum; S 3-4: From the focused velocity spectrum, the offset velocity is obtained using the automatic picking algorithm; The velocity extension method described in S3-2 is implemented in the frequency-wavenumber domain based on the following formula: ， in, is the initial velocity, is the migration velocity to be selected, is the Fourier transform result of the initial wave field, is the wave number, is the frequency, is the wave field corresponding to the migration velocity to be selected. In the calculation, the migration velocity to be selected is solved by using the fast Fourier forward and inverse transform and the phase shift operator. The set of the two-dimensional diffraction wave migration images corresponding to all the migration velocities to be selected constitutes a three-dimensional velocity spectrum. The negative entropy described in S3-3 is realized based on the following formula: ， ， in, is the pixel number of the time dimension in the sliding window, is the pixel number of the spatial dimension in the sliding window, and is the size of the sliding window, That is, the total number of pixels in the sliding window, is the Hilbert transform result of the original wave field amplitude, It is a sliding window The Hilbert transform result of the original wave field amplitude at the position, yes The normalized amplitude within the sliding window, Represents negative entropy, and the calculation formula of negative entropy is simplified to: ， ， in, Represents AGC automatic gain, is the sampling point number of each data channel, yes At sampling point The amplitude of the velocity spectrum is calculated, and the negative entropy of each point in the three-dimensional velocity spectrum is obtained to obtain the focused velocity spectrum; The calculation method of the automatic picking algorithm in S3-4 includes: extracting the negative entropy corresponding to each channel from the focusing velocity spectrum, and obtaining the optimal offset velocity that can maximize the following formula : ， in, is the optimal offset speed The corresponding trace in the corresponding two-dimensional diffraction wave migration image is The negative entropy of the moment, i.e. the degree of focus, Indicates the start time of radar recording, Indicates the end time of radar recording. The optimal displacement speed The integral of the corresponding focusing degree in the time dimension, that is, the sum of the focusing degrees of each channel, finds the optimal migration velocity corresponding to the highest focusing value at all times for each channel data, and combines the optimal migration velocities of each channel to form the initial velocity profile; The initial velocity profile is smoothed in time and space domains using a shaping regularization window to remove outliers with excessively large gradients, thereby obtaining an optimal migration velocity profile that is heterogeneous in time and space but changes smoothly. The degree of smoothing is adjusted by the smoothing radius.

2. The method for inverting the dielectric constant of the lunar subsurface layer based on diffraction wave velocity analysis according to claim 1, characterized in that: The plane wave destruction method described in S3-1 is based on the following formula: ， in, is the spatial position of the radar antenna, For walking time, For TRON, For and The local inclination angle of the event axis is the local inclination angle of the reflected wave, which is a constant, and the local inclination angle of the diffracted wave is nonlinear. By calculating the local inclination angle of the event axis, the reflected wave and the diffracted wave can be separated.

3. The method for inverting the dielectric constant of the lunar subsurface layer based on diffraction wave velocity analysis according to claim 1, characterized in that: The Dix conversion described in step 4 converts the optimal migration velocity into the interval velocity according to the following formula: ， in, For the The layer velocity of the layer, To the The two-way travel time of the layer is , Indicates The optimal migration speed of the layer, For the round-trip travel time to the n-1th floor, use Indicates Optimal offset speed for the layer.

4. The method for inverting the dielectric constant of the lunar subsurface layer based on diffraction wave velocity analysis according to claim 3 is characterized in that: The dielectric constant described in step 5 is estimated using the following formula: ， in, is the dielectric constant, is the propagation speed of radar waves in vacuum, When calculating the layer, correspond .

Images