A fracture identification method based on airborne gravity and magnetic data
By combining the Euler deconvolution method and neural network algorithm, combined with fuzzy clustering and particle swarm optimization, the problem of fault identification accuracy under complex geological conditions was solved, and more accurate fault location and scale identification was achieved.
Patent Information
- Application Number
- CN202510399781.4
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- Filing Date
- 2025-04-01
- Publication Date
- 2025-09-12
- Estimated Expiration
- 2045-04-01
AI Technical Summary
Existing technologies have difficulty accurately identifying the location, depth and scale of faults under complex geological conditions. The Euler deconvolution method lacks accuracy and stability in processing gravity and magnetic data in complex geological environments, making it difficult to meet the needs of rapid and accurate identification in large areas.
The Euler deconvolution method is combined to process historical airborne gravity and magnetic data, a neural network algorithm mapping model is established, and a preset fault location prediction model and fuzzy clustering are used. The fault characteristic parameters are optimized through particle swarm intelligence, and a variety of technical means are integrated to identify faults.
The accuracy and reliability of fault identification have been improved, and key information such as the location, depth and scale of the fault can be determined more accurately.
Smart Images

Figure CN120335036B_ABST
Abstract
Description
Technical Field
[0001] The present invention relates to the technical field of fracture identification, and more particularly to a fracture identification method based on aerial gravity and magnetic data. Background Art
[0002] In the fields of Earth science and resource exploration, faults, as key components of geological structures, profoundly influence various geological processes and resource distribution patterns. From vast oil and gas fields on land to deeply buried mineral deposits and even the crustal plates of the ocean, faults not only maintain the stability of geological structures, but their formation and distribution also directly determine the formation and migration paths of oil, gas, and mineral resources. Therefore, accurately identifying the location, morphology, and scale of faults is of immeasurable significance for in-depth research on geological evolution, efficient exploration of potential resources, and effective assessment of geological hazard risks.
[0003] While traditional fault identification methods, such as geological sampling and shallow-sediment profile exploration, can provide relatively detailed geological information in local areas, they generally suffer from inefficiency, high costs, and limited coverage. With the continuous expansion of exploration, traditional methods are no longer able to meet the urgent need for rapid and accurate fault identification over large areas. With the continuous advancement of geophysical exploration technology, airborne geophysical surveying, with its significant advantages of high efficiency and wide coverage, has gradually become a key technical tool in geological exploration. Airborne gravity and magnetic surveying can rapidly acquire Earth's gravity and magnetic field data over large areas. Different geological bodies, due to differences in density and magnetism, produce unique anomalies in gravity and magnetic data. The presence of a fault significantly alters the density and magnetism of the surrounding rocks, resulting in corresponding anomalies in airborne gravity and magnetic data. However, relying solely on these gravity and magnetic anomaly data makes it difficult to accurately determine key parameters such as fault location, depth, and scale, and the accuracy and reliability of fault identification need to be further improved.
[0004] Euler deconvolution, a mature potential field inversion technique, has been widely used in geological exploration. Based on the Euler homogeneous equation, this method rapidly inverts and calculates the location and depth of the field source by utilizing potential field anomalies and the "tectonic index" of geological bodies. In the practical application of fault identification, by rationally setting the structural index associated with the fault and combining it with airborne gravity and magnetic anomaly data, the Euler deconvolution method provides a new technical approach for determining the location and depth of the fault, significantly improving the efficiency and accuracy of fault identification.
[0005] Although a large number of studies have applied various gravity and magnetic anomaly inversion technologies to geological structure detection, in the field of fault identification, the accuracy and stability of existing inversion methods in processing gravity and magnetic data are obviously insufficient when faced with complex and changeable geological environments, making it difficult to cope with fault identification tasks under various complex geological conditions.
[0006] Therefore, how to accurately identify faults under complex geological conditions is an urgent problem that needs to be solved by technicians in this field. Summary of the Invention
[0007] In view of this, the present invention provides a fracture identification method based on airborne gravity and magnetic data to solve the problems existing in the above-mentioned background technology.
[0008] In order to achieve the above object, the present invention adopts the following technical solutions:
[0009] A fault identification method based on airborne gravity and magnetic data, comprising:
[0010] Acquire historical airborne gravity and magnetic data, process the data through Euler deconvolution, and obtain historical field source location and depth data; establish a mapping model between the historical airborne gravity and magnetic data and fault characteristic parameters, field source location, and depth data based on a neural network algorithm;
[0011] Based on the historical source location and depth data in the Euler deconvolution process, the preset fracture location prediction model is used to predict the fracture location and obtain the fracture location data; at the same time, the fracture type data is obtained based on fuzzy clustering;
[0012] Based on the area to be measured, obtaining the airborne gravity and magnetic data to be measured; confirming the initial estimated fracture parameters of the area to be measured, calculating the estimated airborne gravity and magnetic data based on the mapping model, calculating the loss function between the airborne gravity and magnetic data to be measured and the estimated airborne gravity and magnetic data, and determining whether the loss function value exceeds a preset accuracy threshold; if it does not exceed the preset accuracy threshold, optimizing the fracture characteristic parameters through particle swarm intelligence; if so, ending the iteration and outputting the optimal fracture characteristic parameters;
[0013] Outputs the fracture location, fracture type, and fracture parameters output in the final iteration step.
[0014] Preferably, the historical airborne gravity and magnetic data are obtained by gravity and magnetic measurements at two different measurement altitudes, namely the first historical airborne gravity and magnetic data Q1 and the second historical airborne gravity and magnetic data Q2.
[0015] Preferably, the fracture characteristic parameters include fracture length, fracture width, fracture inclination and fracture type.
[0016] Preferably, the mapping model specifically includes:
[0017]
[0018] Where l is the fracture length, w is the fracture width, d is the fracture depth, θ is the fracture inclination, F cis the fault type, t1 and t2 are the measured heights, G1 and G2 are the historical gravity and magnetic data, and (x, y) is the field source position.
[0019] Preferably, the fracture position prediction specifically includes:
[0020] Each particle i represents a set of SVM parameter combinations p i =(γ i ,C i ), where γ i is the kernel parameter, C i is the penalty factor;
[0021]
[0022] Where k is the number of iterations, is the velocity vector of the i-th particle at the k+1-th iteration, corresponding to the velocity components of the kernel parameter γ and the penalty factor C, ω is the inertia weight, c1 and c2 are learning factors, r1 and r2 are random numbers between [0,1], pBest i,γ and pBest i,C is the kernel parameter and penalty factor corresponding to the optimal fitness value achieved by the i-th particle in its own history, gBest γ and gBest C is the global optimal kernel parameter and penalty factor corresponding to the optimal fitness value of all particles in the historical iteration, and is the kernel parameter and penalty factor of the i-th particle at the k-th iteration, m is the number of verification samples, is the SVM parameter combination (γ i ,C i ) The result of predicting the j-th validation sample, y j is the actual fracture position label of the jth verification sample;
[0023] SVM model:
[0024]
[0025] 0≤α i ≤C,i=1,…,n;
[0026] K(x i ,x j )=exp(-γ‖x i -x j ‖ 2 );
[0027]
[0028] Among them, xi 、x j , x are the input sample data, which are gravity gradient, magnetic gradient and depth data respectively; γ is the kernel parameter, α i With α j are all Lagrange multipliers, y i is the category label of the sample, b is the bias term;
[0029] Initialize the particle swarm and randomly generate a set of initial positions of particles and speed
[0030] For each particle i, use its SVM parameter combination to train the SVM model and solve the Lagrange multiplier α i , α j And bias term b, calculate the fitness value Fitness i ;
[0031] Update pBest of each particle i,γ 、pBest i,C and gBest γ 、gBest C ;
[0032] Update the particle's speed and position according to the particle speed and position update formula;
[0033] Repeat the iterative update until the maximum number of iterations is reached to obtain the final fracture location prediction model and perform fracture location prediction.
[0034] Preferably, obtaining fracture type data based on fuzzy clustering specifically includes:
[0035] Determine the number of clusters c of the fracture type and select the fuzzy factor m, m∈(1,+∞), set the iteration stop threshold ε; construct the n×c membership matrix U=[u ij ],u ij For sample q i The degree of membership to the jth class,
[0036] Calculate the cluster center v of each class according to the membership matrix U j , the calculation formula is:
[0037]
[0038] Among them, v j is the cluster center of the jth class, which is obtained by weighted averaging the sample data, with the weight determined by the mth power of the membership degree;
[0039] Using the current cluster center v j Update the membership matrix U, and the update formula is:
[0040]
[0041] Among them, ||q i -v j || represents sample x i With cluster center v j distance;
[0042] Define the objective function J of fuzzy clustering, the formula is:
[0043]
[0044] Calculate the objective function value J of this iteration (t) and the objective function value J of the previous iteration (t-1) , if J is satisfied (t) -J (t -1) ≤ε, the iteration is considered to have converged and the iteration is stopped; otherwise, the iteration is continued, and the cluster center and membership matrix are continuously updated until the termination condition is met; when the iteration is completed, the fracture type of each sample is determined according to the final membership matrix U, and for each sample q i , the fracture type is the category with the largest membership.
[0045] Preferably, the loss function specifically includes:
[0046]
[0047] Among them, Loss is the loss function value, n represents the number of data samples, and are the measured gravity data and measured magnetic data of the i-th sample, and are the historical gravity data and historical magnetic data of the i-th sample, respectively, which are the data calculated based on the established mapping model. The historical gravity and magnetic data include historical gravity data and historical magnetic data.
[0048] As can be seen from the above technical solutions, compared with existing technologies, the present invention provides a method for fault identification based on airborne gravity and magnetic data. This method combines the Euler deconvolution method with historical airborne gravity and magnetic data to obtain source location and depth data. A neural network algorithm is then used to establish a mapping model, taking into account the relationship between fracture characteristic parameters and gravity and magnetic data. Furthermore, a preset fracture location prediction model and fuzzy clustering are used to determine the fracture location and type, respectively. By integrating multiple technical approaches, fault information is analyzed and mined from different perspectives. Compared to single methods, this method significantly improves the accuracy and reliability of fault identification, enabling more precise determination of key information such as the fracture's location, depth, and scale. BRIEF DESCRIPTION OF THE DRAWINGS
[0049] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are merely embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on the provided drawings without paying any creative work.
[0050] Figure 1 A diagram of the steps of the method provided by the present invention. DETAILED DESCRIPTION
[0051] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.
[0052] The embodiment of the present invention discloses a fracture identification method based on aerial gravity and magnetic data. Figure 1 Shown, including:
[0053] Acquire historical airborne gravity and magnetic data, process them using Euler deconvolution, and obtain historical field source location and depth data. Establish a mapping model between historical airborne gravity and magnetic data and fault characteristic parameters, field source location, and depth data based on a neural network algorithm.
[0054] Based on the historical source location and depth data in the Euler deconvolution process, the preset fracture location prediction model is used to predict the fracture location and obtain the fracture location data; at the same time, the fracture type data is obtained based on fuzzy clustering;
[0055] Based on the area to be measured, the airborne gravity and magnetic data to be measured are obtained; the initial estimated fracture parameters of the area to be measured are confirmed, the airborne gravity and magnetic data are estimated based on the mapping model, the loss function between the airborne gravity and magnetic data to be measured and the estimated airborne gravity and magnetic data is calculated, and it is determined whether the loss function value exceeds the preset accuracy threshold. If it does not exceed the preset accuracy threshold, the fracture characteristic parameters are optimized through particle swarm intelligence; if so, the iteration is terminated and the optimal fracture characteristic parameters are output;
[0056] Outputs the fracture location, fracture type, and fracture parameters output in the final iteration step.
[0057] The source body with the center point at (x0, y0, z0) satisfies the following Euler equation:
[0058]
[0059] The parameter N in the homogeneous equation is defined as the tectonic index, which needs to be determined based on the source shape or known information about the anomaly's properties. It has a certain correspondence with the regular shape of the source (Table 1). For example, for a gravity source, the tectonic index of a narrow two-dimensional rock wall is 0, while the tectonic index of a uniform mass sphere is 2. For a given source geometry, a tectonic index is an exponential factor corresponding to the rate field drop-off and distance.
[0060] Table 1 Correspondence between field source and structural index
[0061] Gravity anomaly source Tectonic Index Magnetic anomaly source Tectonic Index rock wall or fault structure 0 Point Pole 1 Uniform mass horizontal cylinder 1 <![CDATA[Inclined thin plate (Z a )]]> 1 Limited inclination steps 1 Uniformly magnetized horizontal cylinder (bipolar line) 2 Uniform mass sphere 2 Uniformly magnetized sphere (dipole) 3
[0062] When performing Euler deconvolution inversion calculations, it is generally necessary to consider the influence of the regional field or background field B. The potential field is regarded as the sum of the point source field and the regional field, and the Euler equation is:
[0063]
[0064] By solving the Euler equation, the source position and depth can be inverted.
[0065] In a specific embodiment, the historical airborne gravity and magnetic data are obtained by gravity and magnetic measurements at two different measurement altitudes, namely, the first historical airborne gravity and magnetic data Q1 and the second historical airborne gravity and magnetic data Q2.
[0066] In a specific embodiment, the fracture characteristic parameters include fracture length, fracture width, fracture inclination and fracture type.
[0067] In a specific embodiment, the mapping model specifically includes:
[0068]
[0069] Where l is the fracture length, w is the fracture width, d is the fracture depth, θ is the fracture inclination, F c is the fault type, t1 and t2 are the measured heights, G1 and G2 are the historical gravity and magnetic data, and (x, y) is the field source position.
[0070] In a specific embodiment, predicting the fracture position specifically includes:
[0071] Each particle i represents a set of SVM parameter combinations p i =(γ i ,C i ), where γ i is the kernel parameter, C i is the penalty factor;
[0072]
[0073] Where k is the number of iterations, is the velocity vector of the i-th particle at the k+1-th iteration, corresponding to the velocity components of the kernel parameter γ and the penalty factor C, ω is the inertia weight, c1 and c2 are learning factors, r1 and r2 are random numbers between [0,1], pBest i,γ and pBest i,C is the kernel parameter and penalty factor corresponding to the optimal fitness value achieved by the i-th particle in its own history, gBest γ and gBest C is the global optimal kernel parameter and penalty factor corresponding to the optimal fitness value of all particles in the historical iteration, and is the kernel parameter and penalty factor of the i-th particle at the k-th iteration, m is the number of verification samples, is the SVM parameter combination (γ i ,C i ) The result of predicting the j-th validation sample, y j is the actual fracture position label of the jth verification sample;
[0074]
[0075] Among them, r is the weight vector, b is the bias term, ξ i It is a slack variable used to deal with the situation of linear inseparability, and C is a penalty factor that controls the degree of penalty for misclassified samples.
[0076] Solving the above optimization problem by the Lagrange multiplier method yields the dual problem:
[0077] SVM model:
[0078]
[0079] 0≤α i ≤C,i=1,…,n;
[0080] K(x i ,x j )=exp(-γ‖x i -x j ‖ 2 );
[0081]
[0082] Among them, x i 、x j , x are the input sample data, which are gravity gradient, magnetic gradient and depth data respectively; γ is the kernel parameter (radial basis function RBF), which controls the width of the kernel function and affects the fitting and generalization ability of the model; α i With αj are all Lagrange multipliers, obtained by solving the dual problem; y i is the category label of the sample, which can encode the fracture position in fracture position prediction; b is the bias term;
[0083] Initialize the particle swarm and randomly generate a set of initial positions of particles and speed
[0084] For each particle i, use its SVM parameter combination to train the SVM model and solve the Lagrange multiplier α i , α j And bias term b, calculate the fitness value Fitness i ;
[0085] Update pBest of each particle i,γ 、pBest i,C and gBest γ 、gBest C ;
[0086] Update the particle's speed and position according to the particle speed and position update formula;
[0087] Repeat the iterative update until the maximum number of iterations is reached to obtain the final fracture location prediction model and perform fracture location prediction.
[0088] In a specific embodiment, the gravity gradient and the magnetic gradient are ▽G and ▽M;
[0089]
[0090] in, is the gradient of the gravity field G in the x, y, and z directions; is the gradient of the magnetic field M in the x, y, and z directions; the aerogravimetric data includes the data of the gravity field and the magnetic field.
[0091] In a specific embodiment, obtaining fracture type data based on fuzzy clustering specifically includes:
[0092] Determine the number of clusters c of the fracture type and select the fuzzy factor m, m∈(1,+∞), which is usually 2; set the iteration stop threshold ε, such as ε=10 -6 , used to determine whether the iteration converges; construct n×c membership matrix U=[u ij ],u ij For sample q i The degree of membership to the jth class, Sample q i Including gravity gradient, magnetic field gradient, and depth data;
[0093] Calculate the cluster center v of each class according to the membership matrix U j , the calculation formula is:
[0094]
[0095] Among them, v j is the cluster center of the jth class, which is obtained by weighted averaging the sample data, with the weight determined by the mth power of the membership degree;
[0096] Using the current cluster center v j Update the membership matrix U, and the update formula is:
[0097]
[0098] Among them, ||q i -v j || represents sample x i With cluster center v j distance;
[0099] Define the objective function J of fuzzy clustering, the formula is:
[0100]
[0101] Calculate the objective function value J of this iteration (t) and the objective function value J of the previous iteration (t-1) , if J is satisfied (t) -J (t -1) ≤ε, the iteration is considered to have converged and the iteration is stopped; otherwise, the iteration is continued, and the cluster center and membership matrix are continuously updated until the termination condition is met; when the iteration is completed, the fracture type of each sample is determined according to the final membership matrix U, and for each sample q i , the fault type is the category with the largest membership, that is In this way, all samples are divided into corresponding fault types, and the fault type data acquisition based on fuzzy clustering is completed.
[0102] In a specific embodiment, the loss function specifically includes:
[0103]
[0104] Among them, Loss is the loss function value, n represents the number of data samples, and are the measured gravity data and measured magnetic data of the i-th sample, and are the historical gravity data and historical magnetic data of the i-th sample, respectively, which are the data calculated based on the established mapping model. The historical gravity and magnetic data include historical gravity data and historical magnetic data.
[0105] The various embodiments in this specification are described in a progressive manner, with each embodiment focusing on the differences from other embodiments. Reference can be made to the common and similar parts between the various embodiments. For the devices disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple, and the relevant parts can be referred to the method description.
[0106] The above description of the disclosed embodiments is intended to enable one skilled in the art to implement or use the present invention. Various modifications to these embodiments will be readily apparent to one skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention is not limited to the embodiments shown herein but is intended to conform to the widest scope consistent with the principles and novel features disclosed herein.
Claims
1. A fracture identification method based on airborne gravity and magnetic data, characterized in that: include: Acquire historical airborne gravity and magnetic data, process the data through Euler deconvolution, and obtain historical field source location and depth data; establish a mapping model between the historical airborne gravity and magnetic data and fault characteristic parameters, field source location, and depth data based on a neural network algorithm; Based on the historical source location and depth data in the Euler deconvolution process, the preset fracture location prediction model is used to predict the fracture location and obtain the fracture location data; at the same time, the fracture type data is obtained based on fuzzy clustering; Based on the area to be measured, obtaining the airborne gravity and magnetic data to be measured; confirming the initial estimated fracture characteristic parameters of the area to be measured, calculating the estimated airborne gravity and magnetic data based on the mapping model, calculating the loss function between the airborne gravity and magnetic data to be measured and the estimated airborne gravity and magnetic data, and determining whether the loss function value exceeds a preset accuracy threshold; if it does not exceed the preset accuracy threshold, optimizing the fracture characteristic parameters through particle swarm intelligence; if so, ending the iteration and outputting the optimal fracture characteristic parameters; Output the fracture position and fracture characteristic parameters output in the final iteration step; The fracture characteristic parameters include fracture length, fracture width, fracture dip and fracture type; The mapping model specifically includes: Where l is the fracture length, w is the fracture width, d is the fracture depth, θ is the fracture inclination, F c is the fault type, t1 and t2 are the measured heights, G1 and G2 are the historical airborne gravity and magnetic data, and (x, y) is the source location.
2. A method for identifying fractures based on aerial gravity and magnetic data according to claim 1, characterized in that: The historical airborne gravity and magnetic data are obtained by gravity and magnetic measurements at two different measurement altitudes, namely the first historical airborne gravity and magnetic data Q1 and the second historical airborne gravity and magnetic data Q2.
3. The method for identifying fractures based on aerial gravity and magnetic data according to claim 1, characterized in that: The fracture position prediction specifically includes: Each particle i represents a set of SVM parameter combinations p i =(γ i ,C i ), where γ i is the kernel parameter, C i is the penalty factor; Where k is the number of iterations, is the velocity vector of the i-th particle at the k+1-th iteration, corresponding to the velocity components of the kernel parameter γ and the penalty factor C, ω is the inertia weight, c1 and c2 are learning factors, r1 and r2 are random numbers between [0,1], pBest i,γ and pBest i,C is the kernel parameter and penalty factor corresponding to the optimal fitness value achieved by the i-th particle in its own history, gBest γ and gBest C is the global optimal kernel parameter and penalty factor corresponding to the optimal fitness value of all particles in the historical iteration, and is the kernel parameter and penalty factor of the i-th particle at the k-th iteration, m is the number of verification samples, is the SVM parameter combination (γ i ,C i ) The result of predicting the j-th validation sample, y j is the actual fracture position label of the jth verification sample; SVM model: 0≤α i ≤C,i=1,…,n; K(x i ,x j )=exp(-γ‖x i -x j ‖ 2 ); Among them, x i 、x j , x are the input sample data, which are gravity gradient, magnetic gradient and depth data respectively; γ is the kernel parameter, α i With α j are all Lagrange multipliers, y i is the category label of the sample, b is the bias term, and n is the number of sample data; Initialize the particle swarm and randomly generate a set of initial positions of particles and speed For each particle i, use its SVM parameter combination to train the SVM model and solve the Lagrange multiplier α i , α j And bias term b, calculate the fitness value Fitness i ; Update pBest of each particle i,γ 、pBest i,C and gBest γ 、gBest C ; Update the particle's speed and position according to the particle speed and position update formula; Repeat the iterative update until the maximum number of iterations is reached to obtain the final fracture location prediction model and perform fracture location prediction.
4. The method for identifying fractures based on aerial gravity and magnetic data according to claim 1, characterized in that: The acquisition of fracture type data based on fuzzy clustering specifically includes: Determine the number of clusters c of the fracture type and select the fuzzy factor m, m∈(1,+∞), set the iteration stop threshold ε; construct the n×c membership matrix U=[u ij ],u ij For sample q i The degree of membership to the jth class, Calculate the cluster center v of each class according to the membership matrix U j , the calculation formula is: Among them, v j is the cluster center of the jth class, which is obtained by weighted averaging the sample data, with the weight determined by the mth power of the membership degree; Using the current cluster center v j Update the membership matrix U, and the update formula is: Among them, ||q i -v j || represents sample x i With cluster center v j distance; Define the objective function J of fuzzy clustering, the formula is: Calculate the objective function value J of this iteration (t) and the objective function value J of the previous iteration (t-1) , if J is satisfied (t) -J (t-1) ≤ε, the iteration is considered to have converged and the iteration is stopped; otherwise, the iteration is continued, and the cluster center and membership matrix are continuously updated until the termination condition is met; when the iteration is completed, the fracture type of each sample is determined according to the final membership matrix U, and for each sample q i , the fracture type is the category with the largest membership.
5. The method for identifying fractures based on aerial gravity and magnetic data according to claim 1, characterized in that: The loss function specifically includes: Among them, Loss is the loss function value, n represents the number of data samples, and are respectively the measured airborne gravity data and the measured airborne magnetic data of the i-th sample, and are the historical airborne gravity data and historical airborne magnetic data of the i-th sample, which are the data calculated based on the established mapping model. The historical airborne gravity and magnetic data include historical airborne gravity data and historical airborne magnetic data.
Citation Information
Patent Citations
SVM classifier parameter optimization method based on improved particle swarm algorithm
CN108875788A
Intelligent identification method and identification system based on support vector machine and civil aviation engine
CN111582510A