A method for calibrating inertial sensors in a measurement while drilling system

Abstract

The invention discloses an inertial sensor calibration method in a while drilling measurement system, comprising: introducing a Tent chaotic map to generate an initial population, an adaptive factor and a logarithmic function, statically improving a traditional northern goshawk optimization algorithm, thereby statically calibrating an error parameter of the inertial sensor; loading the static calibration error parameter value as initial prior knowledge, simulating an explosion phenomenon, and dynamically improving the traditional northern goshawk optimization algorithm, thereby dynamically calibrating and updating the error parameter of the inertial sensor when drilling is stopped at a working measuring point of the while drilling measurement system; the invention respectively designs a static calibration step of the inertial sensor and a dynamic calibration step of the inertial sensor, avoids the cumbersome steps of the traditional calibration scheme of the inertial sensor, and proposes static improvement and dynamic improvement for the traditional northern goshawk optimization algorithm, thereby improving the speed and accuracy of the error calibration.

Description

Technical Field

[0001] The invention belongs to the technical field of inertial navigation inertial sensor error calibration, and in particular relates to an inertial sensor calibration method in a while drilling measurement system. Background Art

[0002] Measurement While Driling (MWD) is an important technical branch of underground directional drilling technology in coal mines. It is mainly used for monitoring while drilling during the construction of medium-deep hole directional drilling. As the name suggests, it measures various parameters of the borehole during the drilling process, including (pitch angle, azimuth, roll angle, well depth) trajectory parameters, (resistivity, natural gamma, porosity) geological information parameters and (torque, drilling pressure, temperature at the bottom of the well) engineering parameters. The drill bit is positioned and the trajectory is drawn based on the axis trace information obtained by the MWD short section, so that the driller can understand the drilling construction situation at any time, and adjust the drill tool combination and drilling process parameters in time according to the measurement results, so that the borehole can extend in the predetermined direction as much as possible. This technology was first developed and applied to the oil and gas extraction industry abroad, and then gradually applied to the coal industry. It is currently mostly used in geological exploration, trenchless and other drilling technology fields.

[0003] The complex working environment in coal mines makes it difficult for other navigation technologies to play a role. The inertial navigation system has become the main means of navigation in coal mines because it does not rely on external information and does not radiate energy to the outside. However, due to the process, materials, working environment and other reasons, inertial sensors usually need to be statically calibrated before use. The impact of the dynamic environment on the error terms is analyzed to establish a dynamic compensation model for dynamic compensation of sensor measurements. The current inertial sensor calibration technology has the following problems:

[0004] (1) Most inertial sensor calibrations are performed in a laboratory environment, and the calibrated error parameters are not updated during system operation. However, the measurement while drilling system is affected by factors such as vibration and temperature changes during operation, which can cause changes in the zero bias value and proportional factor error of the inertial sensor. Therefore, it is unreliable to rely solely on the initial static laboratory calibration.

[0005] (2) In order to solve the first problem mentioned above, many scholars have studied the specific impact of factors such as vibration and temperature changes on sensor error parameters. They have established a functional relationship between temperature and error parameters through parameter identification, polynomial fitting, and neural network prediction, and used the established functional relationship to perform real-time dynamic compensation of error parameters. However, since there are many factors that cause error parameters to change and the compensation steps are usually complicated, the dynamic compensation solution is not only costly, but also difficult to compensate for the influence of all factors.

[0006] It can be seen that the current solution of static calibration + dynamic compensation of inertial sensors still has a lot of room for improvement. Summary of the invention

[0007] In order to overcome the defects of the prior art, the present invention proposes a dynamic calibration method for inertial sensors in a measurement while drilling system. The method utilizes the characteristic of the measurement while drilling system that drilling is stopped at a measuring point to calibrate and update the error terms, thereby avoiding cumbersome dynamic compensation steps for the error terms and improving the working efficiency of the inertial sensor.

[0008] The purpose of the present invention can be achieved by the following technical solutions:

[0009] A method for calibrating an inertial sensor in a measurement while drilling system comprises statically calibrating an error parameter of the inertial sensor, and dynamically calibrating and updating the error parameter of the inertial sensor when drilling is stopped at a working measuring point of the measurement while drilling system;

[0010] The specific steps of the static calibration of the error parameters include:

[0011] S111: Establish an inertial sensor error model, including an accelerometer error model and a gyroscope error model, introduce gravity acceleration and the earth's rotation speed as calibration benchmarks, read accelerometer measurement data, gyroscope measurement data and calibration benchmark information, and establish a static calibration objective function;

[0012] S112: Based on the objective function, analyzing the traditional Northern Goshawk optimization algorithm, and statically improving the traditional Northern Goshawk optimization algorithm, wherein the static improvement includes:

[0013] The method of generating initial population by using Tent chaotic mapping is introduced to enhance the diversity of initial individuals in the solution space of static improved optimization algorithm.

[0014] Adaptive factors and logarithmic functions are introduced to balance the relationship between global search and local search of the static improvement optimization algorithm, and to enhance the speed and accuracy of static improvement optimization algorithm in static calibration of error parameters.

[0015] S113: setting initialization parameters of the static improved optimization algorithm, generating an initial population of the static improved optimization algorithm, and calculating an objective function value of the initial population of the static improved optimization algorithm;

[0016] S114: Iteratively optimize the objective function value using the statically improved Northern Goshawk optimization algorithm until the objective function value satisfies the constraint condition, and output the static calibration error parameter value;

[0017] The specific steps of the error parameter dynamic calibration and error parameter update include:

[0018] S121: Establish an inertial sensor error model, including an accelerometer error model and a gyroscope error model, introduce gravity acceleration and the earth's rotation speed as calibration benchmarks, read accelerometer measurement data, gyroscope measurement data and calibration benchmark information, and establish a dynamic calibration objective function;

[0019] S122: Based on the objective function, analyzing the traditional northern goshawk optimization algorithm, and dynamically improving the traditional northern goshawk optimization algorithm, wherein the dynamic improvement includes:

[0020] Loading the static calibration error parameter value as initial prior knowledge to obtain the next dynamic calibration prediction error parameter value;

[0021] Centering on the next dynamic calibration prediction error parameter value, simulating the explosion phenomenon, generating the initial population of the dynamic improvement optimization algorithm in the solution space, and calculating the objective function value of the dynamic improvement optimization algorithm;

[0022] Adaptive factors and logarithmic functions are introduced to balance the relationship between global search and local search of the dynamic improvement optimization algorithm, and to enhance the speed and accuracy of dynamic calibration of error parameters of the dynamic improvement optimization algorithm.

[0023] S123: setting initialization parameters of the dynamic improvement optimization algorithm, generating an initial population of the dynamic improvement optimization algorithm, and calculating an objective function value of the initial population of the dynamic improvement optimization algorithm;

[0024] S124: Iteratively optimize the objective function value using the dynamically improved Northern Goshawk optimization algorithm until the objective function value satisfies the constraint condition, and output the dynamic calibration error parameter value;

[0025] S125: using the dynamic calibration error parameter value to update the static calibration error parameter value;

[0026] S126: Repeat S123 to S124, and use the dynamic calibration error parameter value to update the last calibration error parameter value, until a stop instruction is issued.

[0027] Furthermore, the accelerometer error model is shown in formula (S-1):

[0028]

[0029] In formula (S-1), represents the accelerometer measurement value, represents the true value of the accelerometer, b a is the accelerometer bias error, A abThe transformation matrix from the carrier coordinate system b to the accelerometer coordinate system a is defined as follows:

[0030]

[0031] In formula (S-2), s and t are the intermediate temporary coordinate systems constructed, S at is the scale factor error matrix, as shown in formula (S-3); is a non-orthogonal error matrix, as shown in formula (S-4); is the misalignment error matrix, as shown in formula (S-5):

[0032]

[0033] In formula (S-3), They represent the accelerometer scale factor error parameters respectively;

[0034]

[0035] In formula (S-4), They represent the non-orthogonal error parameters of the accelerometer respectively;

[0036]

[0037] In formula (S-5), They represent the misalignment error parameters between the carrier coordinate system and the accelerometer coordinate system respectively.

[0038] Furthermore, the gyroscope error model is shown in formula (S-6):

[0039]

[0040] In formula (S-6), represents the gyroscope measurement value, represents the true value of the gyroscope, b g represents the gyro bias error, A gb represents the transformation matrix from the carrier coordinate system b to the gyroscope coordinate system g, which is defined as formula (S-7):

[0041]

[0042] In formula (S-7), S gt is the scale factor error matrix, as shown in formula (S-8), is a non-orthogonal error matrix, as shown in formula (S-9), is the misalignment error matrix, as shown in formula (S-10);

[0043]

[0044] In formula (S-8), They represent the gyroscope scale factor error parameters respectively;

[0045]

[0046] In formula (S-9), They represent the gyroscope non-orthogonal error parameters respectively;

[0047]

[0048] In formula (S-10), represent the misalignment error parameters between the carrier coordinate system and the gyroscope coordinate system respectively.

[0049] Furthermore, the static calibration objective function includes an accelerometer objective function and a gyroscope objective function. The accelerometer objective function is as shown in formula (S-11), and the gyroscope objective function is as shown in formula (S-12):

[0050]

[0051] In formulas (S-11) to (S-12), G represents the local gravity acceleration vector, ω represents the modulus of the earth's rotation velocity, l represents the number of collected data points, and n represents the total number of collected data points;

[0052] When the measurement while drilling system is working, it includes the first measurement while drilling feature and the second measurement while drilling feature:

[0053] The first characteristic of MWD is that the non-magnetic drill collar will stop drilling after drilling a certain distance. At this time, the MWD system will perform geological measurement, drill bit attitude estimation and drilling parameter measurement.

[0054] The second measurement while drilling feature includes that no matter what attitude and position the inertial measurement element is in, the accelerometer measurement output modulus is always equal to the local gravity acceleration modulus, and the gyroscope measurement output modulus is always equal to the earth's rotation rate.

[0055] based on and The calculations are as shown in formula (S-13) and formula (S-14);

[0056]

[0057] According to the second measurement while drilling characteristics, the non-alignment error parameters between the carrier coordinate system and the accelerometer coordinate system, the accelerometer zero bias error parameters, the accelerometer scale factor error parameters and the accelerometer non-orthogonal error parameters are used to construct a solution vector of the accelerometer static calibration objective function, as shown in formula (S-15):

[0058]

[0059] According to the misalignment error parameter between the carrier coordinate system and the gyroscope coordinate system, the gyroscope zero bias error parameter, the gyroscope scale factor error parameter and the gyroscope non-orthogonal error parameter, a solution vector of the gyroscope dynamic calibration objective function is constructed, as shown in formula (S-16):

[0060]

[0061] Based on the above analysis, the inertial sensor error calibration problem can be transformed into the problem of finding the minimum value of J(a) and J(g). The solution vector is adjusted through the intelligent optimization algorithm to continuously approach J(a) to 0, thereby obtaining the sensor error parameters to compensate the sensor measurement output.

[0062] The traditional northern goshawk optimization algorithm is as follows:

[0063]

[0064] In formula (S-17), X represents the population of northern goshawk, X i represents the i-th solution in the population, x i,j represents the jth variable value in the i-th solution individual in the population, N represents the number of individuals in the population, and m represents the number of parameters to be optimized in the problem;

[0065]

[0066] In formula (S-18), F(X) is the solution vector of the objective function, F i represents the objective function value of the ith individual. The hunting process simulated by the northern goshawk optimization algorithm can be divided into two stages. The first stage is the global prey search stage, which simulates the northern goshawk's prey recognition and attack behavior; the second stage is the local search stage, which simulates the northern goshawk's attack and rapid approach to the prey, and the prey tries to escape, and the northern goshawk pursues it.

[0067] The detailed expression of the global prey search stage is shown in equations (S-19) to (S-21):

[0068] P i =X k ,i=1,2,…,N;k=1,2,…,i-1,i,i+1,…,N(S-19)

[0069] In formula (S-19), P i is the prey position of the i-th northern goshawk, k is a random natural number between [1, N], indicating that the prey position is randomly selected in the initial population;

[0070]

[0071] In formula (S-20), It represents the jth dimension of the i-th new solution, is the corresponding objective function value, x i,j represents the jth dimension of the i-th solution of the initial population in the last iteration, p i,j Represents the jth dimension of the location of the prey selected by the northern goshawk this time; r is a random number between [0,1], the value of I is randomly selected between 1 or 2, and r and I are used to randomly generate NGO behavior patterns.

[0072]

[0073] In formula (S-21), represents the new state of the ith solution, It represents the objective function value of the NGO global search for prey phase corresponding to the new state. This formula describes the greedy strategy of the northern goshawk optimization algorithm, that is, retaining the optimal individual in the current stage.

[0074] The detailed description of the local search stage is shown in Equation (S-22) to Equation (S-24):

[0075]

[0076] In formula (S-22), represents the value of its jth dimension, and R is the assumed range radius of this hunting activity;

[0077]

[0078] In formula (S-23), t represents the iteration counter and T represents the maximum number of iterations.

[0079]

[0080] In formula (S-24), represents the new state of the ith solution, It represents the objective function value of the NGO local search phase corresponding to the new state. This formula also describes the greedy strategy of the Northern Goshawk optimization algorithm.

[0081] Furthermore, the method of introducing the Tent chaotic map to generate the initial population is as shown in formula (S-25):

[0082]

[0083] In formula (S-25), Z k represents the kth chaotic value, whose initial value is Z 0 ∈[0,1], β∈[0,1], the value of β affects the distribution of chaotic values;

[0084] The adaptive factor and logarithmic function are introduced to balance the relationship between the global search and the local search of the static improved optimization algorithm, as shown in formula (S-26):

[0085]

[0086] In formula (S-26), w represents the adaptive factor. The larger the value, the slower the parameter R decays in the early stage of the algorithm. In engineering, a suitable adaptive factor w should be selected to speed up the global search of the optimization algorithm and enhance the fine search capability of the optimization algorithm in the local search stage. The introduction of the logarithmic function can further enhance the calibration performance of the algorithm.

[0087] Formula (S-26) is used instead of the Northern Goshawk optimization algorithm formula (S-23), so that the parameter R is no longer adjusted in the form of a linear function. Instead, a larger value is taken in the early stage of the algorithm to ensure the global optimization ability of the algorithm. A smaller value is taken in the later stage of the algorithm to enable the algorithm to perform a more detailed search in a better area and avoid the algorithm falling into a local optimum.

[0088] Furthermore, since the non-orthogonal error and non-alignment error of the inertial sensor have been obtained after static calibration and can be used as a known error matrix, and do not change with time and environment, the dynamic calibration of the error parameters includes calibrating the accelerometer zero bias error and the accelerometer scale factor error; it also includes calibrating the gyroscope zero bias error and the gyroscope scale factor error.

[0089] Furthermore, the dynamic calibration objective function is consistent with the static calibration objective function, as shown in equations (S-11) to (S-12);

[0090] The solution vector of the accelerometer dynamic calibration objective function is constructed using the accelerometer zero bias error parameter and the accelerometer scale factor error parameter, as shown in formula (S-27):

[0091]

[0092] Similarly, the solution vector of the gyroscope dynamic calibration objective function is constructed, as shown in formula (S-28):

[0093]

[0094] The static calibration error parameter value is loaded as the initial prior knowledge to obtain the next dynamic calibration prediction error parameter value, as shown in formula (S-29):

[0095]

[0096] In formula (S-28), represents the difference vector of the error parameter values ​​of the previous two calibrations, k represents the current kth dynamic calibration of the measuring point, X k-1,Best , X k-2,Best They represent the error parameter values ​​of the k-1th and k-2th calibrations respectively;

[0097] according to The k-th calibration error parameter value is predicted and generated, as shown in formula (S-30):

[0098]

[0099] Furthermore, the explosion phenomenon is described by explosively generating the initial population of the k-th calibration algorithm, as shown in formula (S-31):

[0100] X=normrnd(aver,sigma,m,n) (S-31)

[0101] In formula (S-31), the normrnd function represents the generation of a set of random numbers that conform to the normal distribution, which contains 4 parameters, among which: aver represents the mean of the normal distribution, that is, the value of each dimension of the prediction error parameter value, sigma represents the standard deviation of the normal distribution, which determines the range of the explosion. The larger the value, the larger the explosion range, and the greater the probability of particles falling far away, otherwise it is smaller, and m,n represents the generation of an m×n matrix.

[0102] Furthermore, the constraint conditions include a convergence accuracy constraint and a maximum number of iterations constraint.

[0103] Compared with the prior art, the present invention has the following technical effects:

[0104] The inertial sensor calibration method in the measurement while drilling system proposed in the present invention respectively designs a static calibration step and a dynamic calibration step of the inertial sensor, thereby avoiding the cumbersome steps of the traditional calibration scheme of the inertial sensor, and proposes static improvement and dynamic improvement for the traditional Northern Goshawk optimization algorithm, thereby improving the speed and accuracy of error calibration. BRIEF DESCRIPTION OF THE DRAWINGS

[0105] Figure 1 It is a calibration flow chart based on the static improved northern goshawk optimization algorithm of the present invention;

[0106] Figure 2It is a calibration flow chart based on the dynamically improved northern goshawk optimization algorithm of the present invention;

[0107] Figure 3 This is an effect diagram of the present invention introducing an adaptive factor and a logarithmic function;

[0108] Figure 4 It is a comparison diagram of the iteration curves of the present invention. DETAILED DESCRIPTION

[0109] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments in the invention, all other embodiments obtained by ordinary technicians in this field without creative work, any modifications, equivalent substitutions, improvements, etc., should be included in the protection scope of the present invention.

[0110] A method for calibrating an inertial sensor in a measurement while drilling system comprises statically calibrating an error parameter of the inertial sensor, and dynamically calibrating and updating the error parameter of the inertial sensor when drilling is stopped at a working measuring point of the measurement while drilling system;

[0111] like Figure 1 As shown, the specific steps of the static calibration of the error parameters include:

[0112] S111: Establish an inertial sensor error model, including an accelerometer error model and a gyroscope error model, introduce gravity acceleration and the earth's rotation speed as calibration benchmarks, read accelerometer measurement data, gyroscope measurement data and calibration benchmark information, and establish a static calibration objective function;

[0113] S112: Based on the objective function, a traditional Northern Goshawk Optimization algorithm (NGO) is analyzed, and a static improvement is performed on the traditional Northern Goshawk Optimization algorithm to obtain an Adaptive Northern Goshawk Optimization algorithm (ANGO), wherein the static improvement includes:

[0114] The method of generating initial population by using Tent chaotic mapping is introduced to enhance the diversity of initial individuals in the solution space of static improved optimization algorithm.

[0115] Adaptive factors and logarithmic functions are introduced to balance the relationship between global search and local search of the static improvement optimization algorithm, and to enhance the speed and accuracy of static improvement optimization algorithm in static calibration of error parameters.

[0116] S113: setting initialization parameters of the static improved optimization algorithm, generating an initial population of the static improved optimization algorithm, and calculating an objective function value of the initial population of the static improved optimization algorithm;

[0117] S114: Iteratively optimize the objective function value using the statically improved Northern Goshawk optimization algorithm until the objective function value satisfies the constraint condition, and output the static calibration error parameter value;

[0118] like Figure 2 As shown, the specific steps of the error parameter dynamic calibration and error parameter update include:

[0119] S121: Establish an inertial sensor error model, including an accelerometer error model and a gyroscope error model, introduce gravity acceleration and the earth's rotation speed as calibration benchmarks, read accelerometer measurement data, gyroscope measurement data and calibration benchmark information, and establish a dynamic calibration objective function;

[0120] S122: Based on the objective function, the traditional Northern Goshawk optimization algorithm is analyzed, and the traditional Northern Goshawk optimization algorithm is dynamically improved to obtain a Dynamic Adaptive Northern Goshawk Optimization algorithm (DANGO), wherein the dynamic improvement includes:

[0121] Loading the static calibration error parameter value as initial prior knowledge to obtain the next dynamic calibration prediction error parameter value;

[0122] Centering on the next dynamic calibration prediction error parameter value, simulating the explosion phenomenon, generating the initial population of the dynamic improvement optimization algorithm in the solution space, and calculating the objective function value of the dynamic improvement optimization algorithm;

[0123] Adaptive factors and logarithmic functions are introduced to balance the relationship between global search and local search of the dynamic improvement optimization algorithm, and to enhance the speed and accuracy of dynamic calibration of error parameters of the dynamic improvement optimization algorithm.

[0124] S123: setting initialization parameters of the dynamic improvement optimization algorithm, generating an initial population of the dynamic improvement optimization algorithm, and calculating an objective function value of the initial population of the dynamic improvement optimization algorithm;

[0125] S124: Iteratively optimize the objective function value using the dynamically improved Northern Goshawk optimization algorithm until the objective function value satisfies the constraint condition, and output the dynamic calibration error parameter value;

[0126] S125: using the dynamic calibration error parameter value to update the static calibration error parameter value;

[0127] S126: Repeat S123 to S124, and use the dynamic calibration error parameter value to update the last calibration error parameter value, until a stop instruction is issued.

[0128] The accelerometer error model is shown in formula (S-1):

[0129]

[0130] In formula (S-1), represents the accelerometer measurement value, represents the true value of the accelerometer, b a is the accelerometer bias error, A ab The transformation matrix from the carrier coordinate system b to the accelerometer coordinate system a is defined as follows:

[0131]

[0132] In formula (S-2), s and t are the intermediate temporary coordinate systems constructed, S at is the scale factor error matrix, as shown in formula (S-3); is a non-orthogonal error matrix, as shown in formula (S-4); is the misalignment error matrix, as shown in formula (S-5):

[0133]

[0134] In formula (S-3), They represent the accelerometer scale factor error parameters respectively;

[0135]

[0136] In formula (S-4), They represent the non-orthogonal error parameters of the accelerometer respectively;

[0137]

[0138] In formula (S-5), They represent the misalignment error parameters between the carrier coordinate system and the accelerometer coordinate system respectively.

[0139] The gyroscope error model is shown in formula (S-6):

[0140]

[0141] In formula (S-6), represents the gyroscope measurement value, represents the true value of the gyroscope, b g represents the gyro bias error, A gb represents the transformation matrix from the carrier coordinate system b to the gyroscope coordinate system g, which is defined as formula (S-7):

[0142]

[0143] In formula (S-7), S gt is the scale factor error matrix, as shown in formula (S-8), is a non-orthogonal error matrix, as shown in formula (S-9), is the misalignment error matrix, as shown in formula (S-10);

[0144]

[0145] In formula (S-8), They represent the gyroscope scale factor error parameters respectively;

[0146]

[0147] In formula (S-9), They represent the gyroscope non-orthogonal error parameters respectively;

[0148]

[0149] In formula (S-10), represent the misalignment error parameters between the carrier coordinate system and the gyroscope coordinate system respectively.

[0150] The static calibration objective function includes an accelerometer objective function and a gyroscope objective function. The accelerometer objective function is as shown in formula (S-11), and the gyroscope objective function is as shown in formula (S-12):

[0151]

[0152] In formulas (S-11) to (S-12), G represents the local gravity acceleration vector, ω represents the modulus of the earth's rotation velocity, l represents the number of collected data points, and n represents the total number of collected data points;

[0153] When the measurement while drilling system is working, it includes the first measurement while drilling feature and the second measurement while drilling feature:

[0154] The first characteristic of MWD is that the non-magnetic drill collar will stop drilling after drilling a certain distance. At this time, the MWD system will perform geological measurement, drill bit attitude estimation and drilling parameter measurement.

[0155] The second measurement while drilling feature includes that no matter what attitude and position the inertial measurement element is in, the accelerometer measurement output modulus is always equal to the local gravity acceleration modulus, and the gyroscope measurement output modulus is always equal to the earth's rotation rate.

[0156] based on and The calculations are as shown in formula (S-13) and formula (S-14);

[0157]

[0158] According to the second measurement while drilling characteristics, the non-alignment error parameters between the carrier coordinate system and the accelerometer coordinate system, the accelerometer zero bias error parameters, the accelerometer scale factor error parameters and the accelerometer non-orthogonal error parameters are used to construct a solution vector of the accelerometer static calibration objective function, as shown in formula (S-15):

[0159]

[0160] According to the misalignment error parameter between the carrier coordinate system and the gyroscope coordinate system, the gyroscope zero bias error parameter, the gyroscope scale factor error parameter and the gyroscope non-orthogonal error parameter, a solution vector of the gyroscope dynamic calibration objective function is constructed, as shown in formula (S-16):

[0161]

[0162] Based on the above analysis, the inertial sensor error calibration problem can be transformed into the problem of finding the minimum value of J(a) and J(g). The solution vector is adjusted through the intelligent optimization algorithm to continuously approach J(a) to 0, thereby obtaining the sensor error parameters to compensate the sensor measurement output.

[0163] The traditional northern goshawk optimization algorithm is as follows:

[0164]

[0165] In formula (S-17), X represents the population of northern goshawk, X i represents the i-th solution in the population, x i,j represents the jth variable value in the i-th solution individual in the population, N represents the number of individuals in the population, and m represents the number of parameters to be optimized in the problem;

[0166]

[0167] In formula (S-18), F(X) is the solution vector of the objective function, F irepresents the objective function value of the ith individual. The hunting process simulated by the northern goshawk optimization algorithm can be divided into two stages. The first stage is the global prey search stage, which simulates the northern goshawk's prey recognition and attack behavior; the second stage is the local search stage, which simulates the northern goshawk's attack and rapid approach to the prey, and the prey tries to escape, and the northern goshawk pursues it.

[0168] The detailed expression of the global prey search stage is shown in equations (S-19) to (S-21):

[0169] P i =X k ,i=1,2,…,N;k=1,2,…,i-1,i,i+1,…,N(S-19)

[0170] In formula (S-19), P i is the prey position of the i-th northern goshawk, k is a random natural number between [1, N], indicating that the prey position is randomly selected in the initial population;

[0171]

[0172] In formula (S-20), It represents the jth dimension of the i-th new solution, is the corresponding objective function value, x i,j represents the jth dimension of the i-th solution of the initial population in the last iteration, p i,j Represents the jth dimension of the location of the prey selected by the northern goshawk this time; r is a random number between [0,1], the value of I is randomly selected between 1 or 2, and r and I are used to randomly generate NGO behavior patterns.

[0173]

[0174] In formula (S-21), represents the new state of the ith solution, It represents the objective function value of the NGO global search for prey phase corresponding to the new state. This formula describes the greedy strategy of the northern goshawk optimization algorithm, that is, retaining the optimal individual in the current stage.

[0175] The detailed description of the local search stage is shown in Equation (S-22) to Equation (S-24):

[0176]

[0177] In formula (S-22), represents the value of its Jth dimension, and R is the assumed range radius of this hunting activity;

[0178]

[0179] In formula (S-23), t represents the iteration counter and T represents the maximum number of iterations.

[0180]

[0181] In formula (S-24), represents the new state of the ith solution, It represents the objective function value of the NGO local search phase corresponding to the new state. This formula also describes the greedy strategy of the Northern Goshawk optimization algorithm.

[0182] The method of introducing Tent chaotic mapping to generate the initial population is as shown in formula (S-25):

[0183]

[0184] In formula (S-25), Z k represents the kth chaotic value, whose initial value is Z 0 ∈[0,1], β∈[0,1], the value of β affects the distribution of chaotic values; after several experiments, the present invention preferably takes β=0.4;

[0185] The adaptive factor and logarithmic function are introduced to balance the relationship between the global search and the local search of the static improved optimization algorithm, as shown in formula (S-26):

[0186]

[0187] In formula (S-26), w represents the adaptive factor. The larger the value, the slower the parameter R decays in the early stage of the algorithm. In engineering, a suitable adaptive factor w should be selected to speed up the global search of the optimization algorithm and enhance the fine search ability of the optimization algorithm in the local search stage. The introduction of the logarithmic function can further enhance the calibration performance of the algorithm. After several experiments, w = 1.2 is preferably selected.

[0188] The present invention adopts Matlab_R2021b software, takes the static calibration and dynamic calibration of the accelerometer as examples, and takes the accelerometer error model established by S111 and S121 as the object to carry out simulation experiments.

[0189] Firstly, 1000 points on the unit sphere are selected as gravity vectors at different positions using Matlab software to simulate the gravitational acceleration, that is, the true value; and the error parameter values ​​are set as shown in Table 1.

[0190] Table 1 Accelerometer error parameter settings

[0191]

[0192] The initial parameter settings of the algorithm are shown in Table 2.

[0193] Table 2 Initial parameter settings of static calibration algorithm

[0194]

[0195] Formula (S-26) is used instead of the Northern Goshawk optimization algorithm formula (S-23), so that the parameter R is no longer adjusted in the form of a linear function. Instead, a larger value is taken in the early stage of the algorithm to ensure the global optimization ability of the algorithm. A smaller value is taken in the later stage of the algorithm to enable the algorithm to perform a more detailed search in a better area and avoid the algorithm falling into a local optimum.

[0196] like Figure 3 As shown in the figure, by comparing the R value change curve when the adaptive factor and logarithmic function are introduced and the R value change curve when they are not introduced, it can be seen from the figure that the R value is larger and decays faster in the early stage of iteration, which is beneficial to expand the global search range of the optimization algorithm. In the later stage of iteration, the R value is smaller and decays slower, which is beneficial for the optimization algorithm to perform a more refined search during local search and prevent the algorithm from falling into the local optimum.

[0197] By setting the initial parameters of the algorithm, the algorithm stops iterating when it reaches the specified convergence accuracy or the number of iterations. As shown in Tables 3 and 4, the preset error parameters, the error parameter values ​​obtained by static calibration of the Northern Goshawk Optimization Algorithm (NGO), and the error parameter values ​​obtained by static calibration of the Adaptive Northern Goshawk Optimization Algorithm (ANGO) are compared.

[0198] Table 3 Comparison of calibration errors between the original algorithm and the proposed algorithm 1

[0199]

[0200] Table 4 Comparison of calibration errors between the original algorithm and the proposed algorithm 2

[0201]

[0202] It can be seen that the statically improved ANGO has higher calibration accuracy than NGO. Figure 4 It can be seen that ANGO calibration is faster.

[0203] Since the non-orthogonality error and non-alignment error of the inertial sensor have been obtained after static calibration, they can be used as a known error matrix and do not change with time and environment. Therefore, the dynamic calibration of the error parameters includes calibrating the accelerometer zero bias error and the accelerometer scale factor error; it also includes calibrating the gyroscope zero bias error and the gyroscope scale factor error.

[0204] The dynamic calibration objective function is consistent with the static calibration objective function, as shown in formulas (S-11) to (S-12);

[0205] The solution vector of the accelerometer dynamic calibration objective function is constructed using the accelerometer zero bias error parameter and the accelerometer scale factor error parameter, as shown in formula (S-27):

[0206]

[0207] Similarly, the solution vector of the gyroscope dynamic calibration objective function is constructed, as shown in formula (S-28):

[0208]

[0209] The static calibration error parameter value is loaded as the initial prior knowledge to obtain the next dynamic calibration prediction error parameter value, as shown in formula (S-29):

[0210]

[0211] In formula (S-29), represents the difference vector of the error parameter values ​​of the previous two calibrations, k represents the current kth dynamic calibration of the measuring point, X k-1,Best , X k-2,Best They represent the error parameter values ​​of the k-1th and k-2th calibrations respectively;

[0212] according to The k-th calibration error parameter value is predicted and generated, as shown in formula (S-30):

[0213]

[0214] The explosion phenomenon is described by explosively generating the k-th calibration algorithm initial population, as shown in formula (S-31):

[0215] X=normrnd(aver,sigma,m,n) (S-31)

[0216] In formula (S-31), the normrnd function represents the generation of a set of random numbers that conform to the normal distribution, which contains 4 parameters, among which: aver represents the mean of the normal distribution, that is, the value of each dimension of the prediction error parameter value, sigma represents the standard deviation of the normal distribution, which determines the range of the explosion. The larger the value, the larger the explosion range and the greater the probability of particles falling far away, and vice versa. m,n represents the generation of an m×n matrix;

[0217] After several experiments, the optimal values ​​are: sigma=0.2, m=1, n=6.

[0218] The constraint conditions include a convergence accuracy constraint and a maximum iteration number constraint.

[0219] Keep the initial error parameters unchanged as shown in Table 1, and set the second and third calibration error parameters as shown in Table 5 and Table 6 to simulate the dynamic calibration process.

[0220] Table 5 Accelerometer error parameter settings for the second measuring point

[0221]

[0222] Table 6 Accelerometer error parameter settings for the third measuring point

[0223]

[0224] The DANGO related parameter settings are shown in Table 7.

[0225] Table 7 DANGO initial parameter settings

[0226]

[0227] When the dynamic calibration reaches the same accuracy requirement as the static calibration, record the results of the DANGO iteration curve, ANGO iteration curve, and NGO iteration curve, such as Figure 4 As shown, the NGO curve represents the iterative curve of the traditional northern goshawk optimization algorithm during static calibration; ANGO represents the iterative curve of the northern goshawk optimization algorithm after static improvement; DANGO represents the iterative curve of the northern goshawk optimization algorithm after dynamic improvement. It can be seen that the dynamically improved algorithm has a significant improvement compared with the statically improved algorithm and the traditional algorithm, which proves the superiority of the dynamically improved northern goshawk optimization algorithm proposed in the present invention.

Claims

1. A method for calibrating inertial sensors in a measurement while drilling system. It is characterized in that The calibration method includes static calibration of error parameters of the inertial sensor, and dynamic calibration and error parameter update of the inertial sensor when drilling is stopped at a working measuring point of the measurement while drilling system; The specific steps of the static calibration of the error parameters include: S111: Establish an inertial sensor error model, including an accelerometer error model and a gyroscope error model, introduce gravity acceleration and the earth's rotation speed as calibration benchmarks, read accelerometer measurement data, gyroscope measurement data and calibration benchmark information, and establish a static calibration objective function; S112: Based on the objective function, analyzing the traditional Northern Goshawk optimization algorithm, and statically improving the traditional Northern Goshawk optimization algorithm, wherein the static improvement includes: The method of generating initial population by using Tent chaotic mapping is introduced to enhance the diversity of initial individuals in the solution space of static improved optimization algorithm. Adaptive factors and logarithmic functions are introduced to balance the relationship between global search and local search of the static improvement optimization algorithm, and to enhance the speed and accuracy of static improvement optimization algorithm in static calibration of error parameters. S113: setting initialization parameters of the static improved optimization algorithm, generating an initial population of the static improved optimization algorithm, and calculating an objective function value of the initial population of the static improved optimization algorithm; S114: Iteratively optimize the objective function value using the statically improved Northern Goshawk optimization algorithm until the objective function value satisfies the constraint condition, and output the static calibration error parameter value; The specific steps of the error parameter dynamic calibration and error parameter update include: S121: Establish an inertial sensor error model, including an accelerometer error model and a gyroscope error model, introduce gravity acceleration and the earth's rotation speed as calibration benchmarks, read accelerometer measurement data, gyroscope measurement data and calibration benchmark information, and establish a dynamic calibration objective function; S122: Based on the objective function, analyzing the traditional northern goshawk optimization algorithm, and dynamically improving the traditional northern goshawk optimization algorithm, wherein the dynamic improvement includes: Loading the static calibration error parameter value as initial prior knowledge to obtain the next dynamic calibration prediction error parameter value; Centering on the next dynamic calibration prediction error parameter value, simulating the explosion phenomenon, generating the initial population of the dynamic improvement optimization algorithm in the solution space, and calculating the objective function value of the dynamic improvement optimization algorithm; Adaptive factors and logarithmic functions are introduced to balance the relationship between global search and local search of the dynamic improvement optimization algorithm, and to enhance the speed and accuracy of dynamic calibration of error parameters of the dynamic improvement optimization algorithm. S123: setting initialization parameters of the dynamic improvement optimization algorithm, generating an initial population of the dynamic improvement optimization algorithm, and calculating an objective function value of the initial population of the dynamic improvement optimization algorithm; S124: Iteratively optimize the objective function value using the dynamically improved Northern Goshawk optimization algorithm until the objective function value satisfies the constraint condition, and output the dynamic calibration error parameter value; S125: using the dynamic calibration error parameter value to update the static calibration error parameter value; S126: Repeat S123 to S124, and use the dynamic calibration error parameter value to update the last calibration error parameter value, until a stop instruction is issued.

2. The inertial sensor calibration method in the measurement while drilling system according to claim 1, It is characterized in that The accelerometer error model is shown in formula (1): In formula (1), represents the accelerometer measurement value, represents the true value of the accelerometer, b a is the accelerometer bias error, A ab The transformation matrix from the carrier coordinate system b to the accelerometer coordinate system a is defined as follows: In formula (2), s and t are the intermediate temporary coordinate systems constructed, S at is the scale factor error matrix, as shown in formula (3); is a non-orthogonal error matrix, as shown in formula (4); is the misalignment error matrix, as shown in formula (5): In formula (3), They represent the accelerometer scale factor error parameters respectively; In formula (4), They represent the non-orthogonal error parameters of the accelerometer respectively; In formula (5), They represent the misalignment error parameters between the carrier coordinate system and the accelerometer coordinate system respectively.

3. The inertial sensor calibration method in the measurement while drilling system according to claim 2, It is characterized in that The gyroscope error model is shown in formula (6): In formula (6), represents the gyroscope measurement value, represents the true value of the gyroscope, b g represents the gyro bias error, A gb represents the transformation matrix from the carrier coordinate system b to the gyroscope coordinate system g, which is defined as formula (7): In formula (7), S gt is the scale factor error matrix, as shown in formula (8), is a non-orthogonal error matrix, as shown in formula (9), is the misalignment error matrix, as shown in formula (10); In formula (8), They represent the gyroscope scale factor error parameters respectively; In formula (9), They represent the gyroscope non-orthogonal error parameters respectively; In formula (10), represent the misalignment error parameters between the carrier coordinate system and the gyroscope coordinate system respectively.

4. The inertial sensor calibration method in the measurement while drilling system according to claim 3, It is characterized in that The static calibration objective function includes an accelerometer objective function and a gyroscope objective function. The accelerometer objective function is as shown in formula (11), and the gyroscope objective function is as shown in formula (12): In formulas (11) and (12), G represents the local gravity acceleration vector, ω represents the modulus of the earth's rotation velocity, l represents the number of collected data points, and n represents the total number of collected data points. and The calculations are as shown in formula (13) and formula (14); According to the misalignment error parameter between the carrier coordinate system and the accelerometer coordinate system, the accelerometer zero bias error parameter, the accelerometer scale factor error parameter and the accelerometer non-orthogonal error parameter, a solution vector of the accelerometer static calibration objective function is constructed, as shown in formula (15): According to the misalignment error parameter between the carrier coordinate system and the gyroscope coordinate system, the gyroscope zero bias error parameter, the gyroscope scale factor error parameter and the gyroscope non-orthogonal error parameter, a solution vector of the gyroscope dynamic calibration objective function is constructed, as shown in formula (16):

5. The inertial sensor calibration method in the measurement while drilling system according to claim 4, It is characterized in that The method of introducing Tent chaotic mapping to generate the initial population is as shown in formula (17): In formula (17), Z k represents the kth chaotic value, whose initial value is Z 0 ∈[0,1], β∈[0,1], the value of β affects the distribution of chaotic values; The adaptive factor and logarithmic function are introduced to balance the relationship between the global search and the local search of the static improved optimization algorithm, as shown in formula (18): In formula (18), w represents the adaptive factor.

6. The inertial sensor calibration method in the measurement while drilling system according to claim 5, It is characterized in that The error parameter dynamic calibration includes calibrating the accelerometer bias error and the accelerometer proportional factor error; and also includes calibrating the gyroscope bias error and the gyroscope proportional factor error.

7. The inertial sensor calibration method in the measurement while drilling system according to claim 6, It is characterized in that The dynamic calibration objective function is consistent with the static calibration objective function, as shown in equations (11) to (12); The solution vector of the accelerometer dynamic calibration objective function is constructed using the accelerometer zero bias error parameter and the accelerometer scale factor error parameter, as shown in formula (19): Similarly, the solution vector of the gyroscope dynamic calibration objective function is constructed, as shown in formula (20): The static calibration error parameter value is loaded as the initial prior knowledge to obtain the next dynamic calibration prediction error parameter value, as shown in formula (21): In formula (21), represents the difference vector of the error parameter values ​​of the previous two calibrations, k represents the current kth dynamic calibration of the measuring point, X k-1,Best , X k-2,Best They represent the error parameter values ​​of the k-1th and k-2th calibrations respectively; according to The k-th calibration error parameter value is predicted and generated as shown in formula (22):

8. The inertial sensor calibration method in the measurement while drilling system according to claim 7, It is characterized in that The explosion phenomenon is described by explosively generating the initial population of the k-th calibration algorithm, as shown in formula (23): X=normrnd(aver,sigma,m,n) (23) In formula (23), the normrnd function represents the generation of a set of random numbers that conform to the normal distribution, which contains 4 parameters, among which: aver represents the mean of the normal distribution, that is, the value of each dimension of the prediction error parameter value, sigma represents the standard deviation of the normal distribution, and m,n represents the generation of an m×n matrix.

9. The inertial sensor calibration method in the measurement while drilling system according to claim 8, It is characterized in that The constraint conditions include a convergence accuracy constraint and a maximum iteration number constraint.

Images