An ice radar imaging method based on ice surface curvature

Abstract

An ice radar imaging method based on ice surface curvature belongs to the field of ice radar data processing. It includes: a flat ice surface propagation distance generation module, a high curvature ice surface propagation distance generation module, a propagation distance difference generation module, a secondary phase error calculation module, a subaperture length calculation module and an ice radar data imaging module. In order to solve the imaging quality of data collected by airborne ice radar in a high curvature ice surface environment, an ice radar imaging method based on ice surface curvature is proposed. This method improves the imaging quality by correcting the distance parameters of the ice radar imaging method in a high curvature ice surface environment. Thereby, more accurate results are obtained for glacier analysis.

Description

Technical Field

[0001] The invention relates to the field of ice radar data processing, and in particular to an ice radar imaging method based on ice surface curvature. Background Art

[0002] Glaciers are an important part of the Earth's environmental system, reflecting and influencing changes in global climate and sea level. The ice sheets and ice shelves in the polar regions are the most sensitive to global climate change and respond most dramatically, acting as "indicators" and "amplifiers" of climate change. Taking the Antarctic region as an example, its ice reserves account for about 90% of the world's total glaciers and about 70% of the world's freshwater resources. If the Antarctic ice sheet melts completely, its contribution to sea level rise could reach nearly 58 meters.

[0003] The analysis of glaciers is mainly carried out through glacier stereotomography. Due to its strong penetration and high-precision ranging capabilities, ice radar has become the core technology of glacier stereotomography. Ice radar (Ice-penetrating radar), also known as radio-echo sounding radar (Radio-echo sounder) or ice radar (Ices sounding radar), is a geophysical detection system based on electromagnetic wave theory that studies the characteristics of ice and snow media through radar echoes. The weak attenuation of electromagnetic waves propagating inside the ice sheet and the good stratification and homogeneity of the ice sheet are the theoretical basis for ice radar to detect polar ice sheets.

[0004] In terms of frequency-domain focused ice radar imaging algorithms, the earliest SAR method implemented was the matched filter method proposed by Legarsky et al., which correlates the echo response of a reference point target with radar data from the Greenland ice sheet. This method has recently been applied to ice sheet data from West Antarctica. Alley et al. proposed the fk ice radar imaging processing algorithm, which is based on the algorithm of inverting radar echoes to their target sources in the ice. This method has also been applied to data from Greenland and Antarctica.

[0005] In terms of time-domain focused ice radar imaging algorithms, Kusk and Dall applied the time-domain back-projection method to ice radar data. This method can process non-uniformly sampled data, but the calculation time is long. Wu et al. used the time-domain multi-subaperture back-projection algorithm for imaging, providing accurate strip imaging of the ice sheet surface and bedrock.

[0006] However, airborne ice radars usually face the problem of large fluctuations in the ice surface in the detection area. Existing ice radar imaging algorithms do not fully consider this problem faced by airborne ice radars.

[0007] In order to solve the imaging quality problem of data collected by airborne ice radar in high-curvature ice environment, an ice radar imaging method based on ice surface curvature is proposed. This method improves the imaging quality by correcting the distance parameters of the ice radar imaging method in high-curvature ice environment, thereby achieving more accurate results for glacier analysis. Summary of the invention

[0008] In order to overcome the error problem caused by the ice surface curvature in the existing airborne ice radar imaging method, the present invention proposes an ice radar data imaging method based on ice surface curvature. The advantage of this method is that it uses the ice surface curvature data in the area as a reference to correct the propagation distance error caused by the ice surface curvature in the area, and improves the imaging quality of ice radar data in high curvature ice environments. The overall flow chart of the method is as follows Figure 1 As shown, it is mainly divided into six modules, namely, flat ice surface propagation distance generation module, high curvature ice surface propagation distance generation module, propagation distance difference generation module, secondary phase error calculation module, sub-aperture length calculation module and ice radar data imaging module.

[0009] Flat ice surface propagation distance generation module: Since this method needs to use the distance difference between the flat ice surface model and the high curvature ice surface model to calculate the sub-aperture length, it is necessary to input the aircraft flight altitude h, radar detection depth d, and the dielectric constant ε of the ice medium. i Calculate the distance from the radar to the point target under the flat ice model. The flowchart of this module is as follows: Figure 2 As shown in Figure 1, it includes three parts: calculating the distance the signal propagates in the air, the distance the signal propagates in the ice medium, and the refraction angle of the ice surface. In the process of generating the flat ice surface propagation distance, the propagation distance R f The propagation distance R in the air fair and the propagation distance R in the ice medium fice The sum of (as shown in Formula 1):

[0010] R f =R fair +R fice *η ice (Formula 1)

[0011] Where η ice is the refractive index of the air-ice interface. When refraction occurs on the ice surface, the incident angle θ i and the refraction angle θ o Satisfies Snell's refraction law (as shown in formula 2)

[0012]

[0013] where η iceis the refractive index of the ice medium, and the refractive index and the relative dielectric constant ε of the air-ice interface satisfy the following relationship (as shown in Formula 3):

[0014] η= (Formula 3)

[0015] High curvature ice surface propagation distance generation module: Different from the flat ice surface model propagation distance calculation, when calculating the high curvature plane model propagation distance, in addition to the input aircraft flight altitude, radar detection depth, and dielectric constant of the ice medium calculated in the flat ice surface model, the slope angle θ of the target area is also required to calculate the distance from the radar to the point target. The flowchart of this module is as follows Figure 3 As shown in Figure 1, it includes three steps: calculating the distance the signal propagates in the air, the distance the signal propagates in the ice medium, and the refraction angle of the ice surface. In the process of generating the propagation distance of the high curvature ice surface, the propagation distance R c The propagation distance R in the air cair and the propagation distance R in the ice medium cice The sum of is expressed as (as shown in Formula 4):

[0016] Rc=R cair +R cice *η ice (Formula 4)

[0017] When refraction occurs on the ice surface, the incident angle θ i and the refraction angle θ o It also satisfies Snell's law of refraction, where η ice is the refractive index of the ice medium.

[0018] When calculating the signal propagation distance under the high curvature ice surface model, the propagation distance in the air and ice medium will change due to the curvature of the ice surface. Figure 4 As shown, at this time, for the propagation distance R in the air cair , propagation distance R in ice medium cice , angle of incidence θ i , slope angle θ, refraction angle θ o , the distance Δ of the refraction point from the ground, the aircraft flight altitude h and the radar detection depth d should satisfy the following relationship (as shown in Formula 5):

[0019]

[0020] The propagation distance R in a high curvature ice environment is obtained by formula 5 c .

[0021] Propagation distance difference generation module: This module only has a part of calculating the distance difference. The flow chart is as follows Figure 5As shown, the distance difference ΔR (as shown in Formula 6) is obtained by subtracting the flat ice surface signal propagation distance obtained in the above module from the high curvature ice surface signal propagation distance:

[0022] ΔR=R c -R f (Formula 6)

[0023] Thus, the difference in signal propagation distances under the two ground models is obtained.

[0024] Secondary phase error calculation module: Since the secondary phase error calculation method is required when analyzing the sub-aperture length, this module uses the propagation distance difference obtained in the previous module to calculate the size of the secondary phase error. The flowchart of this module is as follows: Figure 6 As shown in , it includes two parts: calculating the Taylor expansion and calculating the quadratic phase error. We define the azimuth distance from the sub-ice point target to the sky and earth point directly below the radar as X, as Figure 4 As shown in a The relationship is as follows (as shown in Formula 7):

[0025] L a =2X (Formula 7)

[0026] The obtained propagation distance difference ΔR is Taylor expanded to be a Taylor expansion about X (as shown in Formula 8):

[0027] ΔR(X)≈k1X+k2X 2 +k3X 3 (Formula 8)

[0028] Where k1, k2, and k3 represent the constant coefficients of the first, second, and third terms in the Taylor expansion, respectively. The values ​​of k1 and k3 are almost 0, so ΔR is mainly determined by k2. Therefore, the quadratic phase error can be expressed as (as shown in Formula 9):

[0029]

[0030] At this point, the quadratic phase error expression of the distance difference is obtained, where f c represents the center frequency of the radar, c represents the speed of light, and k2 is the second term of the Taylor expansion constant.

[0031] Sub-aperture length calculation module: This module uses the quadratic phase error obtained in the quadratic phase error calculation module as input, and analyzes the quadratic phase error obtained by formula 9 through the quadratic phase error discrimination method. This module only has the sub-aperture length calculation part. The module flow chart is as follows: Figure 7 As shown, the magnitude of the quadratic phase error is set to be less than π / 4 (as shown in Formula 10):

[0032]

[0033] When this condition is met, the distance error caused by the quadratic phase error can be ignored. That is to say, the distance X obtained by the quadratic phase error at this time satisfies the requirement of ignoring the propagation distance difference between the flat ice surface and the high curvature ice surface. According to formula 7, the sub-aperture length L that meets the requirements at this time can be calculated. a .

[0034] Ice radar data imaging module: This module uses the subaperture length L obtained in the subaperture length calculation module a As input, imaging is performed through the fast back-projection method, including two parts: sub-image calculation module and final image calculation module. The algorithm flow chart is as follows: Figure 8 As shown, the radar signal S after pulse compression out It can be expressed as (as shown in formula 11):

[0035]

[0036] Where a0 represents the original signal envelope, t f For fast time, t s is the slow time, f c represents the center frequency of the radar, and R represents the distance from the target to the radar.

[0037] The sub-image calculation process can be expressed by the following formula: out exist arrive The sub-image is obtained by summing within the range of (as shown in Formula 12):

[0038]

[0039] Where Lsub and L obtained by formula 7 a The same value has the same meaning, both representing the subaperture length, L f sub represents the fastest time acquired by the radar when generating the sub-image, t s sub represents the slow time of radar acquisition when generating sub-images, t s sub′ Indicates that in the process of generating sub-images, when t f sub When the value remains unchanged, arrive All t within the range s sub The value of S out is the radar signal after pulse compression calculated by formula 11.

[0040] The final image calculation process can be expressed by the following formula: sub exist arrive The final image is obtained by summing within the range of (as shown in Formula 13):

[0041]

[0042] Where t f For fast time, t s is the slow time, t s ′ means that in the process of generating the image, when t f When the value remains unchanged, arrive All t in the range s The value of L represents the full aperture length, S sub is the sub-image result calculated by formula 7, is the input signal of formula 12, f c Indicates the center frequency

[0043] The value of L is given by Equation 14:

[0044] L=2*(h*tan(θ i )+d*tan(θ o )) (Formula 14)

[0045] In the formula, h is the distance from the radar to the ice surface, d is the maximum depth of radar detection, and θ is i represents the incident angle at the air-ice interface, θ o represents the refraction angle at the air-ice interface.

[0046] So far, this method has obtained the final imaging result, the imaging result is as follows Fig. 9 shown.

[0047] Beneficial Effects

[0048] The present invention proposes an ice radar data imaging method based on ice surface curvature. The method can obtain more accurate ice radar data imaging results and improve the imaging quality in high-curvature ice surface environments. BRIEF DESCRIPTION OF THE DRAWINGS

[0049] Figure 1 The present invention is implemented in an overall flow chart.

[0050] Figure 2 This is a flow chart of the flat ice surface propagation distance generation module of the present invention.

[0051] Figure 3 This is a flow chart of the high-curvature ice surface propagation distance generation module of the present invention.

[0052] Figure 4 This is a signal propagation model diagram of the present invention.

[0053] Figure 5 This is a flow chart of the propagation distance difference generation module of the present invention.

[0054] Figure 6 This is a flow chart of the secondary phase error calculation module of the present invention.

[0055] Figure 7 This is a flow chart of the sub-aperture length calculation module of the present invention.

[0056] Figure 8 This is a flow chart of the ice radar data imaging module of the present invention.

[0057] Fig. 9 Schematic diagram of the final imaging result. DETAILED DESCRIPTION

[0058] First, the original radar data and its ground elevation model (DEM) parameters are obtained. At this time, the following information can be obtained: the slope of the ice surface θ on the flight path of the airborne radar, the aircraft flight altitude h, the maximum depth measured ice thickness d, and the azimuth distance X and the incident angle θ i , use the above parameters to substitute into formula 5 for calculation, and calculate the distance Rc from the target to the radar under the inclined ice surface model. Also use the above parameters, assuming that the ice surface slope θ is 0, and use formula 5 again to calculate the distance R from the target to the radar under the current position flat ice surface model f , the R obtained in the above steps c and R f Substitute into formula 6 to calculate the distance difference ΔR, and expand ΔR into a Taylor expansion with respect to X. In the obtained Taylor expansion, the values ​​of the constant terms k1 and k3 are almost 0, so the size of QPE is mainly determined by k2, which is the quadratic term.

[0059] Next, the radar center frequency f c , the speed of light c and the Taylor expansion constant k2 just obtained are substituted into Formula 9 to calculate the QPE size at this time, and then use Formula 10 to obtain At this time, the size of X can be inversely solved. Here, the X obtained is the required azimuth distance. Substitute X into formula 7 to calculate the required sub-aperture length L a .

[0060] Using the resulting subaperture length L a Substituting into formula 12, Lsub is equivalent to L a , both represent the subaperture length, that is, the input signal is summed and integrated along the azimuth direction within a subaperture length interval, where t sIndicates slow time, i.e. azimuth time, t f represents fast time, i.e., distance time, is the center frequency of the radar system, and the input Sout is the original signal after pulse compression. The sub-image result can be obtained by the above calculation.

[0061] Substitute the obtained sub-image result into formula 13. Similar to formula 12, the above result is summed and integrated in the azimuth direction within the full aperture L, so that the sub-images can be synthesized into the final imaging result. In the formula, t s Indicates slow time, i.e. azimuth time, t f represents fast time, i.e. distance time, f c is the center frequency of the radar system, the input S sub The sub-image result obtained by formula 12 can be obtained after calculation by formula 13. Fig. 9 The final image result is shown.

[0062] It should be understood that although this specification is described based on implementation methods, not every implementation method contains only one independent technical solution. This narrative method of the specification is only for the sake of clarity. Those skilled in the art should regard the specification as a whole. The technical solutions in each implementation method may also be appropriately combined and implemented according to the understanding of those skilled in the art.

[0063] The series of detailed descriptions listed above are only specific descriptions of feasible implementation modes of the present invention. They are not used to limit the protection scope of the present invention. All equivalent implementation modes or changes that do not deviate from the spirit of the invention should be included in the protection scope of the present invention.

Claims

1. An ice radar imaging method based on ice surface curvature, characterized in that: The following steps are involved: Step 1: Calculation of propagation distance difference: (1) According to the propagation distance of the radar to the point target under the flat ice model and the propagation distance of the radar to the point target under the high curvature ice model, the flat ice model propagation distance R is calculated respectively. f And the propagation distance R of the high curvature ice surface model c ; (2) Calculate the propagation distance R of the flat ice surface model f Propagation distance R with high curvature ice surface model c The distance difference ΔR is obtained by expanding the distance difference ΔR into a Taylor expansion of the azimuth distance X from the radar to the point target; Step 2: Calculation of the secondary phase error QPE: (1) Calculate the quadratic phase error QPE according to the Taylor expansion of the azimuth distance X; (2) Let the quadratic phase error QPE be Step 3: Subaperture length calculation: (1) The subaperture length L is obtained by inversely solving the quadratic phase error QPE a ; Step 4: Ice radar data imaging: (1) Through the sub-aperture length L a Get a sub-image; (2) The obtained sub-images are synthesized into the final imaging result.

2. The ice radar imaging method based on ice surface curvature according to claim 1, characterized in that: The step 1 of calculating the propagation distance difference comprises: The propagation distance R of the flat ice surface model is obtained by formula 1 and formula 2 f (X) and the propagation distance R of the high curvature ice surface model c (X): R f (X)=R fair (X)+R fice (X)*η ice (official1) R c (X)=R cair (X)+R cice (X)*η ice (official2) In formula 1 and formula 2, X represents the azimuth distance from the radar to the point target, R fair (X) represents the propagation distance in the air under the flat ice surface model, R fice (X) represents the propagation distance in the ice under the flat ice surface model, R cair (X) represents the propagation distance in the air under the high curvature ice surface model, R cice (X) represents the propagation distance in the ice under the high curvature ice surface model, η ice represents the refractive index of the air-to-ice interface, which is a known quantity; Propagate the flat ice surface model by a distance R f (X) and the propagation distance R of the high curvature ice surface model c (X) Subtract and calculate the distance difference ΔR of the propagation distance under the two ice surface models; ΔR=R f (X)- R c (X) (Formula 3) The distance difference ΔR is expanded into a Taylor expansion about the azimuth distance X, which is expressed by formula 4: ΔR(X) ≈ k1X + k2X 2 +k3X 3 (Formula 4) In formula 4, k1, k2, and k3 represent X, X in Taylor expansion respectively. 2 , X 3 The coefficients of k1 and k3 are almost 0, so ΔR is mainly determined by k2; The second step of calculating the secondary phase error QPE comprises: The quadratic phase error QPE is calculated according to Formula 5: In formula 5, f c represents the center frequency, c represents the speed of light, and k2 is the X in the Taylor expansion. 2 The coefficient of The third step of calculating the subaperture length comprises: Let the quadratic phase error The azimuth distance X is obtained by inverse solution and substituted into formula 6 to obtain the subaperture length L a , L a =2X (Formula 6) The fourth step of ice radar data imaging comprises: The sub-aperture length La when generating the sub-image is set to Lsub. arrive Within the range of all radar signals S after pulse compression out Perform sum operation to obtain sub-image S sub , the process is expressed by formula 7: In formula 7, t f sub represents the fastest time acquired by the radar when generating the sub-image, t s sub represents the slow time of radar acquisition when generating sub-images, t s sub′ Indicates that in the process of generating sub-images, when t f sub When the value remains unchanged, arrive All t in the range s sub The value of S out represents the input signal obtained from the radar, f c Indicates the center frequency of the radar when it is working; When generating the final image, the sub-image result S sub Substitute into formula 8 and convert the sub-image result S sub exist arrive The final imaging result S is obtained by summing within the range final , where L represents the full aperture length, In formula 8, t f For fast time, t s is the slow time, t s ′ means that in the process of generating the image, when t f When the value remains unchanged, arrive All t in the range s The value of Ssub is the sub-image result calculated by formula 7, and is also the input signal of formula 8. c represents the center frequency, The value of L is given by Equation 9: L = 2 * (h * tan(θ i ) + d * tan(θ o )) (Equation 9) In formula 9, h represents the distance from the radar to the ice surface, d represents the maximum depth of radar detection, and θ i represents the incident angle at the air-ice interface, θ o represents the refraction angle at the air-ice interface.

Images