Optimization method for antiaircraft gun trajectory simulation and shooting parameter calculation process

Abstract

The invention relates to the field of antiaircraft guns, and provides an antiaircraft gun trajectory simulation and shooting parameter calculation process optimization method. The invention aims to solve the problems of low shooting precision, non-ideal hail suppression effect and low operation efficiency caused by the fact that an antiaircraft gun depends on manual aiming in manual hail suppression operation. According to the main scheme, the method comprises the following steps: 1) defining a mass point ballistic differential equation under a non-standard condition, including differential expressions of speed and position; 2) fitting a projectile air resistance law curve to obtain a resistance coefficient expression; 3) converting the target coordinate, the initial speed of the shell and the wind speed into a launching coordinate system; 4) carrying out ballistic integration by using an Adams forecast-correction method, and calculating the position and the speed of the shell; 5) predicting the future position of the target according to the trajectory data; (6) correcting the shooting angle, including judging whether a target is hit or not, and correcting the shooting angle if the target is not hit; and 7) controlling the number of iterations to optimize the calculation efficiency. The device is used for improving the shooting precision and the operation efficiency of the antiaircraft gun in the manual hail suppression operation.

Description

Technical Field

[0001] The invention relates to the field of anti-aircraft guns and provides an optimization method for anti-aircraft gun trajectory simulation and shooting parameter solution flow. Background Art

[0002] Hailstorms are one of the most common and highly destructive meteorological disasters in the world, often causing serious losses to agriculture, construction, power facilities and transportation systems. In order to cope with frequent hailstorms, artificial weather modification technology has gradually become an important means of hailstorm prevention and control. Among them, antiaircraft guns can effectively suppress the development of hail clouds due to their explosion effect, and have become a key tool for hail prevention operations. However, despite the large scale of hail prevention operations, their protection effect is still not ideal. The reasons for this problem are: on the one hand, the aging of technical personnel and the high mobility of personnel have led to insufficient technical level of hail prevention operations; on the other hand, the current antiaircraft gun operations mainly rely on manual aiming and shooting, which has problems of insufficient accuracy and poor environmental adaptability, which directly affects the hail prevention effect. Although the strong center azimuth and distance of the hail cloud can be read by radar, it becomes extremely difficult to visually judge the strong center of the cloud body and related operation information under complex weather conditions. Therefore, improving the accuracy of the antiaircraft gun trajectory and providing clear operation information for the gun point have become one of the key issues to improve the efficiency of artificial hail prevention operations.

[0003] Based on the above background, this paper proposes an optimization scheme for the trajectory simulation and shooting parameter solution process of anti-aircraft guns, aiming to improve the efficiency and accuracy of hail prevention operations. First, combined with the actual situation of current artificial hail prevention operations, the trajectory simulation is carried out by using a moderately simplified particle trajectory differential equation under non-standard conditions while ensuring the calculation accuracy and efficiency. Secondly, the 1943 resistance law curve is re-fitted with a polynomial using the least squares method to make it more consistent with the actual table value. In addition, this paper uses the fourth-order Runge-Kutta method and the Adams forecast-correction method to solve the ballistic differential equation, and verifies the calculation accuracy of the two methods by comparing with the measured trajectory data. Finally, according to the actual meteorological conditions of artificial hail prevention, the solution process of shooting parameters is optimized. Considering the influence of wind on the trajectory of anti-aircraft guns, the launch azimuth, elevation angle and ballistic flight time of anti-aircraft guns can be accurately solved at the same time. Through the above research, this paper provides a theoretical basis and practical reference for improving the scientificity and effectiveness of artificial hail prevention operations. Summary of the invention

[0004] The purpose of this paper is to solve the problem that the existing artificial hail prevention operations rely on manual aiming, resulting in low shooting accuracy, unsatisfactory hail prevention effects and low operating efficiency.

[0005] The present invention provides an optimization method for antiaircraft gun trajectory simulation and firing parameter solution process, which is characterized by comprising the following steps:

[0006] Step 1: Define the simplified nonstandard center-of-mass motion differential equation as:

[0007]

[0008]

[0009]

[0010]

[0011] where v x 、v y and v z are the speeds of the projectile in each direction of the coordinate axis, and w x 、w y and w z are the wind speeds in each direction, c is the ballistic coefficient, which depends on the shape, diameter and weight of the projectile, H(y) is the air density, G(v r , c s ) is the drag function, g represents the neutral acceleration;

[0012] Step 2: Fit the air resistance law curve of the projectile to obtain the resistance coefficient c s expression;

[0013] Step 3: According to the altitude angle θ and azimuth The target coordinates in the northeast sky coordinate system and the initial velocity of the projectile are v 0 and wind speed are converted into the launch coordinate system to obtain the coordinates of the target in the launch coordinate system (X L , Y L , Z L ), the initial velocity of the projectile v 0 The velocity component of the component in the emission coordinate system (v x0 , v yo , v z0 ), the wind speed component in the launch coordinate system (w x , w y , w z );

[0014] Step 4: Trajectory calculation: Use the Adams prediction-correction method to perform trajectory integration and obtain the position (x, y, z) and velocity (v x , v y , v z ), and determine the coordinates of the trajectory endpoint based on the trajectory data

[0015] Step 5: Based on the trajectory flight time T and the target's motion parameters, predict the target's position (x m,y m , z m );

[0016] Step 6: Shooting Angle Correction:

[0017] Step 6.1: Calculate the coordinates of the trajectory endpoint The distance between the predicted position and the target position (x m ,y m , z m ) is less than the preset threshold e d :

[0018] If D is less than the preset threshold e d , it is considered that the target is hit, and the current elevation angle θ is output i , azimuth and ballistic flight time T;

[0019] If D is greater than the preset threshold, then according to the trajectory end point coordinates and the target predicted position (x m ,y m , z m ), correct the current elevation angle θ i and azimuth Get the corrected elevation angle θ i+1 and azimuth

[0020] Where k is the correction factor, is the current shooting angle, is the altitude and azimuth corresponding to the end point of the trajectory, is the altitude and azimuth corresponding to the target's future position;

[0021] Step 6.2: Set the corrected elevation angle θ i+1 , azimuth Target predicted position (x m ,y m , z m ), return to step 4.

[0022] In the above scheme, the resistance coefficient c in step 2 is s The expression is as follows, where Ma is the Mach number:

[0023] c s =0.157(0≤Ma<0.7)

[0024] c s =-2.78195+11.74974Ma-15.66703Ma 2 +6.96900Ma 3 (0.7≤Ma<0.9)

[0025] c s =101.46937-321.91614Ma+339.45900Ma 2 -118.68720Ma 3 (0.9≤Ma<1)

[0026]

[0027] p 1 =0.154944442142445

[0028] p 2 =-2.64577102571359

[0029] p 3 =-0.389894040604623

[0030] p 4 =2.74288518588992

[0031] p 5 =0.381002068711428

[0032] p 6 =-1.49903292197855

[0033] p 7 =-0.237775531351884

[0034] p 8 =0.47141060301032

[0035] p 9 =0.114120378908769.

[0037] In the above scheme, in step 3, according to the azimuth and the elevation angle θ to convert the target's initial coordinates in the northeast sky coordinate system Use the transformation matrix to transform the target coordinates into the coordinates of the launch coordinate system (X L , Y L , Z L ).

[0038] Transformation Matrix:

[0039]

[0040] get:

[0041]

[0042]

[0043]

[0044] In the above scheme, in step 3, according to the azimuth The initial velocity component of the projectile is converted into the component in the launch coordinate system by using the height angle θ:

[0045]

[0046]

[0047]

[0048] In the above scheme, in step 3, according to the azimuth The conversion matrix of the wind speed in the northeast sky coordinate system to the wind speed in the launch coordinate system using the altitude angle θ is as follows:

[0049]

[0050] get:

[0051]

[0052] w y =wind U

[0053] .

[0055] In the above scheme, step 4 includes the following steps:

[0056] Step 4.1: Use the Runge-Kutta method to calculate the ballistic data within the time step Δt as the starting value of the Adams method, and calculate the coordinates and velocity values ​​of the first four steps.

[0057] Coordinates of the first four steps:

[0058] (x 1 ,y 1 , z 1 ), (x 2 ,y 2 , z 2 ), (x 3 ,y 3 , z 3 ), (x 4 ,y 4 , z 4 )

[0059] The speed values ​​for the first four steps are:

[0060] (v x1 , v y1 , v z1 ), (v x2 , v y2 , v z2 ), (v x3 , v y3 , v z3 ), (v x4 , v y4 , v z4 )

[0061] Step 4.1.1: According to the initial velocity v of the projectile 0 The velocity component of the component in the emission coordinate system (v x0 , v yo , v z0 ), wind speed component (w x , w y , w z ), time step Δt, calculate k 1 :

[0062] k 1x =Δt·(-cH(y)G(v r , c s )(v x0 -w x ))

[0063] k 1y =Δt·(-cH(y)G(v r , c s )(v y0 -w y )-g)

[0064] k 1z =Δt·(-cH(y)G(v r , c s )(v z0 -w z ))

[0065] Step 4.1.2: Calculate k 2 :

[0066]

[0067]

[0068]

[0069] Step 4.1.3: Calculate k 3 :

[0070]

[0071]

[0072]

[0073] Step 4.1.4: Calculate k 4 :

[0074] k 4x =Δt·(-cH(y)G(v r , c s )(v x0 +k 3x -w x ))

[0075] k 4y =Δt·(-cH(y)G(v r , c s )(v y0 +k 3y -w y )-g)

[0076] k 4z =Δt·(-cH(y)G(v r , c s )(v z0 +k 3z -w z ))

[0077] Step 4.1.5: Update position and velocity:

[0078]

[0079]

[0080]

[0081]

[0082]

[0083]

[0084] Repeat steps 4.1.2-4.1.5: to calculate the coordinate and velocity values ​​of the first four steps:

[0085] Use (x 1 ,y 1 , z 1 ) and (v x1 , v y1 , v z1) as the new initial condition, repeat steps 4.1.2-4.1.5, and calculate the coordinate and velocity values ​​(x 2 ,y 2 , z 2 ) and (v x2 , v y2 , v z2 );

[0086] Use (x 2 ,y 2 , z 2 ) and (v x2 , v y2 , v z2 ) as the new initial condition, repeat steps 4.1.2-4.1.5, and calculate the coordinate and velocity values ​​(x 3 ,y 3 , z 3 ) and (v x3 , v y3 , v z3 );

[0087] Use (x 3 ,y 3 , z 3 ) and (v x3 , v y3 , v z3 ) as the new initial condition, repeat steps 4.1.2-4.1.5, and calculate the coordinate and velocity values ​​(x 4 ,y 4 , z 4 ) and (v x4 , v y4 , v z4 );

[0088] Step 4.2: According to the Adams prediction formula, predict the next trajectory data (x n+1 ,y n+1 , z n+1 ) and (v x(n+1) , v y(n+1) , v z(n+1) );

[0089] Step 4.2.1: Get known data

[0090] Get the coordinates and speed values ​​of the first four steps:

[0091] Coordinates: (x 1 ,y 1 , z 1 ), (x 2 ,y 2 , z 2 ), (x 3,y 3 , z 3 ), (x 4 ,y 4 , z 4 )

[0092] Speed: (v x1 , v y1 , v z1 ), (v x2 , v y2 , v z2 ), (v x3 , v y3 , v z3 ), (v x4 , v y4 , v z4 )

[0093] Get the current time step: Δt;

[0094] Initialize n=4;

[0095] Step 4.2.2: Calculate the forecast value

[0096] Use the Adams forecast formula to calculate the coordinates and velocity forecast values ​​for the next step:

[0097] Coordinate prediction value:

[0098]

[0099]

[0100]

[0101] Speed ​​prediction value:

[0102]

[0103]

[0104]

[0105] Step 4.3: Correct the predicted value according to the Adams correction formula to obtain more accurate trajectory data;

[0106] Step 4.3.1: Substitute the Adams correction formula:

[0107] Substitute the corrected acceleration value calculated in step 1 into the Adams correction formula to calculate the corrected coordinate and velocity values:

[0108] Coordinate correction value:

[0109]

[0110]

[0111]

[0112] Speed ​​correction value:

[0113]

[0114]

[0115]

[0116] Step 4.3.2: Repeat steps 4.2 and 4.3 until the entire trajectory data calculation is completed;

[0117] Step 4.5: Determine the coordinates of the trajectory endpoint based on the trajectory data and speed

[0118] In the above scheme, step 5.1 includes: The correction formula can be expressed as:

[0119]

[0120] .

[0122] In step 6.2, the number of iterations is also determined before returning to step 4. If the number of iterations is greater than the threshold, the process ends. Otherwise, the corrected high and low angles θ are i+1 , azimuth Target predicted position (x m ,y m , z m ), return to step 4.

[0123] Because the present invention adopts the above technical means, it has the following beneficial effects:

[0124] 1. Simplification of particle trajectory differential equations under non-standard conditions: In view of the differences between actual conditions and standard conditions in anti-aircraft gun shooting, the particle trajectory differential equations under non-standard conditions are appropriately simplified to eliminate some complex parameters with less influence. This simplification method not only improves the calculation efficiency, but also ensures the accuracy of the calculation results, thus better adapting to the actual operating environment of artificial hail suppression operations.

[0125] 2. Refitting of the 1943 drag law curve: The 1943 drag law curve was compiled by the Soviet Union and given in the form of tabular values. In order to meet modern computing needs, it is essential to convert the tabular law into an analytical form. The data below Mach 1 were refitted by polynomial curve fitting based on the least squares method, which significantly improved the accuracy in the low Mach number region.

[0126] 3. Selection of numerical solution for trajectory: In trajectory calculation, the fourth-order Runge-Kutta method is combined with the Adams prediction-correction method to balance the calculation accuracy and efficiency. The specific analysis is as follows:

[0127] Problems with single technology:

[0128] Fourth-order Runge-Kutta method: The fourth-order Runge-Kutta method is an improved method based on Taylor series expansion, which is used to approximate the solution of differential equations. It solves differential equations by stepwise approximation, and the calculation formula is simple and easy to implement. However, in practical applications, the calculation efficiency is relatively low because the fourth-order right-hand side function needs to be calculated.

[0129] Adams forecast-correction method: The accuracy is greatly affected by initial conditions and the efficiency is high.

[0130] Combined effect:

[0131] The trajectory data of the first four steps are calculated by the fourth-order Runge-Kutta method as the starting value of the Adams prediction-correction method, which reduces the amount of calculation and improves the calculation efficiency.

[0132] By combining the prediction and correction steps, the Adams method can ensure a certain degree of calculation accuracy while ensuring calculation efficiency. Especially at a larger step size (0.04 seconds to 0.1 seconds), the error fluctuation of the Adams method is not large, and the performance is excellent.

[0133] 4. Optimization of the shooting parameter solution process: Among the traditional shooting parameter solution methods, the iterative correction method is widely used due to its faster convergence speed. However, this method has shortcomings in overall algorithm efficiency, and the errors in the high and low angle iterations are easily transferred to the position iterations. In order to solve this problem, the method was improved, the target was regarded as always in motion, and the correction of the target's high and low angles and coordinates was introduced simultaneously during the iteration process, which reduced the number of iteration layers and improved the calculation efficiency.

[0134] 5. Precise control of the iterative solution process: In order to ensure that the fire control solution can be completed within the specified time, an upper limit on the number of iterations is set. This precise control not only avoids the low calculation efficiency caused by too many iterations, but also ensures the accuracy of the calculation results and avoids unnecessary waste of computing resources. BRIEF DESCRIPTION OF THE DRAWINGS

[0135] Figure 1 It is a simplified flow chart of the present invention;

[0136] Figure 2 Fitting curve for the resistance law in 1943;

[0137] Figure 3 is the refitted curve;

[0138] Figure 4 is the error diagram of the fourth-order Runge-Kutta method;

[0139] Figure 5 This is the error map of the Adams forecast-correction method;

[0140] Figure 6 Comparison of height accuracy under different integration step sizes;

[0141] Figure 7 Comparison of distance accuracy under different integration step sizes;

[0142] Figure 8 Calculate the initial parameter table for the anti-aircraft gun firing parameters;

[0143] Fig. 9 Iteration process table

[0144] Fig.10 are the projectile trajectory and target trajectory in the last iteration. DETAILED DESCRIPTION

[0145] The following is a detailed description of the embodiments of the present invention. Although the present invention will be described and illustrated in conjunction with some specific embodiments, it should be noted that the present invention is not limited to these embodiments. On the contrary, modifications or equivalent substitutions made to the present invention should all be included in the scope of the claims of the present invention.

[0146] In addition, in order to better illustrate the present invention, numerous specific details are given in the following specific embodiments. It will be understood by those skilled in the art that the present invention can also be implemented without these specific details.

[0147] Fitting optimization of ballistic characteristics and drag law

[0148] Simplifying the differential motion equations of particle trajectory under non-standard conditions

[0149] In exterior ballistics, when studying the motion of the center of mass of an arrow, the effect of the arrow's rotation around the center of mass can be ignored, and the angle of attack is assumed to be zero. In addition, it is assumed that the shape and mass distribution of the arrow are symmetrical about the longitudinal axis, the surface is flat, the gravitational acceleration is constant and vertically downward, the Coriolis force is ignored, and the meteorological conditions are assumed to be standard windless and rainless. The center of mass motion equation established under these assumptions can reveal the basic laws of the center of mass motion of the arrow and simplify the calculation. However, in anti-aircraft gun shooting, the actual conditions are different from the standard conditions in the shooting table. In order to more accurately correct the errors under these non-standard conditions, it is necessary to establish the differential equation of the center of mass motion of the arrow that is applicable to non-standard conditions.

[0150] Since the range of anti-aircraft guns is short, the surface of the earth can be approximated as a plane, and the influence of the Coriolis force can be ignored, and the acceleration of gravity can also be regarded as a constant. Therefore, the differential equation of the center of mass motion of the projectile under non-standard conditions can be appropriately simplified. The simplified differential equation of the non-standard center of mass motion is:

[0151] Formula (1):

[0152]

[0153]

[0154]

[0155]

[0156]

[0157]

[0158]

[0159]

[0160] where v x 、v y and v z are the speeds of the projectile in each direction of the coordinate axis, and w x 、w y and w z are the wind speeds in each direction. c is the ballistic coefficient, which depends on the shape, diameter and weight of the projectile. H(y) is the air density, which is related to air pressure and virtual temperature. G(v r , c s ) is the drag function, which is related to the velocity v of the projectile relative to the air r And the standard elastic resistance coefficient c s Related, c s It can be found from the table of resistance law of 1943. τ(y) is the virtual temperature related to the height y.

[0161] Curve fitting of the law of air resistance of projectile

[0162] In exterior ballistics calculations, the projectile shape coefficient and the law of air resistance are usually used to describe the ballistic characteristics. However, the widely used 1943 law of air resistance for projectiles is given in tabular form. To meet modern computing needs, it is essential to convert the tabular law into an analytical form, which can be achieved through curve fitting.

[0163] Using 1stOpt mathematical software, the Macquart method and the general global optimization method were used to fit the 1943 table value law. By screening a set of fitting coefficients with the largest correlation coefficient, a rational expression suitable for calculation was obtained. The fitting results are shown in Figure 1 The maximum error of the curve occurs at Ma = 0.91, the absolute value of the error at this point is 3.3%, and the variance of the entire curve is 4.184×10 -6 Although the overall fitting effect is good, the fitting accuracy is still insufficient in the low and high Mach number regions. Therefore, we refit the data below Mach.

[0164] For data below Mach 1, a polynomial curve fitting based on the least squares method is used to refit the table data. The polynomial coefficients are calculated using the following formula:

[0165] W=(A T A) -1 A T B formula (2)

[0166] Among them, A is the transposed matrix of the Vandermonde matrix with the number of columns equal to the number of fitting table value data and the number of rows equal to the number of polynomial coefficients; B is the vector containing the table value data, and W is the obtained polynomial coefficient vector.

[0167] The final fitting results are as follows:

[0168] Formula (3):

[0169] c s =0.157(0≤Ma<0.7)

[0170] c s =-2.78195+11.74974Ma-15.66703Ma 2 +6.96900Ma 3 (0.7≤Ma<0.9)

[0171] c s =101.46937-321.91614Ma+339.45900Ma 2 -118.68720Ma 3(0.9≤Ma<1)

[0172]

[0173] p 1 =0.154944442142445

[0174] p 2 =-2.64577102571359

[0175] p 3 =-0.389894040604623

[0176] p 4 =2.74288518588992

[0177] p 5 =0.381002068711428

[0178] p 6 =-1.49903292197855

[0179] p 7 =-0.237775531351884

[0180] p 8 =0.47141060301032

[0181] p 9 =0.114120378908769

[0182] Among them, c s is the drag coefficient, Ma is the Mach number. The drag law fitting curve in 1943 is as follows Figure 2 The refitted curve is shown in Figure 3 As shown:

[0183] Comparative analysis shows that the accuracy of the new fitting curve in the low Mach number region is significantly improved. The maximum error point of the entire curve appears at Ma = 4.9, and the absolute value of the error at this point is 2.48%. The variance of the entire fitting curve is 1.67×10 -6 In the range of Mach numbers from 0.7 to 1, the fitting accuracy is significantly improved compared with previous results.

[0184] Numerical solution and accuracy analysis of trajectory

[0185] Overview of the ballistic solution

[0186] The core goal of the ballistic solution is to solve the nonlinear differential equations that describe the motion of the projectile. Since the equations involve complex factors such as air resistance, wind speed, and gravity, analytical solutions are usually difficult. Therefore, in the calculation of exterior ballistics, ballistic table solutions or numerical solutions are usually used.

[0187] The ballistic table method simplifies the calculation process in practical applications by pre-calculating and compiling tables, and is suitable for scenarios with relatively simple calculation conditions. In contrast, the numerical solution has become the main method of modern exterior ballistics due to its high accuracy and adaptability. The numerical solution solves differential equations by stepwise approximation. Commonly used algorithms include the Runge-Kutta method and the Adams prediction-correction method.

[0188] The Runge-Kutta method is an improved method based on Taylor series expansion, which is used to approximate the solution of differential equations. The most commonly used fourth-order Runge-Kutta method is widely used in ballistics. The calculation formula is as follows:

[0189]

[0190] If the value at n is known (t n ,y 1n ,y 2n , ..., y mn ), then the function value at n+1 can be obtained by this method. Its truncation error is O(h 5 ). The fourth-order Runge-Kutta method is widely used in trajectory calculations due to its high accuracy and flexibility, but its disadvantage is that the right-hand side function needs to be calculated four times in each step.

[0191] The Adams prediction-correction method only needs to calculate the right-hand side function once per step, which is more efficient. However, it depends on the initial conditions and needs to be self-started by the Runge-Kutta method. The calculation process of this method includes two steps: prediction and correction, which are implemented by the prediction formula and the correction formula respectively.

[0192] The prediction formula is used to predict the function value of the next step based on the solutions of the previous steps calculated by the Runge-Kutta method:

[0193]

[0194] The correction formula uses the latest function value to correct the forecast results, thereby improving the accuracy:

[0195]

[0196] By combining the forecast and correction steps, the Adams forecast-correction method can ensure a certain degree of calculation accuracy while ensuring calculation efficiency.

[0197] Accuracy verification of numerical solution

[0198] In order to verify the accuracy of the above numerical solution, this paper uses the ballistic data of Chongqing 152 Factory 37mm anti-aircraft gun provided by Li Zuxian (Li Zuxian, Liu Hongwu, Liao Jun, etc. Human shadow anti-aircraft gun operation safety assessment method based on external ballistic calculation. Meteorological Science and Technology, 2016, 44 (1): 152-156.) and others for analysis. The ballistic coefficient of the shell is 1.89. The data includes the measured results of height and horizontal distance. The ballistic fuze time is set from 8 seconds to 20 seconds, and the shooting angle ranges from 45° to 85°, involving a total of 63 sets of height data and 63 sets of horizontal distance data.

[0199] During the verification process, the azimuth angle was set to 0, the integration step was set to 0.05s, the initial velocity of the shell was 866m / s, the wind speed was set to 0, and other initial meteorological conditions were set according to the Chinese artillery meteorological standards. A total of 7 fuze durations were considered in this verification, and 9 elevation angles were set for each fuze duration. This application calculates the errors between the coordinates calculated by two numerical solutions (the fourth-order Runge-Kutta method and the Adams prediction-correction method) and the measured values. The absolute values ​​and average values ​​of the height and horizontal distance errors for each fuze duration are shown below. Figure 4 and Figure 5 As shown:

[0200] The results show that at all shooting angles and all fuse durations, the errors of the height and horizontal distance calculated by the two methods are less than 0.6%, indicating that under windless conditions, the two numerical solutions have high accuracy and can accurately simulate the exterior trajectory to meet the needs of artificial hail prevention operations.

[0201] Further comparison shows that under different fuze times, the average error of the Adams prediction-correction method is smaller than that of the fourth-order Runge-Kutta method, indicating that its accuracy is slightly better. To further verify this conclusion, this embodiment adjusts the integration step from 0.01 seconds to 0.1 seconds, and calculates the error average of 63 sets of data under different step lengths. The average absolute value of the height and horizontal distance error is as follows Figure 6 and Figure 7 shown.

[0202] As can be seen from the figure, when the step size is between 0.01 and 0.04 seconds, the errors of the two methods are similar, but as the step size increases to more than 0.04 seconds, the error of the Adams method does not fluctuate much, showing better performance than the fourth-order Runge-Kutta method. When the step size increases to 0.1 seconds, the error of the Adams method is lower than that of the Runge-Kutta method.

[0203] In summary, under a smaller step size (0.01 to 0.04 seconds), the accuracy of the two methods is comparable. Under a larger step size (0.04 to 0.1 seconds), the Adams method performs better, especially for long-term or real-time calculation scenarios. Therefore, the Adams prediction-correction method has certain advantages in ballistic calculations and is suitable for hail-proof anti-aircraft artillery fire control systems. By reasonably adjusting the step size, the calculation efficiency can be improved while ensuring accuracy.

[0204] Optimization of various calculation parameters and iterative process for anti-aircraft gun shooting

[0205] Shooting calculation

[0206] In fire control calculation, the hit point refers to the position where the antiaircraft artillery shell meets the target, which is usually determined by solving the intersection of the ballistic trajectory and the target trajectory. According to the uniqueness theorem of air trajectory, the ballistic coefficient, the firing angle (elevation angle) and the initial velocity of the shell can completely and uniquely determine a trajectory. If the target is in motion, the shooting lead time must also be considered. At this time, the aiming vector, the hit vector and the lead vector form a triangle, called the hit triangle.

[0207] The purpose of the shooting calculation is to determine the elevation angle, azimuth angle and flight time required for the projectile to hit the target. This process usually relies on the ballistic point differential equation and requires the target's spatial coordinate information. In ballistic calculations, the launch coordinate system OX is often used. L Y L Z L

[0208] The target coordinates are usually described in the Northeastern Sky coordinate system O-NEU. By appropriate rotation, the launch coordinate system can be converted to the Northeastern Sky coordinate system. Assume that the azimuth of the anti-aircraft gun launch is (In the northeast sky coordinate system, the angle between the launch direction and the N axis), then the transformation matrix from the launch coordinate system to the northeast sky coordinate system is as follows:

[0209]

[0210]

[0211]

[0212]

[0213] Iterative Correction Method of Shooting Parameters

[0214] Among the mainstream methods for solving shooting parameters, the iterative correction method is widely used because of its fast convergence speed. Li Qiang et al. (Li Qiang, Ouyang Pan. Antiaircraft gun fire control algorithm and simulation based on exterior ballistics. Firepower and Command Control, 2014, 39 (5): 131-134.) first regard the target position as a fixed point when solving the shooting parameters, and calculate it by iterating the height angle. This method first sets a height angle threshold and calculates the initial height angle according to the target position. In the launch coordinate system, when the horizontal coordinate of the projectile trajectory calculated by the Runge-Kutta method is consistent with the target horizontal coordinate, the height difference between the two on the vertical coordinate is calculated. If the height difference is less than the preset threshold, it is considered that the projectile hits the target. Subsequently, the flight time of the projectile is calculated, and the future position of the target is predicted based on this time. The iteration continues until the distance between the projectile and the target meets the requirements, thereby solving the shooting parameters of the moving target.

[0215] However, this method has shortcomings in overall algorithm efficiency, and the errors in the high and low angle iterations are easily transferred to the position iterations, which may lead to divergent results, especially in complex environments. To this end, this application improves the method by considering the target to be always in motion and introducing corrections to the target's high and low angles and coordinates during the iteration process, thereby reducing the number of iteration layers and improving computational efficiency.

[0216] Improve the solution process

[0217] Compared with traditional methods, the improved algorithm proposed in this application significantly improves iteration efficiency, reduces calculation time, and can more accurately and synchronously correct the elevation angle and target position, ensuring that the calculation results of shooting parameters are more reliable.

[0218] Furthermore, in the operation of artificial hail-proof artillery, the influence of wind on the trajectory cannot be ignored. Therefore, it is necessary to improve the above algorithm to retain the advantages of synchronous correction of elevation angle and target position, while introducing wind correction on trajectory. This application adopts a simplified non-standard conditional particle trajectory differential equation and uses the Adams prediction-correction method to solve the trajectory. The improved solution process is as follows Figure 1 As shown:

[0219] Assume that the initial coordinates of the target in the northeast sky coordinate system are The initial velocity of the projectile is v 0 , the initial values ​​of the altitude angle and azimuth angle are iteratively solved as and Since the ballistic differential equation is based on the launch coordinate system, it is necessary to convert the initial velocity components of the projectile into components in the launch coordinate system

[0220] The wind speed vector is defined based on the Northeast Sky coordinate system, with the positive direction of each axis as the positive direction of the wind. In order to consider the impact of wind on the trajectory, the wind speed based on the Northeast Sky coordinate system needs to be converted to the wind speed in the launch coordinate system. The conversion matrix from the wind speed in the Northeast Sky coordinate system to the wind speed in the launch coordinate system is as follows:

[0221]

[0222] By solving the ballistic differential equations, the coordinates of the projectile's trajectory endpoint can be obtained. and its flight time T. Based on the coordinates of the end point of the trajectory, the new elevation angle is calculated relative to the gun point. and azimuth Using the flight time T of the projectile, combined with the target's flight direction, speed and acceleration, the target's position after this period of time (x m ,y m , z m ), and calculate the corresponding new high and low angle θ m and azimuth

[0223] Next, calculate the coordinates of the end point of the trajectory and the target position (x m ,y m , z m ), if the distance D is less than a predetermined threshold value e d , it is determined that the target has been hit, and the current altitude, azimuth and ballistic flight time are output.

[0224] If the distance does not meet the threshold, the elevation and azimuth are corrected and the corrected values ​​are used to continue the iterative calculation. Assume that the elevation and azimuth input into the trajectory differential equation are θ i and The corrected elevation angle and azimuth

[0225] In the process of solving shooting parameters, it is unrealistic to use the solution error as the only sign of the end of iteration. Any iterative algorithm can obtain higher accuracy if the time allows and converges. However, once the iteration time exceeds the allowed limit, no matter how high the accuracy is, it is meaningless. Therefore, the above method sets two conditions for the end of iteration:

[0226] (1) If the calculated distance D is less than the predetermined threshold e d , then it is determined that the target is hit and the iteration ends naturally;

[0227] In order to avoid too many iterations due to too small precision setting or divergence in the solution process, and to ensure that the fire control solution can be completed within the specified time, an upper limit on the number of iterations is set. When the number of iterations reaches the set upper limit, the iterative solution process is terminated.

[0228] Verification and result analysis of shooting calculation

[0229] In order to verify the effectiveness of the above method, this paper sets a set of initial parameters for solution, see Figure 8 As shown in Table 1. This table lists the key parameters for simulating antiaircraft gun shooting, including wind speed, muzzle velocity, target position and speed, etc. By setting these parameters, this paper conducts detailed iterative calculations on the shooting parameters.

[0230] The parameter settings in Table 1 reflect the conditions of various wind speeds and target positions, simulating the antiaircraft gun firing process under different conditions. d When the distance is set to 0.5 meters, after 8 iterations, the final result is a shooting result with an elevation angle of 58.824° and an azimuth angle of 22.039°, a bullet flying time of 15.09 seconds, and a miss distance of only 0.32 meters. Fig. 9 As shown in Table 2.

[0231] As can be seen from Table 2, as the number of iterations increases, the difference between the elevation angle and the azimuth angle gradually decreases, and the missile-target distance gradually approaches the target position, indicating that the solution process gradually converges. After 8 iterations, the final solution of the elevation angle and azimuth angle has reached a high accuracy, and the miss distance is only 0.32 meters, which fully meets the accuracy requirements of hail prevention operations. Fig.10 Shown are the projectile and target trajectories from the final iteration.

[0232] It is worth noting that the number of iterations is significantly affected by the target's moving speed and acceleration. As the target's speed increases and the acceleration increases, the number of required iterations will also increase. From the data in the table, it can be seen that if a looser distance threshold is selected, the number of iterations will be significantly reduced. In addition, the choice of step size has an important impact on the convergence of the solution. If the step size is too long, the solution may not converge to the predetermined threshold; if the step size is too small, although the accuracy is improved, the calculation time will also increase. ] Therefore, the choice of step size must ensure convergence while taking into account computational efficiency.

Claims

1. A method for optimizing the trajectory simulation and firing parameter calculation process of an antiaircraft gun, characterized in that: The following steps are involved: Step 1: Define the simplified nonstandard center-of-mass motion differential equation as: where v x 、v y and v z are the speeds of the projectile in each direction of the coordinate axis, and w x 、w y and w z are the wind speed in each direction, c is the ballistic coefficient, which depends on the shape, diameter and weight of the projectile, H(y) is the air density, G(v r ,c s ) is the drag function, g represents the neutral acceleration; Step 2: Fit the air resistance law curve of the projectile to obtain the resistance coefficient c s expression; Step 3: According to the altitude angle θ and azimuth The target coordinates in the northeast sky coordinate system, the initial velocity of the projectile v0 and the wind speed are converted to the launch coordinate system to obtain the coordinates of the target in the launch coordinate system (X L ,Y L ,Z L ), the velocity component of the projectile's initial velocity v0 in the launch coordinate system (v x0 ,v yo ,v z0 ), the wind speed component in the launch coordinate system (w x ,w y ,w z ); Step 4: Trajectory calculation: Use the Adams prediction-correction method to perform trajectory integration and obtain the position (x, y, z) and velocity (v x ,v y ,v z ), and determine the coordinates of the trajectory endpoint based on the trajectory data Step 5: Based on the trajectory flight time T and the target's motion parameters, predict the target's position (x m ,y m ,z m ); Step 6: Shooting Angle Correction: Step 6.1: Calculate the coordinates of the trajectory endpoint The distance between the predicted position and the target position (x m ,y m ,z m ) is less than the preset threshold e d : If D is less than the preset threshold e d , it is considered that the target is hit, and the current elevation angle θ is output i , azimuth and ballistic flight time T; If D is greater than the preset threshold, then according to the coordinates of the end point of the trajectory and the target predicted position (x m ,y m ,z m ), correct the current elevation angle θ i and azimuth Get the corrected elevation angle θ i+1 and azimuth Where k is the correction factor, is the current shooting angle, is the altitude and azimuth corresponding to the end point of the trajectory, is the altitude and azimuth corresponding to the target's future position; Step 6.2: Set the corrected elevation angle θ i+1 , azimuth Target predicted position (x m ,y m ,z m ), return to step 4.

2. The method for optimizing the trajectory simulation and firing parameter calculation process of an antiaircraft gun according to claim 1, characterized in that: The drag coefficient c in step 2 s The expression is as follows, where Ma is the Mach number: c s =0.157(0≤Ma<0.7) <h2 style=";text-align:left;direction:ltr">c<h2 style=";text-align:left;direction:ltr"> s <h2 style=";text-align:left;direction:ltr"> =-2.78195+11.74974Ma-15.66703Ma<h2 style=";text-align:left;direction:ltr"> 2 <h2 style=";text-align:left;direction:ltr"> +6.96900Ma<h2 style=";text-align:left;direction:ltr"> 3 <h2 style=";text-align:left;direction:ltr"> (0.7≤Ma<0.9) c s =101.46937-321.91614Ma+339.45900Ma 2 -118.68720Ma 3 (0.9≤Ma<1) p1＝0.154944442142445 p2＝-2.64577102571359 p3＝-0.389894040604623 p4＝2.74288518588992 p5＝0.381002068711428 p6＝-1.49903292197855 p7＝-0.237775531351884 p8＝0.47141060301032 p9＝0.114120378908769。 3. The method for optimizing the trajectory simulation and firing parameter calculation process of an antiaircraft gun according to claim 1, characterized in that: In step 3, according to the azimuth and the elevation angle θ to convert the target's initial coordinates in the northeast sky coordinate system Use the transformation matrix to transform the target coordinates into the coordinates of the launch coordinate system (X L ,Y L ,Z L ). Transformation Matrix: get:

4. The method for optimizing the trajectory simulation and firing parameter calculation process of an antiaircraft gun according to claim 1, characterized in that: In step 3, according to the azimuth The initial velocity component of the projectile is converted into the component in the launch coordinate system by using the height angle θ: v yo =v0sinθ 5. The method for optimizing the antiaircraft gun trajectory simulation and firing parameter calculation process according to claim 1, characterized in that: In step 3, according to the azimuth The conversion matrix of the wind speed in the northeast sky coordinate system to the wind speed in the launch coordinate system using the altitude angle θ is as follows: get: w y =wind U 。 6. The method for optimizing the trajectory simulation and firing parameter calculation process of an antiaircraft gun according to claim 1, characterized in that: Step 4 includes the following steps: Step 4.1: Use the Runge-Kutta method to calculate the trajectory data within the time step Δt as the starting value of the Adams method, and calculate the coordinates and velocity values ​​of the first four steps. Coordinates of the first four steps: (x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) The speed values ​​for the first four steps are: (v x1 ,v y1 ,v z1 ),(v x2 ,v y2 ,v z2 ),(v x3 ,v y3 ,v z3 ),(v x4 ,v y4 ,v z4 ) Step 4.1.1: Based on the velocity component (v0) of the projectile's initial velocity in the launch coordinate system x0 ,v yo ,v z0 ), wind speed component (w x ,w x ,w z ), time step Δt, calculate k1: k 1x =Δt·(-cH(y)G(v r ,c s )(v x0 -w x )) k 1y =Δt·(-cH(y)G(v r ,c s )(v y0 -w y )-g) k 1z =Δt·(-cH(y)G(v r ,c s )(v z0 -w z )) Step 4.1.2: Calculate k2: Step 4.1.3: Calculate k3: Step 4.1.4: Calculate k4: k 4x =Δt·(-cH(y)G(v r ,c s )(v x0 +k 3x -w x )) k 4y =Δt·(-cH(y)G(v r ,c s )(v y0 +k 3y -w y )-g) k 4z =Δt·(-cH(y)G(v r ,c s )(v z0 +k 3z -w z )) Step 4.1.5: Update position and velocity: Repeat steps 4.1.2-4.1.5: to calculate the coordinate and velocity values ​​of the first four steps: Using (x1, y1, z1) and (v x1 ,v y1 ,v z1 ) as the new initial condition, repeat steps 4.1.2-4.1.5, and calculate the coordinate and velocity values ​​(x2, y2, z2) and 9v of the second step x2 ,v y2 ,v z2 ); Using (x2, y2, z2) and (v x2 ,v y2 ,v z2 ) as the new initial condition, repeat steps 4.1.2-4.1.5, and calculate the coordinate and velocity values ​​(x3, y3, z3) and 9v of the third step x3 ,v y3 ,v z3 ); Using (x3,y3,z3) and (v x3 ,v y3 ,v z3 ) as the new initial condition, repeat steps 4.1.2-4.1.5, and calculate the coordinate and velocity values ​​(x4, y4, z4) and (v x4 ,v y4 ,v z4 ); Step 4.2: According to the Adams prediction formula, predict the next trajectory data (x n+1 ,y n+1 ,z n+1 )and (v x(n+1) ,v y(n+1) ,v z(n+1) ); Step 4.2.1: Get known data Get the coordinates and speed values ​​of the first four steps: Coordinates: (x1,y1,z1),(x2,y2,z2),(x3,y3,z3),(x4,y4,z4) Speed: (v x1 , v y1 , v z1 ), (v x2 , v y2 , v z2 ), (v x3 , v y3 , v z3 ), (v x4 , v y4 , v z4 ) Get the current time step: Δt; Initialize n=4; Step 4.2.2: Calculate the forecast value Use the Adams forecast formula to calculate the coordinates and velocity forecast values ​​for the next step: Coordinate prediction value: Speed ​​prediction value: Step 4.3: Correct the predicted value according to the Adams correction formula to obtain more accurate trajectory data; Step 4.3.1: Substitute the Adams correction formula: Substitute the corrected acceleration value calculated in step 1 into the Adams correction formula to calculate the corrected coordinates and Speed ​​value: Coordinate correction value: Speed ​​correction value: Step 4.3.2: Repeat steps 4.2 and 4.3 until the entire trajectory data calculation is completed; Step 4.5: Determine the coordinates of the trajectory endpoint based on the trajectory data and speed 7. The method for optimizing the trajectory simulation and firing parameter calculation process of an antiaircraft gun according to claim 1, characterized in that: Step 5.1 includes: The correction formula can be expressed as:

8. The method for optimizing the trajectory simulation and firing parameter calculation process of an antiaircraft gun according to claim 1, characterized in that: In step 6.2, the number of iterations is also determined before returning to step 4. If the number of iterations is greater than the threshold, the process ends. Otherwise, the corrected high and low angles θ are i+1 , azimuth Target predicted position (x m ,y m ,z m ), return to step 4.

Images