Active and passive joint inversion method for ice cover internal temperature profile based on Bayesian estimation and product

Abstract

The invention discloses a Bayesian estimation-based active and passive joint inversion method and a Bayesian estimation-based active and passive joint inversion product for an ice cover internal temperature profile, and effectively reduces the uncertainty of ice cover radiation interpretation by introducing snow accumulation layer density fluctuation characteristics and ice cover internal structure information estimated by ice penetrating radar echoes. On the basis, optimizing and updating initial values and probability density distribution of to-be-estimated parameters of an ice cover internal temperature profile by taking ice cover dielectric absorption attenuation characteristics estimated based on an ice penetrating radar echo profile as constraint conditions; in the process of introducing prior information into ice cover internal temperature profile inversion based on broadband radiation brightness temperature observation, constraint information provided actively and passively is fully combined, active and passive joint inversion of the ice cover internal temperature profile is realized, the inversion performance of the ice cover internal temperature profile is effectively improved, and the inversion efficiency of the ice cover internal temperature profile is improved. The dimensionality and uncertainty of the ice cover internal temperature profile inversion problem are reduced, and the reliability and robustness of the ice cover internal temperature profile inversion result are improved.

Description

Technical Field

[0001] The present invention relates to the field of microwave remote sensing, and in particular to an active and passive joint inversion method and product of an ice sheet internal temperature profile based on Bayesian estimation. Background Art

[0002] In the polar ice sheet region, the internal temperature profile of the ice sheet refers to the physical temperature that changes with depth, reflecting the vertical distribution characteristics of the internal temperature of the ice sheet. It is affected by the surface accumulation rate, surface temperature, geothermal flux, ice thickness, heat exchange and stress in the ice, and is the main factor determining the dynamic and thermal processes inside the ice sheet. In order to meet the research needs of the evolution of the material balance of the ice sheet, it is necessary to effectively estimate and measure the internal temperature profile of the ice sheet. At present, the internal temperature profile of the ice sheet is mainly measured on-site through drilling technology. However, due to the limitations of the measurement environment and the cost of technical implementation, the field observation of the internal temperature profile of the ice sheet is extremely sparse and does not have the ability to provide large-scale spatial coverage. Therefore, the internal temperature profile of the ice sheet is still the most poorly understood internal state characteristic of the ice sheet in the study of ice sheet dynamics and thermodynamics. It is of great significance to develop microwave remote sensing detection methods for the internal temperature profile of the ice sheet.

[0003] Previously, a method of using passive remote sensing to detect the temperature profile inside the ice sheet has been proposed. Based on the fact that the radiation brightness temperature of different frequency channels in the UHF band is sensitive to the relative changes in temperature at different depths of the ice sheet, the temperature profile inside the ice sheet can be inverted based on the ultra-wideband microwave radiation brightness temperature in the UHF band. However, after analyzing the acquired data, it was found that in addition to the temperature profile inside the ice sheet, the radiation brightness temperature characteristics of the ice sheet are also affected by a series of factors including the internal layered structure of the ice sheet and the fluctuation characteristics of the snow layer density. In the absence of effective information on the internal layered structure of the ice sheet and the fluctuation characteristics of the snow layer density, since a single passive detection channel cannot effectively constrain the internal structure characteristics of the ice sheet, the interference caused by the internal layered structure of the ice sheet and the fluctuation characteristics of the snow layer density on the interpretation of the radiation brightness temperature will significantly limit the ability to use passive microwave radiation brightness temperature to invert the temperature profile inside the ice sheet.

[0004] In addition, a method has been proposed before to estimate the average temperature inside the ice sheet based on the dielectric absorption attenuation rate of the ice sheet estimated by ice-penetrating radar echoes. However, this method uses many simplified assumptions when estimating the dielectric absorption attenuation rate of the ice sheet using ice-penetrating radar echoes, such as ignoring the vertical changes in the reflectivity of the reflection interface inside the ice sheet. Therefore, only the average temperature inside the ice sheet is estimated at present, and it cannot provide information on the vertical changes in the temperature of the ice sheet. Summary of the invention

[0005] The purpose of the present invention is to overcome the defects of the prior art and propose an active and passive joint inversion method and product of the internal temperature profile of the ice sheet based on Bayesian estimation.

[0006] In view of this, the present invention proposes an active and passive joint inversion method for the internal temperature profile of the ice sheet based on Bayesian estimation, comprising:

[0007] Step 1) Obtain broadband radiation brightness temperature and ice penetrating radar echo profile observations of the ice sheet area to be studied;

[0008] Step 2) extracting the internal structure and thickness parameters of the ice sheet based on the ice-penetrating radar echo profile observations, and estimating the fluctuation characteristics of the snow layer density;

[0009] Step 3) determining the ice sheet internal temperature profile model, and determining the ice sheet internal temperature profile parameters to be estimated and redundant parameters to be estimated, and performing inversion problem description;

[0010] Step 4) Discretize the parameter space to be estimated, combine the snow layer density fluctuation characteristics, the internal structure of the ice cover and the ice thickness parameters, perform active and passive forward simulations on the discrete parameter space to be estimated, and generate a simulated data set of broadband radiation brightness temperature and ice penetrating radar echo profile;

[0011] Step 5) estimating the dielectric absorption attenuation value of the ice sheet according to the ice penetrating radar echo profile observation value, and optimizing the initial value and probability density distribution of the parameters to be estimated of the temperature profile inside the ice sheet based on the dielectric absorption attenuation value of the ice sheet estimated by the ice penetrating radar echo under the Bayesian estimation framework;

[0012] Step 6) Based on the broadband radiation brightness temperature under the Bayesian estimation framework, the estimated parameters and redundant estimated parameters of the internal temperature profile of the ice sheet are inverted. During the inversion, the initial values ​​and probability density distribution of the estimated parameters of the internal temperature profile of the ice sheet optimized in step 5) are used as prior information. According to the parameter inversion results and the internal temperature profile model of the ice sheet in step 3), the internal temperature profile of the ice sheet in the target area is obtained.

[0013] Preferably, the step 1) comprises:

[0014] Along the ice sheet area to be studied, the intervals are divided by active ice-penetrating radar detection and passive broadband radiometer detection segments, and the broadband radiation brightness temperature observation values ​​and ice-penetrating radar echo intensity time series observation values ​​of each interval are obtained. The electromagnetic wave propagation speed v in the ice sheet is determined according to the following formula:

[0015]

[0016] Where c is the speed of light, v is in m / s, and ε′ is the real part of the average relative permittivity of the ice sheet;

[0017] According to v, the time series of ice-penetrating radar echo intensity is converted into the depth distribution of ice-penetrating radar echo intensity, so as to obtain the ice-penetrating radar echo profile observation values ​​in each interval.

[0018] Preferably, the step 2) comprises:

[0019] Step 2-1) extract the reflection target position depth in the ice penetrating radar echo profile observation value of each interval in the ice sheet area to be studied, divide the internal structure of the ice sheet, and calculate the ice sheet thickness in each interval according to the reflection target position depth of the ice bottom and the reflection position of the ice surface;

[0020] Step 2-2) Based on the empirical model of snow layer density fluctuation characteristics and combined with the parameters of the ice-penetrating radar detection system, the ice-penetrating radar echo profile is forward simulated. By adjusting the snow layer density fluctuation characteristic parameters, the ice-penetrating radar echo profile observations of the snow layer area in each interval are fitted, and the empirical model parameters of snow layer density fluctuation characteristics that achieve the best fit to the ice-penetrating radar echo profile observations in each interval are determined.

[0021] Preferably, the step 3) comprises:

[0022] Step 3-1) Based on the borehole measurement data of the internal temperature profile of the ice sheet, establish the internal temperature profile model T(z) of the ice sheet through polynomial fitting, moment method or by solving the heat conduction equation:

[0023] T(z)=g(x t ,z)

[0024] x t =[x t,1 ,x t,2 ,…,x t,p ] T

[0025] Among them, z is the position depth, and the corresponding parameter to be estimated is x t,1 ,x t,2 ,…,x t,p Construct the temperature parameter vector x t , p represents the total number of parameters to be estimated, and the superscript “T” represents transposition;

[0026] Step 3-2) From the redundant parameters to be estimated x n,i Constructing the redundant parameter vector x n :

[0027] x n =[x n,1 ,x n,2 ,…,x n,q ] T

[0028] Where q represents the total number of redundant parameters;

[0029] Step 3-3) The inversion problem of the temperature profile inside the ice sheet is expressed as:

[0030]

[0031] y=f(x)+em

[0032] in, For broadband radiation brightness temperature observation, with N frequency points y is the observation vector, which is composed of broadband radiation brightness temperature observations and has a dimension of N. f(x) represents the simulation result corresponding to the radiation forward model under the model parameter vector x, and em is the observation noise vector;

[0033] Then the likelihood function L(x) is:

[0034]

[0035] Where C is the noise covariance matrix.

[0036] Preferably, the step 4) comprises:

[0037] Step 4-1) According to the effective value range of the parameter to be estimated, the temperature parameter vector x t and the redundant parameter vector x n The parameter space of the model parameter vector x is discretized to generate the corresponding parameter grid.

[0038] Step 4-2) combines the snow layer density fluctuation characteristics, ice cover internal structure and ice thickness parameters obtained in step 2) to perform forward simulation of the ice cover broadband radiation brightness temperature and ice-penetrating radar echo profile at discrete grid points in the model parameter space to generate the ice cover broadband radiation brightness temperature and ice-penetrating radar echo profile grid; and perform interpolation and refinement to generate a simulated data set of broadband radiation brightness temperature and ice-penetrating radar echo profile.

[0039] Preferably, the step 5) comprises:

[0040] Step 5-1) The estimated parameters of the internal temperature profile of the ice sheet and the corresponding temperature parameter vector x determined in step 3) t , combined with the discrete temperature parameter vector x in step 4) t The simulated data set of ice-penetrating radar echo profile corresponding to the parameter grid is used to obtain the temperature parameter vector x t Simulated values ​​of one-way ice sheet dielectric absorption attenuation corresponding to all parameter grid points in parameter space and the simulated value of the total dielectric absorption attenuation of the one-way ice sheet

[0041] Step 5-2) The dielectric absorption attenuation value N of the ice sheet estimated based on the ice penetrating radar echo profile observations in each interval of the target area is a and the total dielectric absorption attenuation of the ice sheet L a,EFZ , estimate the temperature parameter vector x in each intervalt All parameter grid points in parameter space The posterior probability

[0042] Step 5-3) The posterior probability In the temperature parameter vector x t All parameter grid points in parameter space Calculate on to get the temperature parameter vector x t Probability density distribution on parameter space;

[0043] Step 5-4) Calculate the temperature parameter vector x using the minimum mean square error or maximum a posteriori probability method t Estimated value of

[0044] Step 5-5) Used to update the temperature parameter vector x in each interval t All parameter grid points in parameter space The corresponding probability density is the temperature parameter vector x t Estimated value of It is used to update the initial values ​​of the estimated parameters of the internal temperature profile of the ice sheet in each interval. The two are substituted as prior information into the subsequent Bayesian estimation process of the temperature profile parameters inside the ice sheet based on broadband radiation brightness temperature, and combined with the constraint information provided by active and passive methods, the active and passive joint inversion of the temperature profile parameters inside the ice sheet is realized.

[0045] Preferably, the posterior probability of step 5-2) for:

[0046]

[0047] in, is the temperature parameter grid point The corresponding prior probability density distribution;

[0048]

[0049] in, and They correspond to the estimated ice sheet dielectric absorption attenuation value N based on the ice penetrating radar echo profile observation value. a and L a,EFZ The estimated standard deviation of

[0050] Step 5-4) t Estimated value of The calculation is based on the minimum mean square error method:

[0051]

[0052] x t Estimated value of The maximum a posteriori probability method is calculated as:

[0053]

[0054] Preferably, the step 6) comprises:

[0055] Step 6-1) Based on the redundant parameters to be estimated and the corresponding redundant parameter vector x determined in step 3) n , combined with the discrete temperature parameter vector x in step 4) t The broadband radiation brightness temperature simulation data set corresponding to the parameter grid, as well as the broadband radiation brightness temperature observation vector y in each interval of the target area and the temperature parameter vector x obtained in step 5) t Estimated value of Calculate the redundant parameter vector x on each interval n All parameter grid points in parameter space The posterior probability

[0056]

[0057] Step 6-2) Calculate the temperature parameter vector x in each interval t All parameter grid points in parameter space The posterior probability

[0058]

[0059] Posterior probability In the temperature parameter vector x t All parameter grid points in parameter space The calculation result is the temperature parameter vector x after active and passive joint constraints in each interval. t Probability density distribution on parameter space;

[0060] Step 6-3) Calculate the temperature parameter vector x using the minimum mean square error method t Estimated value of

[0061]

[0062] Or use the maximum a posteriori probability method to calculate the temperature parameter vector x t Estimated value of

[0063]

[0064] The corresponding estimated variance is for:

[0065]

[0066] Step 6-4) Calculate the estimated results of the internal temperature profile of the ice sheet in each interval, according to the temperature parameter vector x in each interval t All parameter grid points in parameter space The corresponding posterior probability The minimum mean square error method is used to calculate the estimated value of the ice sheet temperature at any depth z, thereby obtaining the inversion result of the internal temperature profile of the ice sheet in the target area; among them, the minimum mean square error Calculated as:

[0067]

[0068] Corresponding estimated variance for:

[0069]

[0070] in, is the temperature parameter vector x in each interval t All parameter grid points in the parameter space, j represents the jth grid point.

[0071] In another aspect, the present invention provides a computer program product, comprising a computer program, which implements the steps of the above method when executed by a processor.

[0072] Compared with the prior art, the advantages of the present invention are:

[0073] 1. The present invention provides an active and passive joint inversion algorithm for the internal temperature profile of the ice sheet based on Bayesian estimation, which changes the current situation that the internal temperature profile of the ice sheet can only be detected through on-site drilling technology measurement or single passive or active remote sensing means.

[0074] 2. The present invention can constrain the internal structure of the ice sheet and the density fluctuation characteristics of the snow layer through the active ice-penetrating radar echo profile, reduce the interference of the internal structure of the ice sheet and the density fluctuation characteristics of the snow layer in the interpretation of the radiation brightness temperature, and at the same time, based on the correlation between the dielectric absorption attenuation characteristics of the ice sheet estimated by the ice-penetrating radar echo profile and the internal temperature of the ice sheet, the passive broadband radiation brightness temperature is combined to perform active and passive joint constraints on the internal temperature profile of the ice sheet, effectively reducing the dimension and uncertainty of the inversion problem of the internal temperature profile of the ice sheet, and improving the effectiveness and reliability of the inversion of the internal temperature profile of the ice sheet.

[0075] 3. The present invention proposes to optimize the initial values ​​and probability density distribution of the parameters to be estimated of the temperature profile inside the ice sheet based on the dielectric absorption attenuation value of the ice sheet estimated by the ice-penetrating radar echo under the Bayesian estimation framework, and substitute the optimization results as prior information into the inversion process based on broadband radiation brightness temperature observations, so as to realize the active and passive joint inversion method of the temperature profile inside the ice sheet by combining the constraint information provided by the active and passive methods. This joint method can reduce the impact that may be caused by the simplified assumptions used in the process of estimating the dielectric absorption attenuation of the ice sheet based on the ice-penetrating radar echo profile. Compared with the joint inversion method of using the dielectric absorption attenuation value of the ice sheet estimated based on the ice-penetrating radar echo profile and the broadband radiation brightness temperature observation value as observation quantities directly for the temperature profile inside the ice sheet, the joint inversion algorithm proposed by the present invention is more robust. BRIEF DESCRIPTION OF THE DRAWINGS

[0076] Figure 1 It is a schematic diagram of the selected detection area of ​​the Greenland ice sheet and the corresponding active and passive detection flight tracks;

[0077] Figure 2(a) is the radiation brightness temperature observation data corresponding to the 12 frequency channels of the ultra-wideband software-defined microwave radiometer on the selected flight track, and Figure 2(b) is the average value in each interval;

[0078] Figure 3 It is the radiation brightness temperature observation data corresponding to the interval where the NGRIP, NEEM and Camp Century boreholes are located;

[0079] Figure 4(a) is an echo profile of the multi-channel coherent bathymetric radar and accumulation rate radar in the interval where the NGRIP borehole is located, Figure 4(b) is an echo profile of the multi-channel coherent bathymetric radar and accumulation rate radar in the interval where the NEEM borehole is located, and Figure 4(c) is an echo profile of the multi-channel coherent bathymetric radar and accumulation rate radar in the interval where the Camp Century borehole is located;

[0080] Figure 5 is the thickness of the ice cover along the flight track in the selected target area;

[0081] Figure 6 is a graph showing the single density sample generation results and the corresponding density fluctuation standard deviation distribution results corresponding to the estimated values ​​of the snow layer density fluctuation standard deviation parameter in the interval [1, 41, 75, 118]; Figure 6 (a) is interval 1, Figure 6 (b) is interval 41, Figure 6 (c) is interval 75, and Figure 6 (d) is interval 118;

[0082] Figure 7 shows the initial prior probability density distribution of the parameters to be estimated for the temperature profile inside the ice sheet in the intervals where the NGRIP, NEEM and Camp Century boreholes are located, as well as the probability density distribution after optimization using the ice sheet dielectric absorption attenuation value estimated based on the multi-channel coherent bathymetric radar echo profile; the horizontal axis of Figure 7(a) is the annual average surface temperature, the horizontal axis of Figure 7(b) is the annual average accumulation rate, and the horizontal axis of Figure 7(c) is the geothermal flux;

[0083] Figure 8 shows the inversion results of the parameters to be estimated for the temperature profile inside the ice sheet on the flight track from Camp Century to NGRIP, where Figure 8(a) is the annual mean surface temperature, Figure 8(b) is the annual mean accumulation rate, and Figure 8(c) is the geothermal heat flux;

[0084] Figure 9 shows the inversion results of the internal temperature profile of the ice sheet on the flight track from Camp Century to NGRIP and the corresponding estimated standard deviation; Figure 9(a) is the active and passive combination, and Figure 9(b) is the passive alone;

[0085] Fig.10 is the average estimated standard deviation of ice sheet temperature at different depths along the flight track from Camp Century to NGRIP;

[0086] Figure 11 is a comparison of the inversion results of the internal temperature profile of the ice sheet in the interval where the NGRIP, NEEM and Camp Century boreholes are located and the borehole measurement results of the internal temperature profile of the ice sheet fitted by the Robin model, where Figure 11(a) is NGRIP, Figure 11(b) is NEEM, and Figure 11(c) is Camp Century;

[0087] Fig.12 It is a flow chart of the active and passive joint inversion method of the internal temperature profile of the ice sheet based on Bayesian estimation of the present invention. DETAILED DESCRIPTION

[0088] The purpose of the present invention is to overcome the defect that the actual measurement data of the internal temperature profile of the ice sheet is very scarce due to the extremely sparse distribution of borehole measurements, which cannot meet the research needs of the thermodynamic process and material balance evolution of the ice sheet, overcome the interference of the internal structure of the ice sheet and the density fluctuation characteristics of the snow layer in the previous method of inverting the internal temperature profile of the ice sheet using passive remote sensing ultra-wideband radiation brightness temperature, and overcome the limitation that the dielectric absorption attenuation rate of the ice sheet estimated by ice penetrating radar echo can only estimate the average temperature of the ice sheet, and cannot effectively constrain the vertical change characteristics of the internal temperature of the ice sheet, and provide an active and passive joint inversion algorithm for the internal temperature profile of the ice sheet based on Bayesian estimation. On the basis of the sensitivity of the passive broadband radiation brightness temperature observation channel to the internal temperature profile of the ice sheet, the ice penetrating radar echo profile characteristics are further combined, on the one hand, the snow layer density fluctuation characteristic information and the internal structure information of the ice sheet constrained by the reflection characteristics of the ice penetrating radar echo profile are introduced to reduce the uncertainty of radiation interpretation; on the other hand, according to the correlation between the dielectric absorption attenuation of the ice sheet and the internal temperature of the ice sheet, the ice penetrating radar echo profile attenuation characteristics are used to further constrain the internal temperature profile of the ice sheet. Based on Bayesian estimation, the present invention proposes an algorithm for active and passive joint inversion of the temperature profile inside the ice sheet using ice-penetrating radar echoes and broadband radiation brightness temperature. The method proposed in the present invention can effectively reduce or eliminate the interference of snow layer density fluctuations in the inversion of the temperature profile inside the ice sheet under general actual ice sheet conditions. At the same time, the joint constraint information provided by the active and passive detection channels can be used to improve the effectiveness and reliability of the inversion of the temperature profile inside the ice sheet, and effectively reduce the dimensionality and uncertainty of the inversion problem of the temperature profile inside the ice sheet.

[0089] The present invention provides an active and passive joint inversion method for the internal temperature profile of an ice sheet based on Bayesian estimation, comprising:

[0090] Step 1) Select the target area and obtain broadband radiation brightness temperature and ice penetrating radar echo profile observation values;

[0091] Step 2) extracting the internal structure and thickness parameters of the ice sheet according to the ice-penetrating radar echo profile observations in step 1), and estimating the fluctuation characteristics of the snow layer density;

[0092] Step 3) determining the ice sheet internal temperature profile model, and determining the ice sheet internal temperature profile parameters to be estimated and redundant parameters to be estimated, and performing inversion problem description;

[0093] Step 4) discretizes the parameter space to be estimated determined in step 3), combines the snow layer density fluctuation characteristics, ice cover internal structure and ice thickness parameters obtained in step 2), performs active and passive forward simulation on the discrete parameter space to be estimated, and generates a simulated data set of broadband radiation brightness temperature and ice penetrating radar echo profile;

[0094] Step 5) estimating the dielectric absorption attenuation value of the ice sheet according to the ice-penetrating radar echo profile observation value in step 1), and optimizing the initial value and probability density distribution of the parameters to be estimated of the temperature profile inside the ice sheet based on the dielectric absorption attenuation value of the ice sheet estimated by the ice-penetrating radar echo under the Bayesian estimation framework;

[0095] Step 6) Under the Bayesian estimation framework, according to the inversion problem description in step 3), the parameters to be estimated and redundant parameters to be estimated of the internal temperature profile of the ice sheet are inverted based on the broadband radiation brightness temperature, wherein the initial values ​​and probability density distribution of the parameters to be estimated of the internal temperature profile of the ice sheet obtained by optimization in step 5) are used as prior information. The internal temperature profile of the ice sheet in the target area is obtained based on the parameter inversion results and the internal temperature profile model of the ice sheet in step 3).

[0096] The method steps are systematically described below.

[0097] In step 1), the target area is selected to obtain broadband radiation brightness temperature and ice penetrating radar echo profile observation values, which specifically include:

[0098] Step 1-1) selecting a target ice sheet region for which the internal temperature profile of the ice sheet is to be estimated;

[0099] Step 1-2) The active ice-penetrating radar detection and passive broadband radiometer detection segments in the target area are divided into intervals to obtain broadband radiation brightness temperature observations and ice-penetrating radar echo intensity time series observations in each interval. Then, the electromagnetic wave propagation speed in the ice sheet is determined based on the real part of the average relative dielectric constant ε′ of the ice sheet:

[0100]

[0101] Among them, c is the speed of light, and v is the propagation speed of electromagnetic waves in the ice sheet, with the unit of m / s.

[0102] According to the propagation speed of electromagnetic waves in the ice sheet, the time series of ice-penetrating radar echo intensity is converted into the depth distribution of ice-penetrating radar echo intensity, so as to obtain the observation values ​​of ice-penetrating radar echo profile in each interval.

[0103] In step 2), the internal structure and thickness parameters of the ice sheet are extracted according to the ice-penetrating radar echo profile observations in step 1), and the fluctuation characteristics of the snow layer density are estimated, which specifically includes:

[0104] Step 2-1) Extract the reflection target position depth in the ice-penetrating radar echo profile observation value of each interval in the target area, and divide the internal structure of the ice cover in the target area. Calculate the ice cover thickness in each interval based on the reflection target position depth at the ice bottom and the reflection position on the ice surface.

[0105] Step 2-2) Based on the empirical model of snow layer density fluctuation characteristics and combined with the parameters of the ice-penetrating radar detection system, the ice-penetrating radar echo profile is forward simulated. By adjusting the snow layer density fluctuation characteristic parameters, the ice-penetrating radar echo profile observations of the snow layer area in each interval are fitted, and the empirical model parameters of snow layer density fluctuation characteristics that achieve the best fit to the ice-penetrating radar echo profile observations in each interval are determined to describe the snow layer density fluctuation characteristics in each interval of the target area.

[0106] In step 3), the ice sheet internal temperature profile model is determined, and the parameters to be estimated and the redundant parameters to be estimated of the ice sheet internal temperature profile are determined, specifically including:

[0107] Step 3-1) Based on the borehole measurement data of the internal temperature profile of the ice sheet, the internal temperature profile of the ice sheet can be modeled by polynomial fitting, moment method, etc., or by solving the heat conduction equation to establish the internal temperature profile model T(z) of the ice sheet, where z is the location depth and the corresponding parameter to be estimated x is t,i (i=1,…,p) constitutes the temperature parameter vector x t , the superscript “T” indicates transpose.

[0108] T(z)=g(x t ,z)

[0109] x t =[x t,1 ,x t,2 ,…,x t,p ] T

[0110] Taking the solution of the one-dimensional heat conduction equation as an example, when the ice sheet is in a thermally stable state and horizontal heat convection can be ignored, the internal temperature profile model T(z) of the ice sheet can be expressed as

[0111]

[0112] This model is also called the Robin model. Among them, erf is the error function, T s is the annual average surface temperature in K, G is the geothermal flux in W / m 2 , M is the average annual accumulation rate, in m / year, H is the ice thickness, in m, k c =2.7W / m / K is the thermal conductivity of ice, k d =45m 2 / year is the ice thermal diffusion rate. Figure 1 The fitting results of the internal temperature profiles of the Greenland and Antarctic ice sheets obtained by borehole measurements in some areas of the ice sheets are given using the Robin model. When the ice sheet thickness H can be obtained through existing products or active ice-penetrating radar data, the parameters to be estimated in this ice sheet internal temperature profile model include Ts , M, G, the corresponding temperature parameter vector is

[0113] x t =[T s ,M,G] T

[0114] Step 3-2) In the inversion process based on broadband radiation brightness temperature fitting, in addition to inverting the parameters to be estimated for the temperature profile inside the ice sheet, it is also necessary to determine other unknown parameters that affect the radiation brightness temperature fitting. These unknown parameters are not the parameters to be estimated for the inversion target of the temperature profile inside the ice sheet, but they need to be inverted together with the parameters to be estimated for the temperature profile inside the ice sheet during the inversion process. These unknown parameters are called redundant parameters to be estimated. Redundant parameters to be estimated x n,i (i=1,…,q) constitutes the redundant parameter vector x n .

[0115] x n =[x n,1 ,x n,2 ,…,x n,q ] T

[0116] Taking the snow layer two-scale density fluctuation correlation length parameters l1, l2 and small-scale density fluctuation proportional factor α as redundant parameters to be estimated, the redundant parameter vector x n for

[0117] x n =[l1,l2,α] T

[0118] In step 3), the inversion problem description specifically includes:

[0119] Step 3-3) In the inversion problem of the internal temperature profile of the ice sheet, the internal structure of the ice sheet and the ice sheet thickness parameters extracted by the ice-penetrating radar echo profile in step 2) are used as known parameters. In addition, the density fluctuation characteristics of the snow layer estimated by the ice-penetrating radar echo profile in step 2) are also used as known parameters; the ice-penetrating radar echo profile observations are used as constraints to optimize the initial values ​​and probability density distribution functions of the parameters to be estimated of the internal temperature profile of the ice sheet, which serve as prior information in the inversion problem; finally, the broadband radiation brightness temperature observations are used as observation vectors in the inversion problem to realize Bayesian estimation of the parameters to be estimated and redundant parameters to be estimated of the internal temperature profile of the ice sheet.

[0120] In the specific description, it is assumed that the broadband radiation brightness temperature observation With N frequency points The observation vector is defined as y, which is composed of broadband radiation brightness temperature observations and has a dimension of N. During the inversion process, the radiation brightness temperature simulation value T BThe temperature parameter vector x corresponding to the temperature profile of the interior of the ice sheet is determined by the parameters to be estimated and the redundant parameters to be estimated. t The redundant parameter vector x corresponding to the redundant estimated parameters n Together they constitute the model parameter vector x in the inversion problem of the temperature profile inside the ice sheet.

[0121]

[0122] T B =f(x)

[0123] Temperature parameter vector x t The parameter x to be estimated from the temperature profile inside the ice sheet t,i (i=1,…,n), the redundant parameter vector x n From the redundant parameter x to be estimated n,i (i=1,…,m). The observation noise vector em is independent of the model parameter vector x. f(x) represents the simulation result corresponding to the radiation forward model under the model parameter vector x. On this basis, the inversion problem of the temperature profile inside the ice sheet can be expressed as:

[0124]

[0125] y=f(x)+em

[0126] In this way, the likelihood function can be expressed as:

[0127]

[0128] Where C is the noise covariance matrix.

[0129] In step 4), the parameter space to be estimated determined in step 3) is discretized, and combined with the snow layer density fluctuation characteristics, ice cover internal structure and ice thickness parameters obtained in step 2), active and passive forward simulations are performed on the discrete parameter space to be estimated to generate a simulated data set of broadband radiation brightness temperature and ice penetrating radar echo profile. Specifically, it includes:

[0130] Step 4-1) According to the effective value range of the parameter to be estimated, the temperature parameter vector x t and the redundant parameter vector x n The parameter space of the model parameter vector x is discretized respectively to achieve the discretization of the parameter space and generate the corresponding parameter grid.

[0131] Step 4-2) Combined with the snow layer density fluctuation characteristics, ice cover internal structure and ice thickness parameters obtained in step 2), forward simulation of ice cover broadband radiation brightness temperature and ice penetrating radar echo profile is performed on discrete grid points in the model parameter space to generate ice cover broadband radiation brightness temperature and ice penetrating radar echo profile grids. The ice cover broadband radiation brightness temperature and ice penetrating radar echo profile grids can be interpolated and refined as needed to generate simulated data sets of broadband radiation brightness temperature and ice penetrating radar echo profiles.

[0132] In step 5), the dielectric absorption attenuation of the ice sheet is estimated according to the ice-penetrating radar echo profile observation value in step 1), specifically including:

[0133] Step 5-1) The ice-penetrating radar echo profile corresponds to the depth distribution of the ice-penetrating radar echo intensity from the ice sheet surface to the ice bottom, which is mainly determined by the radar system parameters, the reflectivity of the interface from the ice sheet surface to the ice bottom, the dielectric absorption attenuation of the ice sheet, and the spatial attenuation during the propagation of electromagnetic waves. Under the quasi-specular reflection approximation, the radar echo intensity can be expressed as

[0134]

[0135] in, is the intensity of the i-th reflected echo, P t is the transmission power, G t and G R are the transmitting and receiving gains respectively, λ is the radar wavelength, L a,i is the ice sheet one-way dielectric absorption attenuation experienced by the i-th reflected echo, Γ i is the total one-way transmittance of the ice sheet interface experienced by the i-th reflected echo, corresponding to the one-way reflection loss of the ice sheet surface and internal layered interface, L sys is the system loss, G p is the system gain, R i The reflectivity corresponding to the i-th reflected echo, h is the height of the observation platform, ε′ ice is the real part of the average relative permittivity of the ice sheet, z i The depth of the position corresponding to the i-th reflected echo.

[0136] Since the dielectric absorption attenuation of the ice sheet is a function of the ice sheet temperature, estimating the dielectric absorption attenuation of the ice sheet through the ice-penetrating radar echo profile can provide effective constraint information on the temperature profile inside the ice sheet.

[0137] At the platform height h, according to the reflection depth position z and the real part of the average relative dielectric constant ε′ of the ice sheet, ice , the spatial attenuation in the radar echo intensity can be compensated, and the intensity of the i-th reflected echo after compensation is

[0138]

[0139] Relative to surface echo intensity Reflection echo intensity Normalize and compensate for the radar system parameter P t , G t , G R ,λ,L sys and G p At this time, the normalized i-th reflection echo intensity L is attenuated only by the dielectric absorption of the ice sheet a,i , ice sheet reflection loss Γ i and the reflectivity R at the echo position i The decision can be expressed in the logarithmic domain as

[0140]

[0141] Step 5-2) The dielectric absorption attenuation rate of the ice sheet can be estimated by different algorithms such as the empirical method, the adaptive method and the echo gradient method. Taking the echo gradient method as an example, the dielectric absorption attenuation rate of the ice sheet above the non-reflective area of ​​the ice sheet can be estimated by normalizing the radar echo intensity. To avoid the influence of the reflection characteristics of the snow layer, the shallowest position z in the depth range of the gradient solution is min Not less than the thickness of the snow layer. In addition, considering that the echo intensity of the ice-penetrating radar near the bottom of the ice is weak and the signal-to-noise ratio is low, it may introduce deviations in the estimation of the dielectric absorption attenuation rate of the ice sheet. Therefore, the deepest position z in the depth range of the gradient solution is max The depth position of the upper boundary of the non-reflection zone above the ice bottom interface should not exceed that of the upper boundary of the non-reflection zone above the ice bottom interface. The depth gradient can be expressed as

[0142]

[0143] The one-way reflection loss Γ on the ice sheet surface and internal stratified interface is mainly generated by the shallow snow layer area. min and z max In the ice sheet region between Γ and Γ, Γ can be considered as a constant, so d[Γ] dB / dz can be ignored. In addition, the internal reflection below the snow layer of the ice sheet is mainly caused by the sudden change of conductivity or the change of the ice crystal axis orientation. The corresponding reflectivity R does not change significantly with depth, and its contribution to the depth gradient of radar echo intensity can be almost ignored. Therefore, d[R] dB The / dz term is also approximately equal to 0. In this way, the normalized corrected radar echo intensity The depth gradient is approximately equal to twice the depth gradient of the one-way dielectric absorption attenuation of the ice sheet. The depth gradient of the one-way dielectric absorption attenuation of the ice sheet is the one-way dielectric absorption attenuation rate of the ice sheet N. a .

[0144]

[0145] Ice sheet one-way dielectric absorption attenuation rate N a is a function of ice sheet temperature and ice chemistry

[0146]

[0147] Among them, N a The unit is dB / km, σ ∞ is the conductivity of the ice sheet in μS / m, f is the detection frequency, ε0 ​​is the vacuum dielectric constant, ε″ ice is the imaginary part of the relative dielectric constant of the ice sheet, c is the speed of light in a vacuum, k is the Boltzmann constant, T is the temperature of the ice sheet, T r is the reference temperature, [H + ]、[Cl - ]、 is the molar concentration of the corresponding ion, is the molar conductivity of the corresponding ions, is the activation energy of the corresponding ion.

[0148] Based on the prior information of ice sheet ion concentration distribution, the dielectric absorption attenuation rate of the ice sheet above the non-reflection zone is estimated by using the intensity gradient of the ice-penetrating radar echo profile, which can provide effective constraints on the internal temperature profile of the ice sheet above the non-reflection zone.

[0149] Then, the one-way dielectric absorption attenuation rate N above the ice sheet non-reflection area is obtained by using the ice penetrating radar echo intensity gradient. a After that, the total dielectric absorption attenuation of the ice sheet in the non-reflection area can be obtained: a In determining the ice bottom reflectivity R base After the effective range is determined, the radar echo intensity can be corrected according to the ice bottom normalization. The one-way reflection loss Γ on the ice sheet surface and internal layered interface, and the one-way total dielectric absorption attenuation L on the ice sheet above the non-reflection zone a , the two-way total dielectric absorption attenuation from the non-reflection area to the ice bottom is L a,EFZ Make an estimate

[0150]

[0151] The dielectric absorption attenuation of the ice sheet from the non-reflection zone to the ice bottom estimated using the ice-penetrating radar echo profile can provide effective constraints on the internal temperature profile of the ice sheet from the non-reflection zone to the ice bottom.

[0152] In step 5), the initial values ​​and probability density distribution of the parameters to be estimated of the temperature profile inside the ice sheet are optimized based on the dielectric absorption attenuation value of the ice sheet estimated by the ice-penetrating radar echo under the Bayesian estimation framework. Specifically, it includes:

[0153] Step 5-3) The estimated parameters of the internal temperature profile of the ice sheet and the corresponding temperature parameter vector x determined in step 3) are t , combined with the discrete temperature parameter vector x in step 4) t The simulated data set of ice-penetrating radar echo profile corresponding to the parameter grid is used to obtain the temperature parameter vector x t Simulated values ​​of ice sheet dielectric absorption attenuation corresponding to all parameter grid points in parameter space and

[0154] The ice sheet dielectric absorption attenuation value N estimated based on the ice penetrating radar echo profile observations in each interval of the target area a and L a,EFZ , estimate the temperature parameter vector x in each interval t All parameter grid points in parameter space The posterior probability

[0155]

[0156] in, is the temperature parameter grid point The corresponding prior probability density distribution.

[0157]

[0158] and They correspond to the estimated ice sheet dielectric absorption attenuation value N based on the ice penetrating radar echo profile observation value. a and L a,EFZ The estimated standard deviation of .

[0159] Posterior probability In the temperature parameter vector x t All parameter grid points in parameter space The calculation result is the temperature parameter vector x after the optimized ice sheet dielectric absorption attenuation value estimated by ice penetrating radar echo profile observation in each interval. t Probability density distribution on parameter space.

[0160] Step 5-4) Calculate the temperature parameter vector x in each interval t The estimated result of x t Estimated value of The calculation is based on the minimum mean square error method:

[0161]

[0162] x t Estimated value of The maximum a posteriori probability method is calculated as:

[0163]

[0164] Step 5-5) Used to update the temperature parameter vector x in each interval t All parameter grid points in parameter space Corresponding probability density; temperature parameter vector x t Estimated value of Used to update the initial values ​​of the parameters to be estimated for the temperature profile inside the ice sheet in each interval. The two are substituted as prior information into the subsequent Bayesian estimation process of the temperature profile parameters inside the ice sheet based on broadband radiation brightness temperature, and combined with the constraint information provided by active and passive methods, the active and passive joint inversion of the temperature profile parameters inside the ice sheet is realized.

[0165] In step 6), under the Bayesian estimation framework, according to the inversion problem description in step 3), the estimated parameters of the internal temperature profile of the ice sheet and the redundant estimated parameters are inverted based on the broadband radiation brightness temperature, wherein the initial values ​​and probability density distributions of the estimated parameters of the internal temperature profile of the ice sheet obtained by optimization in step 5) are used as prior information. Specifically, it includes:

[0166] Step 6-1) Based on the redundant parameters to be estimated and the corresponding redundant parameter vector x determined in step 3) n , combined with the discrete temperature parameter vector x in step 4) t The broadband radiation brightness temperature simulation data set corresponding to the parameter grid, as well as the broadband radiation brightness temperature observation vector y in each interval of the target area and the temperature parameter vector x obtained in step 5) t Estimated value of Calculate the redundant parameter vector x on each interval n All parameter grid points in parameter space The posterior probability

[0167]

[0168] Step 6-2) Calculate the temperature parameter vector x in each interval t All parameter grid points in parameter space The posterior probability

[0169]

[0170] Posterior probability In the temperature parameter vector x t All parameter grid points in parameter space The calculation result is the temperature parameter vector x after active and passive joint constraints in each interval. t Probability density distribution on parameter space.

[0171] Step 6-3) Calculate the temperature parameter vector x in each interval t The temperature parameter vector x is calculated using the minimum mean square error (MMSE) method. t Estimated value of

[0172]

[0173] Or use the maximum a posteriori probability method to calculate the temperature parameter vector x t Estimated value of

[0174]

[0175] The corresponding estimated variance is:

[0176]

[0177] In step 6), the internal temperature profile of the ice sheet in the target area is obtained according to the parameter inversion results and the internal temperature profile model of the ice sheet in step 3), which specifically includes:

[0178] Step 6-4) Calculate the estimated results of the internal temperature profile of the ice sheet in each interval, according to the temperature parameter vector x in each interval t All parameter grid points in parameter space The corresponding posterior probability The minimum mean square error (MMSE) and other methods are used to calculate the estimated value of the ice sheet temperature at any depth z, thereby obtaining the inversion result of the internal temperature profile of the ice sheet in the target area. The minimum mean square error (MMSE) calculation method is

[0179]

[0180] The corresponding estimated variance is

[0181]

[0182] The technical solution of the present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0183] Example 1

[0184] Embodiment 1 of the present invention proposes an active and passive joint inversion method of the internal temperature profile of the ice sheet based on Bayesian estimation, comprising the following steps:

[0185] Step 101: Select a target area and obtain broadband radiation brightness temperature and ice-penetrating radar echo profile observations.

[0186] (1) The area of ​​the Greenland ice sheet where the internal temperature profile of the ice sheet has been measured by drilling is selected as the target ice sheet area for inversion of the internal temperature profile of the ice sheet.

[0187] (2) Based on the existing broadband radiation brightness temperature and ice-penetrating radar echo profile detection data in the target area of ​​the Greenland ice sheet, the Camp Century to NGRIP section, which is covered by the ultra-wideband software-defined microwave radiometer airborne ultra-wideband radiation brightness temperature observations in the Greenland ice sheet area and the accumulation rate radar and multi-channel coherent sounding radar echo observations of the US IceBridge Program, was selected to conduct active and passive joint inversion of the temperature profile inside the ice sheet. Among them, the ultra-wideband software-defined microwave radiometer frequency band is 500-2000MHz; the accumulation rate radar frequency band is 600-900MHz; the multi-channel coherent sounding radar frequency band is 180-210MHz, which has effective detection capabilities for the temperature profile inside the ice sheet.

[0188] In the attached figure Figure 1 The selected detection area of ​​the Greenland ice sheet and the corresponding active and passive detection flight tracks are given. It can be seen that the flight track passes through the Camp Century borehole, NEEM borehole and NGRIP borehole areas, and the flight track is developed along the ice divide, corresponding to the dry snow area, and the influence of liquid water in the snow layer can be ignored.

[0189] (3) The flight path from Camp Century to NGRIP is divided into 118 intervals of about 5 km. Figure 2(a) shows the radiation brightness temperature observation data corresponding to the 12 frequency channels of the ultra-wideband software-defined microwave radiometer, and Figure 2(b) shows the average value in each interval. Figure 3 The radiation brightness temperature observation data corresponding to the interval where the NGRIP, NEEM and Camp Century boreholes are located are given, where the dot symbols represent the radiation brightness temperature observation data corresponding to the 12 frequency channels in the interval where the NGRIP, NEEM and Camp Century boreholes are located (interval length 5km), and the solid line represents the interval average value.

[0190] (4) According to the real part of the average relative dielectric constant of the ice sheet, the propagation speed of electromagnetic waves in the ice sheet is determined as follows:

[0191]

[0192] Among them, c is the speed of light, and v is the propagation speed of electromagnetic waves in the ice sheet, with the unit of m / s.

[0193] According to the propagation speed of electromagnetic waves in the ice sheet, the time series of ice-penetrating radar echo intensity is converted into the depth distribution of ice-penetrating radar echo intensity, so as to obtain the ice-penetrating radar echo profile observation values ​​in each interval. Figures 4(a), 4(b) and 4(c) in the attached drawings respectively show the echo profiles of multi-channel coherent sounding radar and accumulation rate radar in the intervals where NGRIP, NEEM and Camp Century boreholes are located.

[0194] Step 102: Extract the internal structure of the ice sheet and the thickness parameters of the ice sheet based on the ice-penetrating radar echo profile observations in step 101, and estimate the fluctuation characteristics of the snow layer density.

[0195] (1) Extract the reflection layer structure from the multi-channel coherent sounding radar and accumulation rate radar echo profile observations in each interval of the target area, and divide the internal structure of the ice sheet in the target area. Calculate the thickness of the ice sheet in each interval based on the depth of the target position reflected by the ice bottom and the position of the ice surface reflection. Figure 5 Results for ice sheet thickness in the target area along the flight track are given.

[0196] (2) Based on the empirical model of snow layer density fluctuation characteristics, combined with the detection system parameters of multi-channel coherent sounding radar and accumulation rate radar, the echo profiles of multi-channel coherent sounding radar and accumulation rate radar are forward simulated. By adjusting the snow layer density fluctuation characteristic parameters, the multi-channel coherent sounding radar and accumulation rate radar echo profile observations of the snow layer area in each interval are fitted, and the empirical model parameters of snow layer density fluctuation characteristics that achieve the best fit for the multi-channel coherent sounding radar and accumulation rate radar echo profile observations in each interval are determined to describe the snow layer density fluctuation characteristics in each interval of the target area. Figure 6 in the attached figure shows the single density sample generation results and the corresponding density fluctuation standard deviation distribution results corresponding to the snow layer density fluctuation standard deviation parameter estimation values ​​in the interval [1,41,75,118], where Figure 6(a) is interval 1, Figure 6(b) is interval 41, Figure 6(c) is interval 75, and Figure 6(d) is interval 118.

[0197] Step 103: Determine the ice sheet internal temperature profile model T(z), determine the parameters to be estimated and redundant parameters to be estimated of the ice sheet internal temperature profile, and describe the inversion problem. t,i (i=1,…,p) constitutes the temperature parameter vector x t , the superscript “T” indicates transpose.

[0198] T(z)=g(x t ,z)

[0199] x t =[x t,1 ,x t,2 ,…,x t,p ]T

[0200] (1) In this embodiment, since the flight track of the selected target area is developed along the ice divide, the horizontal heat convection of the ice sheet in the corresponding area reaches a minimum value. The Robin temperature model based on the solution of the one-dimensional heat conduction equation can meet the requirements of simulating the temperature profile inside the ice sheet.

[0201]

[0202] Where erf is the error function, T s is the annual average surface temperature in K, G is the geothermal flux in W / m 2 , M is the annual average accumulation rate, unit is m / year, H is the ice thickness, unit is m, k c =2.7W / m / K is the thermal conductivity of ice, k d =45m 2 / year is the ice thermal diffusion rate.

[0203] Since the ice sheet thickness H of each interval has been obtained in step 102, when the ice sheet thickness H can be obtained through existing products or active ice-penetrating radar data, the parameters to be estimated in the ice sheet internal temperature profile model include T s , M, G, the corresponding temperature parameter vector is

[0204] x t =[T s ,M,G] T

[0205] The superscript "T" indicates transpose.

[0206] (2) Redundant parameters to be estimated x n,i (i=1,…,q) constitutes the redundant parameter vector x n .

[0207] x n =[x n,1 ,x n,2 ,…,x n,q ] T

[0208] In step 102, the snow layer density fluctuation characteristics have been estimated based on the ice-penetrating radar echo profile observations. However, in the inversion process based on broadband radiation brightness temperature fitting, the correlation length parameters used to describe the density fluctuation scale often need to be further refined and adjusted to better fit the radiation brightness temperature observations. Therefore, in the inversion of the temperature profile inside the ice sheet, the snow layer density fluctuation standard deviation parameters estimated based on the ice-penetrating radar echo profile are used as known parameters, while the snow layer two-scale density fluctuation correlation length parameters l1, l2 and the small-scale density fluctuation scale factor α are used as redundant parameters to be estimated, corresponding to the redundant parameter vector x n for

[0209] x n =[l1,l2,α] T

[0210] (3) In the inversion problem of the internal temperature profile of the ice sheet, the internal structure of the ice sheet and the ice sheet thickness parameters extracted by the ice-penetrating radar echo profile in step 102 are used as known parameters. In addition, the density fluctuation standard deviation parameter in the snow layer density fluctuation characteristics estimated by the ice-penetrating radar echo profile in step 102 is also used as a known parameter; the ice-penetrating radar echo profile observation value is used as a constraint condition to optimize the initial value and probability density distribution function of the parameters to be estimated of the internal temperature profile of the ice sheet, which serves as the prior information in the inversion problem; finally, the broadband radiation brightness temperature observation value is used as the observation vector in the inversion problem to realize the Bayesian estimation of the parameters to be estimated and the redundant parameters to be estimated of the internal temperature profile of the ice sheet.

[0211] In the specific description, it is assumed that the broadband radiation brightness temperature observation With N frequency points The observation vector is defined as y, which is composed of broadband radiation brightness temperature observations and has a dimension of N. During the inversion process, the radiation brightness temperature simulation value T B The temperature parameter vector x corresponding to the temperature profile of the interior of the ice sheet is determined by the parameters to be estimated and the redundant parameters to be estimated. t The redundant parameter vector x corresponding to the redundant estimated parameters n Together they constitute the model parameter vector x in the inversion problem of the temperature profile inside the ice sheet.

[0212]

[0213] T B =f(x)

[0214] Temperature parameter vector x t Including the annual mean surface temperature T s , annual average accumulation rate M and geothermal flux G; redundant parameter vector x nIncluding the correlation length parameters l1, l2 of the two-scale density fluctuation of the snow layer and the proportional factor α of the small-scale density fluctuation; other parameters are assumed to be known parameters. The observation noise vector em is independent of the model parameter vector x and is assumed to be additive Gaussian noise. f(x) represents the simulation result corresponding to the radiation forward model under the model parameter vector x. On this basis, the inversion problem of the temperature profile inside the ice sheet can be expressed as:

[0215]

[0216]

[0217] y=f(x)+em

[0218] In this way, the likelihood function can be expressed as:

[0219]

[0220] Where C is the noise covariance matrix.

[0221] Step 104: Discretize the parameter space to be estimated determined in step 103, and perform active and passive forward simulations on the discrete parameter space to be estimated in combination with the snow layer density fluctuation characteristics, ice cover internal structure and ice thickness parameters obtained in step 2), to generate a simulated data set of broadband radiation brightness temperature and ice-penetrating radar echo profile.

[0222] (1) According to the effective range of the estimated parameters of the internal temperature profile of the ice sheet and the redundant estimated parameters, the temperature parameter vector x t and the redundant parameter vector x n The parameter space of the model parameter vector x is discretized respectively to achieve the discretization of the parameter space and generate the corresponding parameter grid.

[0223] In this embodiment, the temperature parameter vector x t In the parameter space grid division, the annual average surface temperature T s The effective range of the value of is 240 to 249K, the effective range of the annual average accumulation rate M is 0.2 to 0.38m / year, and the effective range of the geothermal flux G is 50 to 95mW / m 2 The annual mean surface temperature T s , the temperature parameter vector x formed by the annual average accumulation rate M and the geothermal flux G t The parameter space is discretized into 17576 grid points (26×26×26).

[0224] In this embodiment, the redundant parameter vector x nIn the parameter space grid division, the effective value range of the large-scale density fluctuation correlation length l1 of the snow layer is 6 to 15 cm, the effective value range of the small-scale density fluctuation correlation length l2 is 1 to 5 cm, and the effective value range of the small-scale density fluctuation proportional factor α is 0.1 to 0.9. The redundant parameter vector x is formed by the two-scale density fluctuation correlation length l1, l2 of the snow layer and the small-scale density fluctuation proportional factor α. n The parameter space is discretized into 4394 grid points (26×13×13).

[0225] (2) Combined with the standard deviation parameters of snow layer density fluctuation, the internal structure of the ice sheet, and the ice thickness parameters obtained in step 102, the broadband radiation brightness temperature simulation values ​​and ice-penetrating radar echo profile simulation values ​​corresponding to all temperature parameter vector grid points and redundant parameter vector grid points in the 118 analysis intervals from Camp Century to NGRIP are calculated, and the simulated data sets of broadband radiation brightness temperature and ice-penetrating radar echo profile in each interval are obtained.

[0226] Step 105: estimating the dielectric absorption attenuation value of the ice sheet according to the ice-penetrating radar echo profile observation value in step 101, and optimizing the initial value and probability density distribution of the parameters to be estimated of the temperature profile inside the ice sheet based on the dielectric absorption attenuation value of the ice sheet estimated by the ice-penetrating radar echo under the Bayesian estimation framework;

[0227] (1) The multi-channel coherent sounding radar selected in this embodiment is in the VHF frequency band and can detect the reflected echo from the bottom of the ice in the target area. Therefore, this embodiment will use the VHF frequency band multi-channel coherent sounding radar echo to estimate the dielectric absorption attenuation value of the ice cover in each interval from Camp Century to the NGRIP target area.

[0228] The ice-penetrating radar echo profile corresponds to the depth distribution of the ice-penetrating radar echo intensity from the ice sheet surface to the ice bottom, which is mainly determined by the radar system parameters, the reflectivity of the interface from the ice sheet surface to the ice bottom, the dielectric absorption attenuation of the ice sheet, and the spatial attenuation during the propagation of electromagnetic waves. Under the quasi-specular reflection approximation, the radar echo intensity can be expressed as

[0229]

[0230] in, is the intensity of the i-th reflected echo, P t is the transmission power, G t and G R are the transmitting and receiving gains respectively, λ is the radar wavelength, L a,i is the ice sheet one-way dielectric absorption attenuation experienced by the i-th reflected echo, Γ i is the total one-way transmittance of the ice sheet interface experienced by the i-th reflected echo, corresponding to the one-way reflection loss of the ice sheet surface and internal layered interface, Lsys is the system loss, G p is the system gain, R i The reflectivity corresponding to the i-th reflected echo, h is the height of the observation platform, ε′ ice is the real part of the average relative permittivity of the ice sheet, z i The depth of the position corresponding to the i-th reflection echo.

[0231] According to the multi-channel coherent bathymetric radar platform height h, reflection depth position z and the real part of the average relative dielectric constant ε′ of the ice sheet ice , the spatial attenuation in the echo intensity of the multi-channel coherent sounding radar is compensated, and the intensity of the i-th reflected echo after compensation is

[0232]

[0233] Relative to surface echo intensity Reflection echo intensity Normalize and compensate for the radar system parameter P t , G t , G R ,λ,L sys and G p At this time, the normalized i-th reflection echo intensity The dielectric absorption of the ice sheet attenuates L a,i , ice sheet reflection loss Γ i and the reflectivity R at the echo position i The decision can be expressed in the logarithmic domain as

[0234]

[0235] (2) In this embodiment, the echo gradient method is used to estimate the dielectric absorption attenuation parameter of the ice sheet. The dielectric absorption attenuation rate of the ice sheet above the ice sheet non-reflection area is corrected by normalizing the radar echo intensity To avoid the influence of the reflection characteristics of the snow layer, the shallowest position z in the depth range of the gradient solution is min Not less than the thickness of the snow layer. In this embodiment, z min It is set to 200m. In addition, considering that the echo intensity of the ice-penetrating radar near the bottom of the ice is weak, the signal-to-noise ratio is low, and it may introduce deviations in the estimation of the dielectric absorption attenuation rate of the ice sheet. Therefore, the deepest position z in the depth range of the gradient solution is max The depth position of the upper boundary of the non-reflection zone above the ice bottom interface should not exceed that of the upper boundary of the non-reflection zone above the ice bottom interface. The depth gradient can be expressed as

[0236]

[0237] The one-way reflection loss Γ on the ice sheet surface and internal stratified interface is mainly generated by the shallow snow layer area. min and z max In the ice sheet region between Γ and Γ, Γ can be considered as a constant, so d[Γ] dB / dz can be ignored. In addition, the internal reflection below the snow layer of the ice sheet is mainly caused by the sudden change of conductivity or the change of the ice crystal axis orientation. The corresponding reflectivity R does not change significantly with depth, and its contribution to the depth gradient of radar echo intensity can be almost ignored. Therefore, d[R] dB The / dz term is also approximately equal to 0. In this way, the normalized corrected radar echo intensity The depth gradient is approximately equal to twice the depth gradient of the one-way dielectric absorption attenuation of the ice sheet. The depth gradient of the one-way dielectric absorption attenuation of the ice sheet is the one-way dielectric absorption attenuation rate of the ice sheet N. a .

[0238]

[0239] Ice sheet one-way dielectric absorption attenuation rate N a is a function of ice sheet temperature and ice chemistry

[0240]

[0241] Among them, N a The unit is dB / km, σ ∞ is the conductivity of the ice sheet in μS / m, f is the detection frequency, ε0 ​​is the vacuum dielectric constant, ε″ ice is the imaginary part of the relative dielectric constant of the ice sheet, c is the speed of light in a vacuum, k is the Boltzmann constant, T is the temperature of the ice sheet, T r is the reference temperature, [H + ]、[Cl - ]、 is the molar concentration of the corresponding ion, is the molar conductivity of the corresponding ions, is the activation energy of the corresponding ion.

[0242] The single-pass dielectric absorption attenuation rate N above the ice sheet non-reflection area is obtained by using the multi-channel coherent sounding radar echo intensity gradient. a After that, the total dielectric absorption attenuation of the ice sheet in the non-reflection area can be obtained: a In determining the ice bottom reflectivity R base After the effective range is determined, the radar echo intensity can be corrected according to the ice bottom normalization. The one-way reflection loss Γ on the ice sheet surface and internal layered interface, and the one-way total dielectric absorption attenuation L on the ice sheet above the non-reflection zone a , the two-way total dielectric absorption attenuation from the non-reflection area to the ice bottom is L a,EFZMake an estimate.

[0243]

[0244] (3) Based on the estimated parameters of the internal temperature profile of the ice sheet and the corresponding temperature parameter vector x determined in step 103 t , combined with the discrete temperature parameter vector x in step 104 t The simulated data set of ice-penetrating radar echo profile corresponding to the parameter grid is used to obtain the temperature parameter vector x t Simulated values ​​of ice sheet dielectric absorption attenuation corresponding to all parameter grid points in parameter space and

[0245] The ice sheet dielectric absorption attenuation value N is estimated based on the multi-channel coherent bathymetric radar echo profile observations in each section from Camp Century to the NGRIP target area. a and L a,EFZ , estimate the temperature parameter vector x in each interval t All parameter grid points in parameter space The posterior probability

[0246]

[0247] in, is the temperature parameter grid point The corresponding prior probability density distribution. In this embodiment, it is assumed that the annual average surface temperature T s The effective value range given in step 104 satisfies the Gaussian distribution, and the ice sheet surface temperature calculated by the atmospheric climate model is taken as T s Prior probability density distribution p(T s ) with a standard deviation of 1K; assuming that the annual average accumulation rate M and geothermal flux G are uniformly distributed within the valid value range given in step 104.

[0248]

[0249] and They correspond to the estimated ice sheet dielectric absorption attenuation value N based on the ice penetrating radar echo profile observation value. a and L a,EFZ The estimated standard deviation of .

[0250] Posterior probability In the temperature parameter vector x t All parameter grid points in parameter space The calculation is performed on the intervals (a total of 17576, 26×26×26), and the calculation result is the temperature parameter vector x after the ice sheet dielectric absorption attenuation value is optimized by ice penetrating radar echo profile observation. t Probability density distribution on parameter space.

[0251] Figure 7 in the attached figure shows the initial prior probability density distribution of the parameters to be estimated for the temperature profile inside the ice sheet in the intervals where the NGRIP, NEEM and Camp Century boreholes are located, as well as the probability density distribution after optimization using the ice sheet dielectric absorption attenuation value estimated based on the multi-channel coherent sounding radar echo profile. The horizontal axis of Figure 7(a) is the annual average surface temperature, the horizontal axis of Figure 7(b) is the annual average accumulation rate, and the horizontal axis of Figure 7(c) is the geothermal flux; it can be seen that due to the annual average surface temperature T s The prior probability density distribution description itself is relatively accurate. Therefore, the ice sheet dielectric absorption attenuation value estimated based on the ice-penetrating radar echo profile mainly improves the constraint ability of the annual average accumulation rate M and the geothermal flux G. Compared with the uniform distribution itself, which does not contain effective constraint information, the ice sheet dielectric absorption attenuation value estimated based on the ice-penetrating radar echo profile can optimize the probability density distribution of the annual average accumulation rate M and the geothermal flux G, thereby improving the constraint ability of the temperature profile inside the ice sheet in the subsequent inversion process based on broadband radiation brightness temperature.

[0252] (4) Calculate the temperature parameter vector x in each interval from Camp Century to NGRIP t The estimated result of x t Estimated value of The calculation is based on the minimum mean square error method:

[0253]

[0254] x t Estimated value of The maximum a posteriori probability method is calculated as:

[0255]

[0256] (5) Used to update the temperature parameter vector x in each interval t All parameter grid points in parameter space Corresponding probability density; temperature parameter vector x t Estimated value of Used to update the initial values ​​of the parameters to be estimated for the temperature profile inside the ice sheet in each interval. The two are substituted as prior information into the subsequent Bayesian estimation process of the temperature profile parameters inside the ice sheet based on broadband radiation brightness temperature, and combined with the constraint information provided by active and passive methods, the active and passive joint inversion of the temperature profile parameters inside the ice sheet is realized.

[0257] Step 106: Under the Bayesian estimation framework, according to the inversion problem description in step 103, the parameters to be estimated and the redundant parameters to be estimated of the internal temperature profile of the ice sheet are inverted based on the broadband radiation brightness temperature, wherein the initial values ​​and probability density distribution of the parameters to be estimated of the internal temperature profile of the ice sheet obtained by optimization in step 105 are used as prior information. The internal temperature profile of the ice sheet in the target area is obtained based on the parameter inversion results and the internal temperature profile model of the ice sheet in step 103.

[0258] (1) Based on the redundant parameters to be estimated and the corresponding redundant parameter vector x determined in step 103 n , combined with the discrete temperature parameter vector x in step 104 t The broadband radiation brightness temperature simulation data set corresponding to the parameter grid, as well as the broadband radiation brightness temperature observation vector y in each interval from Camp Century to the NGRIP target area and the temperature parameter vector x obtained in step 105 t Estimated value of Calculate the redundant parameter vector x on each interval n All parameter grid points in parameter space The posterior probability

[0259]

[0260] (2) Calculate the temperature parameter vector x in each interval from Camp Century to NGRIP t All parameter grid points in parameter space The posterior probability

[0261]

[0262] Posterior probability In the temperature parameter vector x t All parameter grid points in parameter space The calculation result is the temperature parameter vector x after active and passive joint constraints on each section from Camp Century to NGRIP. t Probability density distribution on parameter space.

[0263] (3) Calculate the temperature parameter vector x in each interval from Camp Century to NGRIP t The temperature parameter vector x is calculated using the minimum mean square error (MMSE) method.t Estimated value of

[0264]

[0265] Or use the maximum a posteriori probability method to calculate the temperature parameter vector x t Estimated value of

[0266]

[0267] The corresponding estimated variance is:

[0268]

[0269] Figure 8 in the attached figure shows the inversion results of the estimated parameters of the temperature profile inside the ice sheet on the flight track from Camp Century to NGRIP, where Figure 8(a) is the annual average surface temperature, Figure 8(b) is the annual average accumulation rate, and Figure 8(c) is the geothermal flux. The annual average surface temperature T s The average estimated standard deviation is about 0.65K, and the average estimated standard deviations of the annual accumulation rate M and geothermal flux G are about 0.031m / year and 0.0059W / m 2 . In order to illustrate the significance of the combined active detection channel for improving the accuracy of the inversion of the internal temperature profile of the ice sheet, Figure 8 also shows the estimated parameters of the internal temperature profile of the ice sheet obtained by inverting the broadband radiation brightness temperature detection channel alone. It can be seen that the dielectric absorption attenuation value of the ice sheet estimated based on the ice-penetrating radar echo profile can significantly reduce the estimated standard deviation of the annual average accumulation rate M and geothermal heat flux G in the parameters to be estimated in the internal temperature profile of the ice sheet. Among them, the average estimated standard deviation of the annual average accumulation rate M is reduced from 0.057m / year to 0.031m / year, and the average estimated standard deviation of the geothermal heat flux G is reduced from 0.014W / m 2 Reduced to 0.0059W / m 2 , and can slightly reduce the annual average surface temperature T s The average estimated standard deviation of (from 0.88K to 0.65K).

[0270] (4) Calculate the estimated temperature profile of the ice sheet in each interval from Camp Century to NGRIP, based on the temperature parameter vector x t All parameter grid points in parameter space The corresponding posterior probability The minimum mean square error (MMSE) method is used to calculate the estimated value of the ice sheet temperature at any depth z, thereby obtaining the inversion result of the temperature profile inside the ice sheet in the target area.

[0271]

[0272] The corresponding estimated variance is

[0273]

[0274] Figure 9(a) in the attached figure shows the inversion results of the internal temperature profile of the ice sheet on the flight track from Camp Century to NGRIP under the active and passive joint inversion and the corresponding estimated standard deviation, and Figure 9(b) shows the inversion results using only the passive broadband radiation brightness temperature detection channel and the corresponding estimated standard deviation. It can be seen that in the case of combining the ice sheet dielectric absorption attenuation value estimated based on the ice-penetrating radar echo profile, the inversion accuracy of the ice sheet temperature at different depths in different intervals is better than the inversion results of a single passive channel.

[0275] In the attached figure Fig.10 The average estimated standard deviation of the ice sheet temperature at different depths along the flight track from Camp Century to NGRIP is given. It can be seen that within the depth range of 1.5 km, the average estimated standard deviation of the ice sheet temperature based on the active and passive joint inversion is within 1K. With the increase in depth, the average estimated standard deviation of the ice sheet temperature also increases slightly. When it reaches the depth range of 3 km, the average estimated standard deviation of the ice sheet temperature is about 3.3K, and the average estimated standard deviation of the ice sheet temperature in the study area within the total depth range is about 1.2K. In contrast, the average estimated standard deviation of the ice sheet temperature based on a single passive radiation brightness temperature detection channel is higher than that of the active and passive joint inversion. Especially for ice sheet areas with a depth of more than 1.5 km, the inversion capability of deep ice sheet temperature with large uncertainty can be effectively improved by combining the ice sheet dielectric absorption attenuation constraint information estimated based on the ice penetrating radar echo profile.

[0276] Figure 11 in the attached figure also compares the inversion results of the internal temperature profile of the ice sheet in the interval where the NGRIP, NEEM, and Camp Century boreholes are located with the borehole measurement results of the internal temperature profile of the ice sheet fitted by the Robin model, where Figure 11(a) is NGRIP, Figure 11(b) is NEEM, and Figure 11(c) is Camp Century. It can be seen that the internal temperature profile of the ice sheet obtained by the active and passive joint inversion can well reflect the vertical distribution characteristics of the actual internal temperature of the ice sheet. This shows that the complementary constraint information provided by the active and passive joint detection method can effectively improve the inversion ability of the internal temperature profile of the ice sheet, and while reducing the standard deviation of the temperature estimation at different depths, it can also effectively reduce the absolute deviation between the inversion results of the internal temperature profile of the ice sheet and the actual vertical distribution of the internal temperature of the ice sheet.

[0277] In summary, this embodiment provides an active and passive joint inversion method for the internal temperature profile of the ice sheet based on Bayesian estimation. The specific method flow is as follows: Fig.12 As shown. Based on the sensitivity of the passive broadband radiation brightness temperature observation channel to the internal temperature profile of the ice sheet, the ice-penetrating radar echo profile characteristics are further combined. On the one hand, the snow layer density fluctuation characteristic information and the internal structure information of the ice sheet constrained by the reflection characteristics of the ice-penetrating radar echo profile are introduced to reduce the uncertainty of radiation interpretation; on the other hand, according to the correlation between the dielectric absorption attenuation of the ice sheet and the internal temperature of the ice sheet, the dielectric absorption attenuation value of the ice sheet estimated based on the ice-penetrating radar echo profile is used to further constrain the internal temperature profile of the ice sheet. In the specific inversion algorithm, the dielectric absorption attenuation value of the ice sheet estimated based on the ice-penetrating radar echo profile is used as a constraint condition to optimize the initial value and probability density distribution of the parameters to be estimated for the internal temperature profile of the ice sheet. The optimization result is used as the prior information of the internal temperature profile of the ice sheet and substituted into the inversion process based on broadband radiation brightness temperature observation, so that the constraint information provided by the active and passive can be combined to realize the active and passive joint inversion of the internal temperature profile of the ice sheet on the basis of Bayesian estimation. By using the method provided in this embodiment, it is possible to effectively reduce or eliminate the interference introduced by the snow layer density fluctuations and the internal structure of the ice sheet in the inversion of the temperature profile inside the ice sheet. At the same time, it is possible to utilize the joint constraint information provided by the active and passive detection channels to improve the effectiveness and reliability of the inversion of the temperature profile inside the ice sheet, and effectively reduce the dimensionality and uncertainty of the problem of inversion of the temperature profile inside the ice sheet.

[0278] Example 2

[0279] The embodiment of the present invention may also provide a computer program product, including a computer program / instruction. When the computer program / instruction is executed by a processor, each step in the above method embodiment can be implemented.

[0280] Summary:

[0281] The method of the present invention can obtain the internal temperature profile information of the ice sheet, and based on Bayesian estimation, uses ice-penetrating radar echoes and broadband radiation brightness temperature to perform active and passive joint inversion of the internal temperature profile of the ice sheet. This changes the current situation that the internal temperature profile of the ice sheet can only be detected through on-site drilling technology measurement or single passive or active remote sensing means, overcomes the interference of the internal structure of the ice sheet and the snow layer density fluctuation on the inversion of the internal temperature profile of the ice sheet in single passive remote sensing means, and overcomes the defect that the single active remote sensing means can only estimate the average temperature of the ice sheet and cannot provide information on the vertical change of the internal temperature of the ice sheet.

[0282] The method provided by the present invention makes full use of the advantages of the joint constraint information provided by active and passive detection, and effectively reduces the uncertainty of ice sheet radiation interpretation by introducing the snow layer density fluctuation characteristics and ice sheet internal structure information estimated by ice penetrating radar echo. On this basis, the dielectric absorption attenuation characteristics of the ice sheet based on the ice penetrating radar echo profile are used as constraints to optimize and update the initial values ​​and probability density distribution of the parameters to be estimated of the ice sheet internal temperature profile, and introduce them as prior information in the process of inversion of the ice sheet internal temperature profile based on broadband radiation brightness temperature observation, and fully combine the constraint information provided by active and passive detection to achieve active and passive joint inversion of the ice sheet internal temperature profile, effectively improve the inversion performance of the ice sheet internal temperature profile, reduce the dimension and uncertainty of the ice sheet internal temperature profile inversion problem, and improve the reliability and robustness of the ice sheet internal temperature profile inversion results.

[0283] The method proposed in the present invention can be applied to actual ice sheet areas based on the existing data conditions of the polar ice sheet to invert the temperature profile inside the ice sheet, which can be used to guide the improvement and research direction of subsequent active and passive joint detection methods of the temperature profile inside the ice sheet, and provide an algorithmic basis for promoting the development of remote sensing detection of the temperature profile inside the ice sheet in the future, so as to meet the needs of research on the thermodynamic processes and material balance evolution of the ice sheet, and fill the gap in the information on the temperature profile inside the ice sheet and its spatial distribution characteristics.

[0284] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit the present invention. Although the present invention is described in detail with reference to the embodiments, it should be understood by those skilled in the art that any modification or equivalent replacement of the technical solutions of the present invention does not depart from the spirit and scope of the technical solutions of the present invention and should be included in the scope of the claims of the present invention.

Claims

1. An active and passive joint inversion method for the internal temperature profile of the ice sheet based on Bayesian estimation, comprising: Step 1) Obtain broadband radiation brightness temperature and ice penetrating radar echo profile observations of the ice sheet area to be studied; Step 2) extracting the internal structure and thickness parameters of the ice sheet based on the ice-penetrating radar echo profile observations, and estimating the fluctuation characteristics of the snow layer density; Step 3) determining the ice sheet internal temperature profile model, and determining the ice sheet internal temperature profile parameters to be estimated and redundant parameters to be estimated, and performing inversion problem description; Step 4) Discretize the parameter space to be estimated, combine the snow layer density fluctuation characteristics, the internal structure of the ice cover and the ice thickness parameters, perform active and passive forward simulations on the discrete parameter space to be estimated, and generate a simulated data set of broadband radiation brightness temperature and ice penetrating radar echo profile; Step 5) estimating the dielectric absorption attenuation value of the ice sheet according to the ice penetrating radar echo profile observation value, and optimizing the initial value and probability density distribution of the parameters to be estimated of the temperature profile inside the ice sheet based on the dielectric absorption attenuation value of the ice sheet estimated by the ice penetrating radar echo under the Bayesian estimation framework; Step 6) Based on the broadband radiation brightness temperature under the Bayesian estimation framework, the estimated parameters and redundant estimated parameters of the internal temperature profile of the ice sheet are inverted. During the inversion, the initial values ​​and probability density distribution of the estimated parameters of the internal temperature profile of the ice sheet optimized in step 5) are used as prior information. According to the parameter inversion results and the internal temperature profile model of the ice sheet in step 3), the internal temperature profile of the ice sheet in the target area is obtained.

2. The active and passive joint inversion method of the internal temperature profile of the ice sheet based on Bayesian estimation according to claim 1 is characterized in that: The step 1) comprises: Along the ice sheet area to be studied, the intervals are divided by active ice-penetrating radar detection and passive broadband radiometer detection segments, and the broadband radiation brightness temperature observation values ​​and ice-penetrating radar echo intensity time series observation values ​​of each interval are obtained. The electromagnetic wave propagation speed v in the ice sheet is determined according to the following formula: Where c is the speed of light, v is in m / s, ε ′ is the real part of the average relative permittivity of the ice sheet; According to v, the time series of ice-penetrating radar echo intensity is converted into the depth distribution of ice-penetrating radar echo intensity, so as to obtain the ice-penetrating radar echo profile observation values ​​in each interval.

3. The active and passive joint inversion method of the internal temperature profile of the ice sheet based on Bayesian estimation according to claim 2 is characterized in that: The step 2) comprises: Step 2-1) extract the reflection target position depth in the ice penetrating radar echo profile observation value of each interval in the ice sheet area to be studied, divide the internal structure of the ice sheet, and calculate the ice sheet thickness in each interval according to the reflection target position depth of the ice bottom and the reflection position of the ice surface; Step 2-2) Based on the empirical model of snow layer density fluctuation characteristics and combined with the parameters of the ice-penetrating radar detection system, the ice-penetrating radar echo profile is forward simulated. By adjusting the snow layer density fluctuation characteristic parameters, the ice-penetrating radar echo profile observations of the snow layer area in each interval are fitted, and the empirical model parameters of snow layer density fluctuation characteristics that achieve the best fit to the ice-penetrating radar echo profile observations in each interval are determined.

4. The active and passive joint inversion method of the internal temperature profile of the ice sheet based on Bayesian estimation according to claim 3 is characterized in that: The step 3) comprises: Step 3-1) Based on the borehole measurement data of the internal temperature profile of the ice sheet, establish the internal temperature profile model T(z) of the ice sheet through polynomial fitting, moment method or by solving the heat conduction equation: T(z)=g(x t ,z) x t =[x t,1 ,x t,2 ,…,x t,p ] T Among them, z is the position depth, and the corresponding parameter to be estimated is x t,1 ,x t,2 ,…,x t,p Construct the temperature parameter vector x t , p represents the total number of parameters to be estimated, and the superscript "T" represents transposition; Step 3-2) From the redundant parameters to be estimated x n,i Constructing the redundant parameter vector x n : x n =[x n,1 ,x n,2 ,…,x n,q ] T Where q represents the total number of redundant parameters; Step 3-3) The inversion problem of the temperature profile inside the ice sheet is expressed as: y=f(x)+em in, For broadband radiation brightness temperature observation, with N frequency points y is the observation vector, which is composed of broadband radiation brightness temperature observations and has a dimension of N. f(x) represents the simulation result corresponding to the radiation forward model under the model parameter vector x, and em is the observation noise vector; Then the likelihood function L(x) is: Where C is the noise covariance matrix.

5. The active and passive joint inversion method of the internal temperature profile of the ice sheet based on Bayesian estimation according to claim 4 is characterized in that: The step 4) comprises: Step 4-1) According to the effective value range of the parameter to be estimated, the temperature parameter vector x t and the redundant parameter vector x n The parameter space of the model parameter vector x is discretized to generate the corresponding parameter grid. Step 4-2) combines the snow layer density fluctuation characteristics, ice sheet internal structure and ice thickness parameters obtained in step 2) to perform forward simulation of the ice sheet broadband radiation brightness temperature and ice-penetrating radar echo profile at discrete grid points in the model parameter space to generate the ice sheet broadband radiation brightness temperature and ice-penetrating radar echo profile grid; and perform interpolation and refinement to generate a simulated data set of broadband radiation brightness temperature and ice-penetrating radar echo profile.

6. The active and passive joint inversion method of the internal temperature profile of the ice sheet based on Bayesian estimation according to claim 5 is characterized in that: The step 5) comprises: Step 5-1) The estimated parameters of the internal temperature profile of the ice sheet and the corresponding temperature parameter vector x determined in step 3) t , combined with the discrete temperature parameter vector x in step 4) t The simulated data set of ice-penetrating radar echo profile corresponding to the parameter grid is used to obtain the temperature parameter vector x t Simulated values ​​of one-way ice sheet dielectric absorption attenuation corresponding to all parameter grid points in parameter space and the simulated value of the total dielectric absorption attenuation of the one-way ice sheet Step 5-2) The dielectric absorption attenuation value N of the ice sheet estimated based on the ice penetrating radar echo profile observations in each interval of the target area is a and the total dielectric absorption attenuation of the ice sheet L a,EFZ , estimate the temperature parameter vector x in each interval t All parameter grid points in parameter space The posterior probability Step 5-3) The posterior probability In the temperature parameter vector x t All parameter grid points in parameter space Calculate on to get the temperature parameter vector x t Probability density distribution on parameter space; Step 5-4) Calculate the temperature parameter vector x using the minimum mean square error or maximum a posteriori probability method t Estimated value of Step 5-5) Used to update the temperature parameter vector x in each interval t All parameter grid points in parameter space The corresponding probability density is the temperature parameter vector x t Estimated value of It is used to update the initial values ​​of the estimated parameters of the internal temperature profile of the ice sheet in each interval. The two are substituted as prior information into the subsequent Bayesian estimation process of the temperature profile parameters inside the ice sheet based on broadband radiation brightness temperature, and the constraint information provided by active and passive methods is combined to realize the active and passive joint inversion of the temperature profile parameters inside the ice sheet.

7. The active and passive joint inversion method of the internal temperature profile of the ice sheet based on Bayesian estimation according to claim 6 is characterized in that: Step 5-2) posterior probability for: in, is the temperature parameter grid point The corresponding prior probability density distribution; in, and They correspond to the estimated ice sheet dielectric absorption attenuation value N based on the ice penetrating radar echo profile observation value. a and L a,EFZ The estimated standard deviation of Step 5-4) t Estimated value of The calculation is based on the minimum mean square error method: x t Estimated value of The maximum a posteriori probability method is calculated as:

8. The active and passive joint inversion method of the internal temperature profile of the ice sheet based on Bayesian estimation according to claim 7 is characterized in that: The step 6) comprises: Step 6-1) Based on the redundant parameters to be estimated and the corresponding redundant parameter vector x determined in step 3) n , combined with the discrete temperature parameter vector x in step 4) t The broadband radiation brightness temperature simulation data set corresponding to the parameter grid, as well as the broadband radiation brightness temperature observation vector y in each interval of the target area and the temperature parameter vector x obtained in step 5) t Estimated value of Calculate the redundant parameter vector x on each interval n All parameter grid points in parameter space The posterior probability Step 6-2) Calculate the temperature parameter vector x in each interval t All parameter grid points in parameter space The posterior probability Posterior probability In the temperature parameter vector x t All parameter grid points in parameter space The calculation result is the temperature parameter vector x after active and passive joint constraints in each interval. t Probability density distribution on parameter space; Step 6-3) Calculate the temperature parameter vector x using the minimum mean square error method t Estimated value of Or use the maximum a posteriori probability method to calculate the temperature parameter vector x t Estimated value of Corresponding estimated variance for: Step 6-4) Calculate the estimated results of the internal temperature profile of the ice sheet in each interval, according to the temperature parameter vector x in each interval t All parameter grid points in parameter space The corresponding posterior probability The minimum mean square error method is used to calculate the estimated value of the ice sheet temperature at any depth z, thereby obtaining the inversion result of the internal temperature profile of the ice sheet in the target area; among them, the minimum mean square error Calculated as: Corresponding estimated variance for: in, is the temperature parameter vector x in each interval t All parameter grid points in the parameter space, j represents the jth grid point.

9. A computer program product, comprising a computer program, characterized in that When the computer program is executed by a processor, the steps of the method according to claim 1 are implemented.

Images