Well-ground-air multi-dimensional magnetic anomaly joint inversion method and system

Abstract

The application discloses a kind of well-ground empty multi-dimensional magnetic anomaly joint inversion method and system, system includes data acquisition module, grid division module, matrix generation module, matrix combination module, equation construction module and equation solution module.First, the total magnetization direction of inversion region, geomagnetic direction and observation data are acquired by data acquisition module.Then, grid division module carries out cuboid grid division to underground space.Afterwards, matrix generation module calculates ground data kernel matrix, aviation data kernel matrix, well three-component data kernel matrix and pre-optimization matrix.Then, matrix combination module combines all well three-component data kernel matrix, ground data kernel matrix and aviation data kernel matrix into total kernel matrix, and observation data is combined into total observation data matrix.Afterwards, equation construction module constructs joint inversion target equation based on the foregoing matrix.Last, solve joint inversion target equation by equation solution module.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of seismic data processing, in particular to a well-ground-space multi-dimensional magnetic anomaly joint inversion method and system. BACKGROUND

[0002] According to the different environments during data acquisition, magnetic exploration can be mainly divided into three categories of well, ground and space, wherein the underground magnetic anomaly body information contained in the three categories of well, ground and space magnetic survey data is different, and each has advantages and disadvantages. For ground magnetic survey, the survey network is close to the shallow magnetic anomaly body, which can better identify the shallow detail information, and has high horizontal resolution, but the vertical resolution is seriously insufficient. For aerial magnetic survey, the survey network is far away from the shallow field source and has a larger control range, which can better identify the deep large-scale field source information, but has weak identification ability for shallow small-scale field source information. For three-component magnetic survey in the well, the survey points are distributed along the well axis in the vertical direction and are closer to the deep magnetic anomaly body, which has high vertical resolution, but also has the problem of low horizontal resolution. Since the underground magnetic anomaly body often has different distribution positions and sizes, the conventional single magnetic survey data inversion has the problems of serious inversion result multi-solution and insufficient resolution, and it is difficult to accurately restore the physical property distribution of all ore bodies at the same time. SUMMARY

[0003] In view of the deficiencies of the prior art, the present application provides a well-ground-space multi-dimensional magnetic anomaly joint inversion method and system. The inversion result can have high horizontal resolution and high vertical resolution at the same time, and can identify shallow small-scale information and deep large-scale information at the same time. The specific technical solutions are as follows:

[0004] In a first aspect, a well-ground-space multi-dimensional magnetic anomaly joint inversion method is provided, comprising:

[0005] Obtaining the total magnetization direction, the ground magnetization direction and the observation data of the inversion region, wherein the observation data includes ground observation data, aerial observation data and three-component observation data of multiple drill holes;

[0006] Performing cuboid grid division on the underground space of the inversion region, determining the center coordinates and the corresponding edge length of each cuboid grid unit, and determining the distance parameters including the distances between the center of each grid unit and the drill hole axis, the ground survey network and the aerial survey network;

[0007] According to the total magnetization direction, the ground magnetization direction, the center coordinates and the corresponding edge length of each cuboid grid unit, calculating the ground data kernel matrix, the aerial data kernel matrix corresponding to different measurement heights, and the three-component data kernel matrix corresponding to different drill holes, and calculating the pre-optimization factor corresponding to each grid unit according to the corresponding distance parameters to generate a pre-optimization matrix;

[0008] combining the ground data kernel matrix, all the aerial data kernel matrices and all the borehole three-component data kernel matrices into a total kernel matrix, and combining the observation data into a total observation data matrix;

[0009] establishing a joint inversion objective equation based on the pre-optimization matrix, the total observation data matrix and the total kernel matrix;

[0010] solving the joint inversion objective equation by using a pre-optimization conjugate gradient algorithm to obtain the magnetization intensity distribution of the subsurface space.

[0011] In combination with the first aspect, in a first implementation manner of the first aspect, the ground data kernel matrix, the aerial data kernel matrices corresponding to different measurement heights and the three-component data kernel matrices corresponding to different boreholes are determined, and the determination comprises:

[0012] determining the three-component data kernel matrix according to the total magnetization direction, the coordinates of the center of each cuboid grid unit and the corresponding side length, and the coordinates of each measuring point;

[0013] determining the ground data kernel matrix and the aerial data kernel matrix by using the borehole three-component data kernel matrix and the geomagnetic direction.

[0014] In combination with the first aspect, in a second implementation manner of the first aspect, all the borehole three-component data kernel matrices, the ground data kernel matrix and all the aerial data kernel matrices are multiplied by corresponding weight factors respectively and combined into the total kernel matrix.

[0015] In combination with the first aspect, in a third implementation manner of the first aspect, the obtained observation data are multiplied by corresponding weight factors respectively, and combined into the total observation data matrix according to the combination order of the total kernel matrix.

[0016] In combination with the first aspect, in a fourth implementation manner of the first aspect, the pre-optimization factors corresponding to each grid unit are combined into a pre-optimization matrix in a diagonal matrix form according to the arrangement order of the grid units.

[0017] Secondly, a well-ground-air multi-dimensional magnetic anomaly joint inversion system is provided, comprising:

[0018] a data acquisition module configured to acquire a total magnetization direction of an inversion region, a geomagnetic direction and observation data, the observation data comprising ground observation data, aerial observation data and borehole three-component observation data of a plurality of boreholes;

[0019] a grid division module configured to perform cuboid grid division on a subsurface space of the inversion region, to determine the coordinates of the center of each cuboid grid unit and the corresponding side length, and to determine distance parameters comprising distances between the center of the grid unit and the axis of the borehole, the ground measuring network and the aerial measuring network;

[0020] a matrix generating module, configured to generate a ground data kernel matrix, an aerial data kernel matrix corresponding to different measuring heights, and a borehole three-component data kernel matrix corresponding to different boreholes according to the total magnetization direction, the ground magnetization direction, the coordinates of the center of each cuboid grid unit, and the corresponding edge length, and to generate a preconditioning matrix by calculating a preconditioning factor corresponding to each grid unit according to a corresponding distance parameter;

[0021] a matrix combining module, configured to combine the ground data kernel matrix, all the aerial data kernel matrices, and all the borehole three-component data kernel matrices into a total kernel matrix, and to combine the observation data into a total observation data matrix;

[0022] an equation establishing module, configured to establish a joint inversion target equation based on the preconditioning matrix, the total observation data matrix, and the total kernel matrix;

[0023] an equation solving module, configured to solve the joint inversion target equation by using a preconditioning conjugate gradient algorithm to obtain the magnetization intensity distribution of the subsurface space.

[0024] With reference to the second aspect, in a first implementation manner of the second aspect, the matrix generating module comprises:

[0025] a borehole three-component data matrix generating unit, configured to determine a three-component data kernel matrix according to the total magnetization direction, the coordinates of the center of each cuboid grid unit, and the corresponding edge length, and the coordinates of each measuring point;

[0026] an air-ground data matrix generating unit, configured to determine the ground data kernel matrix and the aerial data kernel matrix by using the borehole three-component data kernel matrix and the ground magnetization direction.

[0027] With reference to the second aspect, in a second implementation manner of the second aspect, the matrix combining module multiplies all the borehole three-component data kernel matrices, the ground data kernel matrix, and all the aerial data kernel matrices by corresponding weight factors respectively and combines them into the total kernel matrix.

[0028] With reference to the second aspect, in a third implementation manner of the second aspect, the matrix combining module multiplies the obtained observation data by corresponding weight factors respectively, and combines them into the total observation data matrix according to the combination order of the total kernel matrix.

[0029] With reference to the second aspect, in a fourth implementation manner of the second aspect, the matrix generating module combines the preconditioning factors corresponding to each grid unit into the preconditioning matrix in the form of a diagonal matrix according to the arrangement order of the grid units.

[0030] Beneficial effects: the ground-air multi-dimensional magnetic anomaly joint inversion method and system can flexibly combine different types of magnetic anomaly data to realize joint inversion, fully utilize different target information contained in different types of magnetic anomaly data, make the inversion result have high horizontal resolution and high vertical resolution, and can simultaneously identify shallow small-scale information and deep large-scale information. BRIEF DESCRIPTION OF DRAWINGS

[0031] In order to more clearly illustrate the specific embodiments of the present application, the drawings required to be used in the specific embodiments will be briefly introduced below. In all the drawings, the elements or parts are not necessarily drawn according to the actual proportion.

[0032] Figure 1 The flow chart of the well-ground-air multi-dimensional magnetic anomaly joint inversion method provided by an embodiment of the present application;

[0033] Figure 2 The system block diagram of the well-ground-air multi-dimensional magnetic anomaly joint inversion system provided by an embodiment of the present application. DETAILED DESCRIPTION

[0034] The embodiments of the technical solutions of the present application will be described in detail below with reference to the drawings. The following embodiments are only used to more clearly illustrate the technical solutions of the present application, and therefore only serve as examples, but cannot limit the protection scope of the present application.

[0035] As shown in the flow chart of the well-ground-air multi-dimensional magnetic anomaly joint inversion method, the inversion method comprises: Figure 1

[0036] Step 1, obtaining the total magnetization direction, the geomagnetization direction and the observation data of the inversion area, the observation data including the ground observation data, the aerial observation data and the three-component observation data of multiple drill holes;

[0037] Step 2, carrying out cuboid grid division on the underground space of the inversion area, determining the center coordinates and the corresponding side length of each cuboid grid unit, and determining the distance parameters including the distances between the center of the grid unit and the drill hole axis, the ground survey network and the aerial survey network;

[0038] Step 3, according to the total magnetization direction, the geomagnetization direction, the center coordinates and the corresponding side length of each cuboid grid unit, determining the ground data kernel matrix, the aerial data kernel matrix corresponding to different measurement heights, and the three-component data kernel matrix corresponding to different drill holes, and calculating the pre-optimization factor corresponding to each grid unit according to the corresponding distance parameters to generate a pre-optimization matrix;

[0039] Step 4, combining the ground data kernel matrix, all aerial data kernel matrices and all three-component data kernel matrices into a total kernel matrix, and combining the observation data into a total observation data matrix;​

[0040] Step 5, establishing a joint inversion objective equation based on the pre-optimized matrix, the total observation data matrix and the total kernel matrix;

[0041] Step 6, solving the joint inversion objective equation by using a pre-optimized conjugate gradient algorithm to obtain the magnetization intensity distribution of the underground space.

[0042] Specifically, first, the total magnetization direction, the geomagnetization direction and the observation data of the inversion region can be obtained. The total magnetization direction includes a total magnetization inclination I and a total magnetization declination D, and the geomagnetization direction includes a geomagnetization inclination I0 and a geomagnetization declination D0. The observation data includes ground observation data obtained by using magnetic exploration, and selected aerial observation data at different measurement heights and selected borehole three-component observation data at different boreholes.

[0043] Then, the underground space of the inversion region can be cuboid grid divided to determine the center coordinates and the corresponding edge lengths of each cuboid grid cell. The determined distance parameters include the horizontal distance between the grid cell center and the borehole axis, the vertical distance between the grid cell center and the ground survey network, and the vertical distance between the grid cell center and each layer of the aerial survey network.

[0044] The horizontal distance between the grid cell center and the borehole axis, the vertical distance between the grid cell center and the ground survey network, and the vertical distance between the grid cell center and the aerial survey network can be calculated by the coordinates of the grid cell center and the coordinates of the borehole axis, the coordinates of the grid cell center and the coordinates of the ground survey network, and the coordinates of the grid cell center and the coordinates of the aerial survey network, respectively.

[0045] Then, the ground data kernel matrix, the aerial data kernel matrix corresponding to different measurement heights and the borehole three-component data kernel matrix corresponding to different boreholes can be determined according to the total magnetization direction, the geomagnetization direction and the center coordinates and the corresponding edge lengths of each cuboid grid cell. In this embodiment, the three-component data kernel matrix can be first determined according to the total magnetization direction, the center coordinates and the corresponding edge lengths of each cuboid grid cell, and the coordinates of each survey point. The specific calculation formula is as follows:

[0046]

[0047]

[0048] Wherein, μ0 is the magnetic permeability of vacuum, usually 4π×10 -7 , (x0, y0, z0) is the center coordinates of the cuboid grid cell, a, b, c are the three edge lengths of the cuboid grid cell, (ξ, η, ζ) is the coordinates of a volume element in the grid cell, and R is the distance from the volume element to the survey point (x, y, z).

[0049] The ground data kernel matrix and the airborne data kernel matrix can be determined by the borehole three-component data kernel matrix and the direction of geomagnetization. The specific calculation formula is as follows:

[0050] g sur = g air = g hax cosI0cosD0+g hay cosI0sinD0+g za sinI0.

[0051] In the embodiment, the pre-optimization factors corresponding to each grid cell can be calculated according to the corresponding distance parameters, and a pre-optimization matrix is generated. Specifically, the pre-optimization factors of each grid cell are calculated according to the horizontal distance between the center of the grid cell and the borehole axis, the vertical distance between the center of the grid cell and the ground survey network, and the vertical distance between the center of the grid cell and each layer of the airborne survey network. The specific calculation formula is as follows:

[0052]

[0053]

[0054] P j denotes the pre-optimization factor corresponding to the jth grid cell; ND is the number of boreholes, K denotes the borehole number, ΔL j is the horizontal distance from the center of the jth grid cell to the borehole axis; ΔH j is the vertical distance from the center of the jth grid cell to the ground survey network; NA is the number of airborne survey networks, k denotes the survey network number, and Δh j is the vertical distance from the center of the jth grid cell to the airborne survey network. β is related to the magnetic anomaly decay rate, generally taken as 2 in two-dimensional inversion and 3 in three-dimensional inversion, and C is a constant used to adjust the value of the pre-optimization factor, generally taken as 1.

[0055] After calculating the pre-optimization factors of each grid cell, the pre-optimization factors corresponding to each grid cell can be combined into a diagonal matrix according to the grid cell arrangement order, and the diagonal matrix obtained by combination is the pre-optimization matrix. The specific form of the pre-optimization matrix is as follows:

[0056]

[0057] Then, all the borehole three-component data kernel matrices, the ground data kernel matrix, and all the airborne data kernel matrices can be multiplied by the corresponding weight factors and combined into the total kernel matrix. The total kernel matrix obtained by combination is as follows:

[0058]

[0059] wherein g sur denotes the ground data kernel matrix, and w​sur g represents the weighting factor for ground data. air1 ...g airn w represents the kernel matrix of aerial data for measuring altitude at n levels. air1 ...w airn g represents the weighting factor for aerial data measuring altitude at n levels. hax1 ...g haxn w represents the kernel matrix of n borehole Hax components. hax1 ...w haxn G represents the weighting factor of the n borehole hax components. hay1 ...g hayn w represents the kernel matrix of n borehole hay components. hay1 ...w hayn G represents the weighting factor of the n borehole hay components; za1 ...g zan w represents the kernel matrix of the z-components of n boreholes. za1 ...w zan This represents the weighting factor of the z-components of the n boreholes.

[0060] The acquired observation data are multiplied by their respective weighting factors and combined in the order of combination with the total kernel matrix to form the total observation data matrix, which is:

[0061]

[0062] in, Represents ground observation data, This represents aerial observation data at the nth floor measurement altitude. This represents the observation data of the hax components of n boreholes. This represents the observation data of the hay components of n boreholes. This represents the observation data of the z-components of n boreholes.

[0063] In practical applications, combinations can be made based on the available data types.

[0064] Then, a joint inversion objective equation can be established based on the total observation data matrix d, the total kernel matrix G, and the pre-optimization matrix P. The joint inversion objective equation is as follows:

[0065] PG T Gm = PG T d;

[0066] Where m is the parameter to be solved, which characterizes the distribution of the predicted magnetization intensity in the underground space.

[0067] Finally, the joint inversion target equation is solved by using the pre-optimization conjugate gradient algorithm, so that the predicted magnetization intensity distribution of the underground space is obtained.

[0068] As shown in the system block diagram of the well-ground-space multi-dimensional magnetic anomaly joint inversion system, the inversion system comprises: Figure 2

[0069] The data acquisition module is configured to acquire the total magnetization direction, the ground magnetization direction and the observation data of the inversion region, wherein the observation data comprises ground observation data, aerial observation data and three-component observation data of multiple drill holes;

[0070] The grid division module is configured to perform cuboid grid division on the underground space of the inversion region, determine the center coordinates and corresponding edge lengths of each cuboid grid unit, and determine distance parameters including distances between the centers of the grid units and the drill hole axes, the ground survey network and the aerial survey network;

[0071] The matrix generation module is configured to generate a ground data kernel matrix, an aerial data kernel matrix corresponding to different measurement heights, and a three-component data kernel matrix corresponding to different drill holes according to the total magnetization direction, the ground magnetization direction, the center coordinates and corresponding edge lengths of each cuboid grid unit, and calculate pre-optimization factors corresponding to each grid unit according to the corresponding distance parameters to generate a pre-optimization matrix;

[0072] The matrix combination module is configured to combine the ground data kernel matrix, all aerial data kernel matrices and all three-component data kernel matrices into a total kernel matrix, and combine the observation data into a total observation data matrix;

[0073] The equation construction module is configured to establish a joint inversion target equation based on the pre-optimization matrix, the total observation data matrix and the total kernel matrix;

[0074] The equation solving module is configured to solve the joint inversion target equation by using a pre-optimization conjugate gradient algorithm to obtain the magnetization intensity distribution of the underground space.

[0075] Specifically, the inversion system is composed of a data acquisition module, a grid division module, a matrix generation module, a matrix combination module, an equation construction module and an equation solving module.

[0076] The data acquisition module can acquire the total magnetization direction, the ground magnetization direction and the observation data of the inversion region. The grid division module can perform cuboid grid division on the underground space of the inversion region, determine the center coordinates and corresponding edge lengths of each cuboid grid unit, and determine distance parameters including distances between the centers of the grid units and the drill hole axes, the ground survey network and the aerial survey network.

[0077] ​The matrix generating module can determine the ground data kernel matrix, the aerial data kernel matrix corresponding to different height layers and the borehole three-component data kernel matrix corresponding to different boreholes according to the total magnetization direction, the geomagnetization direction, the center coordinates and the corresponding side length of each cuboid grid unit.

[0078] In the embodiment, the matrix generating module comprises a borehole three-component data matrix generating unit, a ground-aerial data matrix generating unit and a pre-optimization matrix generating unit. The borehole three-component data matrix generating unit can first determine the borehole three-component data kernel matrix according to the total magnetization direction, the center coordinates and the corresponding side length of each cuboid grid unit and the coordinates of each measuring point. The ground-aerial data matrix generating unit can determine the ground data kernel matrix and the aerial data kernel matrix according to the borehole three-component data kernel matrix and the geomagnetization direction.

[0079] The pre-optimization matrix generating unit can calculate the pre-optimization factor of each grid unit according to the horizontal distance between the center of the grid unit and the borehole axis of each borehole, the vertical distance between the center of the grid unit and the ground measuring network and the vertical distance between the center of the grid unit and each layer of the aerial measuring network. The pre-optimization factor corresponding to each grid unit is combined into a diagonal matrix according to the arrangement order of the grid units, and the diagonal matrix obtained by the combination is the pre-optimization matrix.

[0080] The matrix combination module can multiply all the borehole three-component data kernel matrices, the ground data kernel matrix and all the aerial data kernel matrices by the corresponding weight factors respectively and combine them into the total kernel matrix. The obtained observation data are multiplied by the corresponding weight factors respectively, and combined into the total observation data matrix according to the combination order of the total kernel matrix. In actual application, the matrix combination module can combine the data types owned by it.

[0081] The equation construction module can establish a joint inversion target equation based on the total observation data matrix d, the total kernel matrix G and the pre-optimization matrix P. The equation solving module can solve the joint inversion target equation by using the pre-optimization conjugate gradient algorithm, and the predicted magnetization intensity distribution of the underground space can be obtained.

[0082] The above embodiments are only used to illustrate the technical solutions of the present application, but not limit the same. Although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that they can modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacements to some or all of the technical features; and these modifications or replacements do not make the essence of the corresponding technical solutions deviate from the scope of the technical solutions of the embodiments of the present application, and they should be covered in the scope of the claims and the specification of the present application.

Claims

1. A well-ground multi-dimensional magnetic anomaly joint inversion method, characterized in that, The method comprises the following steps: Obtaining the total magnetization direction, the geomagnetization direction and observation data of an inversion region, wherein the observation data comprises ground observation data, aerial observation data and three-component observation data of multiple drill holes; Carrying out cuboid grid division on the underground space of the inversion region, determining the center coordinates and corresponding side lengths of each cuboid grid unit, and determining distance parameters including the distances between the center of each grid unit and the drill hole axis, the ground survey network and the aerial survey network; According to the total magnetization direction, the geomagnetization direction, the center coordinates and corresponding side lengths of each cuboid grid unit, calculating the ground data kernel matrix, the aerial data kernel matrix corresponding to different measurement heights and the three-component data kernel matrix corresponding to different drill holes, and calculating the preconditioning factor corresponding to each grid unit according to the corresponding distance parameters to generate a preconditioning matrix; The specific calculation formula of the preconditioning factor is as follows: ； in, For the first The pre-optimization factor corresponding to each grid cell For the number of holes, Number the boreholes. For the first The horizontal distance from the center of each grid cell to the borehole axis For the first The vertical distance from the center of each grid cell to the ground survey network. For the number of aerial survey networks, For the measurement network number, For the first The vertical distance from the center of each grid cell to the airborne survey network. It is related to the decay rate of the magnetic anomaly, taking a value of 2 in two-dimensional inversion and 3 in three-dimensional inversion. It is a constant; Combining the preconditioning factor corresponding to each grid unit into a diagonal matrix according to the grid unit arrangement order, and the diagonal matrix obtained by the combination is the preconditioning matrix; Combining the ground data kernel matrix, all aerial data kernel matrices and all three-component data kernel matrices into a total kernel matrix, and combining the observation data into a total observation data matrix; Based on the preconditioning matrix, the total observation data matrix and the total kernel matrix, a joint inversion objective equation is established; The joint inversion objective equation is solved by using a preconditioning conjugate gradient algorithm to obtain the magnetization intensity distribution of the underground space.

2. The well-ground multi-dimensional magnetic anomaly joint inversion method according to claim 1, characterized in that, Determining the ground data kernel matrix, the aerial data kernel matrix and the three-component data kernel matrix comprises: Determining the three-component data kernel matrix of the drill hole according to the total magnetization direction, the center coordinates and corresponding side lengths of each cuboid grid unit and the coordinates of each measuring point; Determining the ground data kernel matrix and the aerial data kernel matrix through the three-component data kernel matrix of the drill hole and the geomagnetization direction.

3. The well-ground multi-dimensional magnetic anomaly joint inversion method according to claim 1, characterized in that, Multiplying all the three-component data kernel matrices, the ground data kernel matrix and all the aerial data kernel matrices by corresponding weight factors respectively and combining them into the total kernel matrix.

4. The well-ground multi-dimensional magnetic anomaly joint inversion method according to claim 1, characterized in that, Multiplying the obtained observation data by corresponding weight factors respectively, and combining them into the total observation data matrix according to the combination order of the total kernel matrix.

5. A well-ground multi-dimensional magnetic anomaly joint inversion system characterized in that, The method comprises the following steps: A data acquisition module is configured to obtain the total magnetization direction, the geomagnetization direction and observation data of an inversion region, wherein the observation data comprises ground observation data, aerial observation data and three-component observation data of multiple drill holes; A grid division module is configured to carry out cuboid grid division on the underground space of the inversion region, determine the center coordinates and corresponding side lengths of each cuboid grid unit, and determine distance parameters including the distances between the center of each grid unit and the drill hole axis, the ground survey network and the aerial survey network; A matrix generation module is configured to generate a ground data kernel matrix, an aerial data kernel matrix corresponding to different measurement heights and a three-component data kernel matrix corresponding to different drill holes according to the total magnetization direction, the geomagnetization direction, the center coordinates and corresponding side lengths of each cuboid grid unit, and generate a preconditioning matrix by calculating the preconditioning factor corresponding to each grid unit according to the corresponding distance parameters, wherein the specific calculation formula of the preconditioning factor is as follows: ； wherein, is the pre-optimization factor corresponding to the th grid cell, is the number of drill holes, is the drill hole number, is the pre-optimization factor corresponding to the th grid cell, is the vertical distance from the center of the th grid cell to the ground survey network, is the number of aerial survey networks, is the survey network number, is the vertical distance from the center of the th grid cell to the aerial survey network, is related to the magnetic anomaly decay velocity, which is taken as 2 in two-dimensional inversion and 3 in three-dimensional inversion, is a constant; The pre-optimization factors corresponding to each grid unit are combined into a diagonal matrix according to the grid unit arrangement order, and the diagonal matrix obtained by the combination is the pre-optimization matrix; The matrix combination module is configured to combine the ground data kernel matrix, all the aviation data kernel matrices and all the borehole three-component data kernel matrices into a total kernel matrix, and combine the observation data into a total observation data matrix; The equation construction module is configured to construct a joint inversion target equation based on the pre-optimization matrix, the total observation data matrix and the total kernel matrix; The equation solving module is configured to solve the joint inversion target equation by using a pre-optimization conjugate gradient algorithm to obtain the magnetic intensity size distribution of the underground space.

6. The well-ground multi-dimensional magnetic anomaly joint inversion system according to claim 5, characterized in that, The matrix generation module comprises: The borehole three-component data matrix generation unit is configured to determine the three-component data kernel matrix according to the total magnetization direction, the coordinates of the center of each cuboid grid unit and the corresponding edge length, and the coordinates of each measuring point; The ground-air data matrix generation unit is configured to determine the ground data kernel matrix and the aviation data kernel matrix by using the borehole three-component data kernel matrix and the ground magnetization direction.

7. The well-ground multi-dimensional magnetic anomaly joint inversion system according to claim 5, characterized in that, The matrix combination module multiplies all the borehole three-component data kernel matrices, the ground data kernel matrix and all the aviation data kernel matrices by corresponding weight factors respectively and combines them into the total kernel matrix.

8. The well-ground multi-dimensional magnetic anomaly joint inversion system according to claim 5, characterized in that, The matrix combination module multiplies the obtained observation data by corresponding weight factors respectively and combines them into the total observation data matrix according to the combination order of the total kernel matrix.

Images