A method for inverting glacier thickness using aviation gravity anomaly

Abstract

The present invention discloses a method for inverting glacier thickness using aerial gravity anomaly, which belongs to the technical field of aerial gravity measurement, and comprises the following steps: obtaining measured aerial gravity data, and performing potential field separation on the aerial gravity data to obtain gravity anomaly data caused by ice-rock interface; according to the measured aerial gravity data, the glacier is divided into a plurality of horizontally juxtaposed rectangular cylinders; the ice-rock interface elevation model on the undulating aerial gravity observation surface is fitted with forward calculation of gravity anomaly, and the gravity anomaly of each observation point is calculated; the objective function of regularized inversion is constructed to obtain the correction amount of the ice-rock interface elevation model, and the ice-rock interface elevation is iteratively updated; the iterative process is repeated until the correction amount of the ice-rock interface elevation is less than a threshold value, and the iteration is stopped to obtain the ice-rock interface elevation. The present invention can solve the problems of model division and glacier thickness calculation when the upper interface is undulating and gravity anomaly calculation on the undulating observation surface.

Description

Technical Field

[0001] The invention relates to the technical field of aerial gravity measurement, and in particular to a method for inverting glacier thickness using aerial gravity anomaly. Background Art

[0002] There are nearly 48,000 glaciers in my country, but less than 50 of them have thickness measurements. Usually, thickness measurements rely on ground-penetrating radar, but since most of my country's glaciers are mountain glaciers, the measurement is difficult and the efficiency is low, and the amount of sampled data is difficult to support the calculation of glacier thickness over a large area. With the continuous development and maturity of aerial gravity measurement technology, it is possible to use aerial gravity anomaly data to invert the ice-rock density interface, and obtain the elevation of the ice-rock interface through inversion calculation, and then obtain the glacier thickness. Elevation refers to the distance from a point along the plumb line to the absolute base surface.

[0003] The traditional spatial domain density interface inversion method can obtain a relatively accurate density interface depth in a sedimentary basin. This method assumes that the upper interface of the sedimentary basin is a plane with the z axis pointing downward. The sedimentary layer is divided into M two-dimensional rectangular prisms to simulate the basement undulation of the sedimentary basin. The width of the two-dimensional prism is known and is a constant, and the thickness is the parameter p. j , the top surface of the two-dimensional rectangular cylinder coincides with the upper interface of the sedimentary basin, and the bottom surface coincides with the basement of the sedimentary basin (e.g. Figure 1 shown).

[0004] The gravity anomaly generated by M two-dimensional juxtaposed rectangular cylinders at the observation point can approximately represent the gravity anomaly caused by the basement undulation of the sedimentary basin at the observation point, that is,

[0005]

[0006] Among them, f i (p) represents the i-th observation point P(x i ) gravity anomaly, ξ j1 , j2 are the horizontal coordinates of the jth cylinder, p j represents the thickness of the jth cylinder, △g(x i ,p j ) represents the gravity anomaly caused by the j-th cylinder at the i-th observation point.

[0007] However, this method cannot obtain reliable glacier thickness when directly applied to density interface inversion based on airborne gravity anomalies, which affects the actual application effect. The reason is that, on the one hand, airborne gravity measurement obtains data on an undulating observation surface at a certain height above the ground, and on the other hand, the surface of mountain glaciers fluctuates greatly with the terrain, that is, the upper interface is not a plane.

[0008] The information disclosed in this background technology section is only intended to enhance the understanding of the overall background of the invention and should not be regarded as an acknowledgment or any form of suggestion that the information constitutes the prior art already known to a person skilled in the art. Summary of the invention

[0009] The purpose of the present invention is to provide a method for inverting glacier thickness using aerial gravity anomalies. It aims to address the technical problem that the spatial domain inversion method in the prior art cannot be applied to aerial gravity data, and can solve the problems of model segmentation and glacier thickness calculation when the upper interface is undulating, as well as gravity anomaly calculation on the undulating observation surface.

[0010] To achieve the above object, the present invention provides a method for inverting glacier thickness using aerial gravity anomaly, comprising the following steps:

[0011] Obtain measured aerial gravity data, perform potential field separation on the aerial gravity data, remove the background field, and obtain gravity anomaly data caused by the ice-rock interface; among them, the aerial gravity data obtained at the observation points on the undulating aerial gravity observation surface;

[0012] According to the measured aerial gravity data, the glacier was divided into multiple horizontally juxtaposed rectangular columns, and the initial elevation of the ice-rock interface was set as the glacier surface elevation.

[0013] The forward calculation of gravity anomaly is performed by fitting the ice-rock interface elevation model on the undulating aerial gravity observation surface, and the gravity anomaly of each observation point is calculated;

[0014] The objective function of regularized inversion is constructed, and the nonlinear objective function is linearized to obtain the correction value of the ice-rock interface elevation model, and the ice-rock interface elevation is iteratively updated;

[0015] Repeat the iterative process until the correction amount of the ice-rock interface elevation is less than the threshold, then stop the iteration and obtain the ice-rock interface elevation;

[0016] The thickness of the glacier is calculated based on the elevation of the glacier surface and the elevation of the ice-rock interface.

[0017] In one embodiment of the present invention, dividing the glacier into a plurality of horizontally juxtaposed rectangular cylinders according to measured aerial gravity data includes: dividing the rectangular cylinders according to the number of observation points, the size of the rectangular cylinders being the distance between the observation point and the adjacent observation point.

[0018] In one embodiment of the present invention, the forward calculation of fitting gravity anomaly to the ice-rock interface elevation model on the undulating aerial gravity observation surface includes:

[0019] The ith observation point P(x i ,zi ) has the following gravity anomaly:

[0020]

[0021] Among them, f i (p) represents the observation point P(x i ,z i ) is the forward modeled ice-rock interface gravity anomaly data, ξ j1 , j2 are the horizontal coordinates of the jth cylinder, p j represents the elevation of the ice-rock interface of the jth glacier column, T j is the glacier surface elevation of the jth glacier column, △g(x i ,z i ,p j ) represents the gravity anomaly caused by the j-th cylinder at the i-th observation point.

[0022] In one embodiment of the present invention, the objective function of the construction regularized inversion specifically includes:

[0023] The objective function of regularized inversion is constructed based on the Tikhonov regularization method:

[0024]

[0025] in, is the data fitting error term; is the model constraint; μ is the regularization parameter with a value range of (0,100]; △g is the observed airborne gravity anomaly data caused by the ice-rock interface obtained by potential field separation; F(p) is the set of fitted gravity anomalies forward calculated at all observation points by the elevation parameters of the ice-rock interface elevation model.

[0026] In one embodiment of the present invention, the model constraint term for:

[0027]

[0028] Among them, p 0 As prior information, the known ice-rock interface elevation information (obtained from geology, earthquakes, drilling, etc.) or the result of the previous inversion can be used. If there is no prior information, p can also be used. 0 is a zero matrix. ; L is a smoothness matrix, which is a (M-2)×M-order matrix, and its specific form is as follows:

[0029]

[0030] In one embodiment of the present invention, the linearization of the nonlinear objective function is to perform Taylor expansion on the objective function line and take the first-order derivative of the objective function ψ(p) as zero.

[0031] In one embodiment of the present invention, by solving the linear equation:

[0032] ψ″(p k )δp k =-ψ′(p k ),

[0033] Get the model correction value δp k .

[0034] In one embodiment of the present invention, the linear equation is solved by singular value decomposition method.

[0035] In one embodiment of the present invention, the iterative formula for the ice-rock interface elevation is:

[0036] p k+1 =p k +δp k , k=0,1,…。

[0037] In one embodiment of the present invention, the correction amount of the ice-rock interface elevation being less than a threshold value comprises:

[0038] Correction δp for ice-rock interface elevation k and the elevation p of the kth iteration k The two-norm ratio of is less than the threshold ε1:

[0039]

[0040] Among them, the threshold ε1 is a preset threshold, and its value range is (0,1].

[0041] Compared with the prior art, a method for inverting glacier thickness using aerial gravity anomalies according to the present invention, based on aerial gravity data, breaks through the idea of ​​traditional spatial domain density interface inversion methods, and newly constructs a glacier thickness calculation method based on aerial gravity, which solves the problems of model segmentation when the interface is undulating and gravity anomaly calculation on the undulating observation surface, so that the high-precision aerial gravity method can be used to measure the thickness of mountain glaciers. BRIEF DESCRIPTION OF THE DRAWINGS

[0042] Figure 1 A schematic diagram of the model division of sedimentary layers using the existing spatial domain interface inversion method;

[0043] Figure 2 A schematic diagram of the dissection of a mountain glacier model provided by the present invention;

[0044] Figure 3A flow chart of a method for inverting glacier thickness using aerial gravity anomaly provided by the present invention;

[0045] Figure 4 A schematic diagram of the steps of a method for inverting glacier thickness using aerial gravity anomaly provided by the present invention;

[0046] Figure 5 A comparison diagram of the ice-rock interface elevation inverted by using aerial gravity anomaly and the result inverted by glacier flow velocity provided by the present invention;

[0047] Figure 6 A comparison diagram of glacier thickness inverted using aerial gravity anomaly and glacier thickness inverted using traditional spatial domain density interface is provided in the present invention. DETAILED DESCRIPTION

[0048] The specific implementation modes of the present invention are described in detail below in conjunction with the accompanying drawings, but it should be understood that the protection scope of the present invention is not limited by the specific implementation modes.

[0049] Unless explicitly stated otherwise, throughout the specification and claims, the term “comprise” or variations such as “include” or “comprising”, etc., will be understood to include the stated elements or components but not to exclude other elements or components.

[0050] like Figures 2 to 6 As shown, a method for inverting glacier thickness using aviation gravity anomaly according to a preferred embodiment of the present invention comprises the following steps:

[0051] S1: Obtain measured aviation gravity data g (gravity from above the ground); wherein the aviation gravity data g is data obtained on an undulating aviation gravity observation surface at a certain height above the ground.

[0052] The specific height of the undulating aerial gravity observation surface from the ground is generally determined according to the actual flight conditions and is not a fixed height. i ,z i ), which is the spatial position of the measured data.

[0053] S2: Perform potential field separation on the aerial gravity data g, remove the background field, and obtain the gravity anomaly data Δg caused by the ice-rock interface.

[0054] The measured aerial gravity data g can be decomposed into the superposition of anomalies caused by the inhomogeneous geological bodies in the underground:

[0055] g=g1+g2+…+△g+…+g n

[0056] Among them, g iis the gravity anomaly data caused by the i-th geological body on the aerial gravity observation surface.

[0057] Therefore, the gravity anomaly data Δg caused by the ice-rock interface can be separated by using the traditional upward continuation method and other potential field separation methods. The gravity anomaly data Δg is the gravity anomaly data g caused by each geological body on the airborne gravity observation surface. i The matrix formed.

[0058] S3: Based on the measured aerial gravity data, the glacier is divided into multiple horizontally juxtaposed rectangular columns, and the initial elevation of the ice-rock interface p is calculated. 0 Set to the glacier surface elevation T.

[0059] like Figure 2 As shown in the figure, the ice layer is adaptively divided into M horizontally juxtaposed rectangular cylinders. The top interface of the ice layer is the undulating terrain, and the bottom interface is the undulating ice-rock interface that needs to be inverted. In order to make full use of the data, the rectangular cylinders are generally divided according to the number of observation points. One observation point corresponds to one divided cylinder, and the size of the cylinder is the distance between the observation point and the adjacent observation point.

[0060] The elevation P is the set of elevations of each rectangular cylinder, P = (p1, p2, ... p j ,…), j=1,2,…,M.

[0061] S4: Perform forward calculation of gravity anomaly fitting to the ice-rock interface elevation model on the undulating aerial gravity observation surface, and calculate the gravity anomaly at each observation point.

[0062] Specifically, Figure 2 As shown in the figure, the forward calculation process of fitting the gravity anomaly to the ice-rock interface elevation model on the undulating airborne gravity observation surface is as follows:

[0063] On the undulating aeronautical gravity observation surface, the i-th observation point P(x i ,z i ) can be approximately expressed as the forward cumulative sum of M horizontally juxtaposed rectangular cylinders:

[0064]

[0065] Among them, f i (p) represents the observation point P(x i ,z i ) is the forward modeled ice-rock interface gravity anomaly data, ξ j1 , j2 is the horizontal coordinate of the jth cylinder, dξ and dζ are the horizontal and vertical integral differentials of the partitioned rectangle, respectively, and p j represents the elevation of the ice-rock interface of the jth glacier column, T jis the glacier surface elevation of the jth glacier column, Δg(x i ,z i ,p j ) represents the gravity anomaly caused by the j-th cylinder at the i-th observation point.

[0066] S5: Based on the Tikhonov regularization method, the objective function of regularized inversion is constructed, and the nonlinear objective function is linearized to obtain the correction value of the ice-rock interface elevation model, and the ice-rock interface elevation p is iteratively updated.

[0067] In this step S5, the objective function ψ(p) of regularized inversion is constructed based on the Tikhonov regularization method, and the ice-rock interface elevation model parameter p that makes the objective function minimal is calculated. The objective function ψ(p) is:

[0068]

[0069] in, is the data fitting error term, is a model constraint term, which is used to introduce additional prior information to constrain the ice-rock interface elevation model parameter p; μ is a regularization parameter, which controls the weight of data fitting error and model constraint deviation in the process of minimizing the objective function, and its value range is generally (0,100]; △g is the observed airborne gravity anomaly data caused by the ice-rock interface obtained by potential field separation in step S2; F(p)≡{f1(p),…,f i (p),…,f N (p),} is the set of fitted gravity anomalies forward calculated at all observation points for the elevation parameter p of the ice-rock interface elevation model in step S4.

[0070] Among them, the model constraints for:

[0071]

[0072] Among them, p 0 As prior information, the known ice-rock interface elevation information (obtained from geology, earthquakes, drilling, etc.) or the result of the previous inversion can be used. If there is no prior information, p can also be used. 0 is a zero matrix. L is a smoothness matrix, which is a (M-2)×M-order matrix, and its specific form is as follows:

[0073]

[0074] In step S5, the objective function line is Taylor expanded, and the first-order derivative of the objective function ψ(p) is taken as zero, that is,

[0075] ψ′(p k )+ψ″(p k )(pp k )=0

[0076] Among them, ψ′(p k ) and ψ″(p k ) are the first-order and second-order derivatives of the objective function ψ(p) at the kth iteration. k The value of .

[0077] δp k =(pp k ),

[0078] Then, the correction value δp of the model parameters is derived k The calculation equation is:

[0079] ψ″(p k )δp k =-ψ′(p k ).

[0080] Solving the above linear equation can obtain the model correction value δp k , and bring the model correction into the iterative formula of the ice-rock interface elevation p. The iterative formula of the ice-rock interface elevation p is:

[0081] p k+1 =p k +δp k , k=0,1,…。

[0082] The above linear equation can be solved by singular value decomposition method, which is an algorithm for solving the equation. There are also many other methods that can be used to solve it.

[0083] S6: Repeat the iterative process from step S4 to S5 until the correction value δp of the ice-rock interface elevation p is k and the elevation p of the kth iteration k When the ratio of the two norms is less than the threshold, the iteration is stopped and the ice-rock interface elevation p is obtained.

[0084] In this step S6, when When , the iteration stops; where ε1 is the preset threshold, and the general value interval is (0,1].

[0085] S7: Calculate the glacier thickness H based on the glacier surface elevation T and the ice-rock interface elevation p.

[0086] In this step S6, the glacier thickness H is H=Tp.

[0087] S8: Output the calculated ice-rock interface elevation p and glacier thickness H.

[0088] The method of inverting glacier thickness using aerial gravity anomaly not only takes into account the undulating glacier surface, but also can be applied to calculate aerial gravity data obtained on an undulating observation surface at a certain height above the ground, thereby improving the reliability of glacier thickness calculation.

[0089] In order to verify the effectiveness and technical advantages of the method, two lateral lines in the Puruogangri Glacier in Qiangtang were selected to perform inversion calculations of glacier thickness based on measured data, and the results were compared with the results of the traditional spatial domain density interface inversion method and the results of glacier flow velocity inversion. Figure 5 and Figure 6 As shown in the figure, it can be seen that this scheme improves the reliability of glacier thickness inversion based on aerial gravity, can obtain a large area and continuous ice-rock interface, and can make up for the shortcomings of ice radar measurement results.

[0090] The foregoing description of specific exemplary embodiments of the present invention is for the purpose of illustration and demonstration. These descriptions are not intended to limit the present invention to the precise form disclosed, and it is clear that many changes and variations can be made based on the above teachings. The purpose of selecting and describing the exemplary embodiments is to explain the specific principles of the present invention and its practical application, so that those skilled in the art can realize and utilize various different exemplary embodiments of the present invention and various different selections and changes. The scope of the present invention is intended to be limited by the claims and their equivalents.

Claims

1. A method for inverting glacier thickness using aerial gravity anomaly, characterized in that: The following steps are involved: Obtain measured aerial gravity data, perform potential field separation on the aerial gravity data, remove the background field, and obtain gravity anomaly data caused by the ice-rock interface; among them, the aerial gravity data obtained at the observation points on the undulating aerial gravity observation surface; According to the measured aerial gravity data, the glacier was divided into multiple horizontally juxtaposed rectangular columns, and the initial elevation of the ice-rock interface was set as the glacier surface elevation. The forward calculation of gravity anomaly is performed by fitting the ice-rock interface elevation model on the undulating aerial gravity observation surface, and the gravity anomaly of each observation point is calculated; The objective function of regularized inversion is constructed, and the nonlinear objective function is linearized to obtain the correction value of the ice-rock interface elevation model, and the ice-rock interface elevation is iteratively updated; Repeat the iterative process until the correction amount of the ice-rock interface elevation is less than the threshold, then stop the iteration and obtain the ice-rock interface elevation; The thickness of the glacier is calculated based on the elevation of the glacier surface and the elevation of the ice-rock interface.

2. The method for inverting glacier thickness using aerial gravity anomaly according to claim 1, characterized in that: The method of dividing the glacier into a plurality of horizontally juxtaposed rectangular columns according to measured aerial gravity data includes: dividing the rectangular columns according to the number of observation points, wherein the size of the rectangular columns is the distance between the observation point and an adjacent observation point.

3. The method for inverting glacier thickness using aerial gravity anomaly according to claim 1, characterized in that: The forward calculation of fitting gravity anomaly to the ice-rock interface elevation model on the undulating aerial gravity observation surface includes: The ith observation point P(x i ,z i ) is as follows: Among them, f i (p) represents the observation point P(x i ,z i ) is the forward modeled ice-rock interface gravity anomaly data, ξ j1 , j2 are the horizontal coordinates of the jth cylinder, p j represents the elevation of the ice-rock interface of the jth glacier column, T j is the glacier surface elevation of the jth glacier column, Δg(x i ,z i ,p j ) represents the gravity anomaly caused by the jth column at the i-th observation point; M is the number of horizontally juxtaposed rectangular columns into which the ice layer is adaptively divided in step S2; z i is the elevation of the ith observation point.

4. The method for inverting glacier thickness using aerial gravity anomaly according to claim 3, characterized in that: The objective function of the structural regularization inversion specifically includes: The objective function of regularized inversion is constructed based on the Tikhonov regularization method: in, is the error term for data fitting; is the model constraint; μ is the regularization parameter with a value range of (0,100]; △g is the observed airborne gravity anomaly data caused by the ice-rock interface obtained by potential field separation; F(p) is the set of fitted gravity anomalies forward calculated at all observation points by the elevation parameters of the ice-rock interface elevation model.

5. The method for inverting glacier thickness using aerial gravity anomaly according to claim 4, characterized in that: The model constraints for: Among them, p 0 As prior information, the known ice-rock interface elevation information or the result of the previous inversion is used; if there is no prior information, p is taken 0 is a zero matrix; L is a smoothness matrix, which is a (M-2)×M-order matrix, and its specific form is as follows:

6. The method for inverting glacier thickness using aerial gravity anomaly according to claim 5, characterized in that: The linearization of the nonlinear objective function is to perform Taylor expansion on the objective function line and take the first-order derivative of the objective function ψ(p) as zero, and obtain ψ′(p k )+ψ″(p k (pp) k )=0 Among them, ψ′(p k ) and ψ″(p k ) are the first-order and second-order derivatives of the objective function ψ(p) at the kth iteration. k The value of .

7. The method for inverting glacier thickness using aerial gravity anomaly according to claim 6, characterized in that: Let δp k =(pp k ), then the model correction value δp is derived k The calculation equation is: ψ″(p k )δp k =-ψ′(p k ), Solving the above linear equation can obtain the model correction value δp k .

8. The method for inverting glacier thickness using aerial gravity anomaly according to claim 7, characterized in that: Solve linear equations via singular value decomposition.

9. The method for inverting glacier thickness using aerial gravity anomaly according to claim 8, characterized in that: The iterative formula for the ice-rock interface elevation is: p k+1 =p k +δp k ,k=0,1,… Among them, p k is the elevation of the kth iteration, p k+1 Elevation at the k+1th iteration.

10. The method for inverting glacier thickness using aerial gravity anomaly according to claim 9, characterized in that: The correction amount of the ice-rock interface elevation being less than a threshold value includes: Correction δp for ice-rock interface elevation k and the elevation p of the kth iteration k The two-norm ratio of is less than the threshold ε1: Among them, the threshold ε1 is a preset threshold, and its value range is (0,1].

Images