Parameter Identification of Pilot Model and Stability Analysis of Human-in-Loop Image Seeker

Abstract

In the human-in-loop (HIL) guidance mode, a pilot quickly identifies and flexibly locks on to a target through a real-time image signal transmitted by the aircraft. Then, the line-of-sight (LOS) angle error in the viewing field is tracked and compensated for in order to improve the guidance and control performance of the image-guided aircraft. Based on the physical structure and device parameters of the image seeker, an appropriate correction network is designed to improve the performance of the seeker stability loop. Aiming at a precise-extended crossover (PEC) pilot model, the structure of the dynamic model is optimized, and the maximum likelihood estimation (MLE) method of the output error structure is used to identify the dynamic parameters. This makes up for the deficiency of the existing modeling. In order to solve the nonlinear optimization problems encountered in the identification process, a hybrid strategy of a genetic algorithm (GA) and Gauss–Newton optimization algorithm is used to improve the probability of finding the global optimal solution. The simplex method is also used to improve the robustness of the algorithm. In addition, a hardware-in-the-loop simulation is designed and multi-round HIL experiment flow is performed. Moreover, based on the adaptability of the pilot to different image signal delays, the effects of different image signal delays on the stability and disturbance rejection rate (DRR) of the seeker control system are studied. The results demonstrate that the hybrid gradient optimization algorithm (HGOA) can find the global optimal value, and the identification model can accurately reflect the dynamic characteristics of the pilot. In the HIL guidance mode, the tracking compensation behavior of the pilot can reduce the influence of image signal delay on the disturbance of the aircraft body isolated by the seeker. The optimized PEC model and the identified dynamic parameters improve the efficiency of pilot training and screening.

1. Introduction

With the development of modern science and technology, image-guided aircraft have been highly valued for their features such as clear imaging, high accuracy, etc. [1,2,3,4,5]. As the core of image precision-guided aircraft, image guidance technology (IGT) performs object detection, identification and tracking, and provides accurate guidance information. With the rapid development of image processing, high-resolution sensors, chips and artificial intelligence, IGT has been widely studied and developed [6,7]. IGT can obtain rich imaging information with multi band, high sensitivity and high resolution, which helps in finding and identifying objects in a complex background [8,9].

The future environment is becoming more complex. In the case of long distance or similar interference, automatic target recognition (ATR) still has many technical difficulties [10,11]. Based on the recognition ability of the pilot, HIL recognition technology is suitable for object recognition under the interference of a complex environment and camouflage, and it has a high flexibility. The pilot can quickly identify and flexibly lock on to the object through the real-time image information sent back by the aircraft, and automatically track the object through the image tracker or manually through the control handle, which is referred to as HIL guidance [12,13]. Due to the limitations of quality and carrying space, optical fiber cannot be used to transmit image information over a longer distance [14]. On the other hand, data link communication using radio has the advantage of a longer transmission distance [15,16]. However, radio signals are susceptible to interference and obstruction, and the transmission and reception of image data require compression and decompression. This process will inevitably lead to image signal delay [17,18]. Therefore, the pilot cannot quickly and accurately identify and capture the object, and the phenomenon of pilot-induced oscillation (PIO) will also be potentially dangerous and paroxysmal [19,20,21].

The pilot model describes the process of the pilot controlling the optical axis of the seeker in order to track the object. When tracking manually, the pilot–seeker is a manual control system. The control loop of the image seeker with pilot dynamics will track and compensate for the error angle of the object in the field of view, so as to improve the control performance and guidance accuracy of the image-guided aircraft. Therefore, it is necessary to clarify the response of the pilot to the tracking error of the image seeker under different image signal delays in order to design the guidance and control system and improve the guidance accuracy. In the mid-20th century, researchers started to model the pilot. For instance, the PEC model [22,23] is the most representative and widely used model, which describes the dynamic characteristics of the pilot in a wider frequency range. However, in an actual control task, in order to better adapt to the dynamics of the control object and the image signal delay, the operation behavior of the pilot is mostly distributed in the higher frequency band, and it shows advanced characteristics. Afterwards, several studies on the identification of pilot model parameters used the PEC model, which contains different types of excitation signals [24]. The used identification methods included the MLE, GA and interval analysis method [25,26,27]. After the 1970s, modern control theory quickly developed. Kleinman proposed the optimal control model [28], and Hess developed the Hess structure model [19], which can accurately characterize the signal processing procedure of the pilot. In recent years, neural network and fuzzy control technologies have been frequently used in pilot modeling [29].

The parameters of the pilot model reflect the detailed characteristics of the pilot response, which can be obtained using experimental data identification. The parameter identification methods are divided into two types. The first one is the pilot response to time-invariant control objects, such as the Fourier coefficient method and MLE method [24,30]. The second one is the pilot response to time-varying control objects, such as the wavelet transform [13,31] and Kalman filtering [32] methods. The MLE method is a classical time-domain parameter identification method, which can identify model parameters directly using fitting time-domain data. Compared with frequency domain methods, time-domain parameter identification methods only use one step to fit the model, thereby reducing the bias and variance of parameter identification. In addition, the time-domain method has fewer design constraints on the excitation signal, and the frequency information of the excitation signal has less influence on the identifiability of the parameters. Zhao J.M. [33] used the auxiliary variable method to identify the parameters of the pilot model. However, he fixed the neural delay link of the pilot model and only tracked the tank target moving in a straight line. Therefore, the response of the pilot cannot be fully stimulated. Tang D.G. [34] studied the performance of the seeker loop and guidance loop in three modes: automatic tracking, manual tracking without delay and manual tracking with delay. However, the parameters of the pilot model remained unchanged in different control environments, while ignoring the adaptability of the pilot to the control environment. Potter [35] studied whether input shaping can improve the tracking performance of the pilot to a flexible system, in a manual control system with a time delay. However, the control object was limited to a low frequency and slightly damped vibration mode. In [36], the authors studied the modeling process of pilot adaptation to time-varying control dynamics in multi axis control tasks. However, the pilot model parameters were not obtained from the experimental data.

This paper focuses on the parameter identification of a pilot model and stability analysis of a human-in-loop image seeker, which mainly solves the following three problems. Firstly, based on a real image seeker model and parameters, the PEC model structure is optimized in order to solve the problem that the PEC model cannot accurately describe the pilot advanced compensation characteristics to image signal delay and the high-frequency distribution of operation behavior, compared to the work in [35]. Secondly, the pilot model dynamics parameters are identified using the HIL hardware-in-the-loop simulation system and experimental flow. A simulation experiment platform and wide-frequency random excitation signal are designed, and multiple rounds of HIL experiments are performed in order to tackle the problem of a single excitation signal and parameters that are not obtained from the parameter identification experiment, compared with the work in [33,36]. The MLE method of output error structure is used to identify the dynamic parameters, and the hybrid strategy of a GA and Gauss–Newton optimization algorithm is used to improve the probability of finding a global optimal solution. The simplex method is used to improve the robustness of the algorithm. Thirdly, the impacts of different image signal delays on the stability and DRR of the image seeker control system are studied, which solves the problem that the pilot dynamics are not considered in DRR research, compared with the work in [34].

The remainder of this paper is organized as follows. Section 2 introduces the image seeker physical model. Section 3 optimizes the PEC model structure, and proposes the method of dynamic parameter identification. Section 4 introduces the HIL hardware-in-the-loop simulation system and experimental flow. In Section 5, the impacts of different image signal delays on the stability and DRR of the seeker control system are studied. Finally, the conclusions are drawn in Section 6.

2. Image Seeker Physical Model

2.1. Seeker Physical Structure and Parameters

An image seeker usually uses an angular rate gyro platform seeker, whose control systems of pitch frame and yaw frame are independent and have a similar structure. The model representing the single channel control loop of the actual image seeker system, including the stabilization loop and tracking loop, is shown in Figure 1. As the inner loop of the tracking loop, the stabilization loop is used to isolate the angular disturbance of the aircraft, stabilize the optical axis of the seeker and improve the velocity and acceleration characteristics of the tracking system. The tracking loop controls the seeker’s optical axis movement according to the target position error information given by the pilot or image tracker under manual or automatic tracking mode, so as to perform the real-time tracking of the target [11,18,34].

When the seeker tracks the target, the image sensor is sensitive to the real-time LOS angle error. The control circuit synthesizes the LOS angle error and the feedback signal from the angular rate gyro in order to generate a control current to drive the torque motor. The latter drives the platform to rotate in order to eliminate the LOS angle error so that the seeker optical axis can track the LOS in real-time. The parameters of a certain platform seeker are given in

Table 1. Platform seeker parameters.

Parameter	$K_1$	$K_2$	L	R	$K_T$	J	${\omega }_n$	$\xi$	$K_E$
Value	13	22	0.007	8	0.24	0.002	85	0.4	0.17

[37].

2.2. Stabilization Loop Correction Network

Because of friction torque, damping torque and other disturbances between the bearings, the aircraft attitude change will affect the direction and output of the seeker in inertial space. In order to eliminate projectile motion coupling, the stabilization loop is usually designed with a higher open-loop gain and wider bandwidth. The stabilization loop, based on the established image seeker control loop, is shown in Figure 2.

Based on these seeker parameters, the hysteresis correction can be designed as follows [38]:

$$G_{\mathrm{c}}(s)=\frac{1+\beta Ts}{1+Ts}=\frac{1+0.91s}{1+1.44s}, \beta <1$$

(1)

The open-loop frequency responses of the stabilization loop before and after correction are presented in Figure 3, and the response parameters are shown in

Table 2. Open-loop frequency response parameters of the stabilization loop.

Parameter	$K$	$K_g$	$\gamma$	${\omega }_c$
Before correction	134	1	15	70
After correction	134	5	58	37

. It can be seen that the stability of the uncorrected system is low, while the addition of hysteresis correction attenuates the amplitude of the system in the high-frequency band, which reduces the amplitude crossing frequency ${\omega }_{\mathrm{c}}$ of the system. In addition, the appropriate selection of the hysteresis correction parameters highly improves the stability of the corrected system, and therefore ensures that the corrected system can meet the requirements of both the stability margin and the decoupling capacity index.

3. Pilot Model Dynamic Parameter Identification

In image-guided weapon systems, the pilot can manipulate the handle to track the target based on the LOS angle error displayed in the image screen, until it hits the target. Therefore, understanding and modelling the pilot’s response to the image display provides insights into the pilot perception process, and plays an important role in the evaluation of aircraft control systems and training simulation equipment.

3.1. Pilot Model Structure Optimization

McRuer describes pilot behavior using the PEC pilot model [22]. The frequency coverage range of the PEC model is extended to higher and lower frequencies, compared with the pilot crossover model whose frequency application range is limited to the vicinity of the amplitude crossing frequency [22]. The pilot response to the LOS angle error consists of two components: the response of the linear transfer function and the residual signal. The residual signal is usually considered a Gaussian signal which represents the non-linear part of the pilot model, because it is difficult to quantitatively measure its magnitude. This results in the inability of quantitatively studying control system metrics. Therefore, in this paper, the effect of the residual signal on the control system is not considered when analyzing the structure of the pilot model [37,39]. In the PEC model, the linear transfer function of the pilot response to the LOS angle error is given by:

$$G_{\mathrm{p}}(s)=K_{\mathrm{p}}\frac{{\tau }_{\mathrm{L}}s+1}{{\tau }_{l}s+1}\frac{{\omega }_{\mathrm{nm}}^2}{s^2+2{\xi }_{\mathrm{nm}}{\omega }_{\mathrm{nm}}s+{\omega }_{\mathrm{nm}}^2}{\mathrm{e}}^{-{\tau }_{\mathrm{d}}s}$$

(2)

As shown in Figure 4, the PEC model consists of a balance term and a physical constraint. The neuromuscular system is represented by the second-order oscillation link, whose dynamics limit the pilot’s response to tracking LOS angle error. The crossover theorem [22] proposed by McRuer shows that, in order to better adapt to the control object dynamics, the pilot will adjust the balance term dynamics so that the frequency characteristics of the system open-loop transfer function at the crossover frequency are close to those of a single integrator system. In actual control tasks, the dynamic form of the pilot’s control object is generally similar, which is mostly distributed in a higher frequency band, and the pilot will balance it near this band [39,40]. Moreover, when the pilot tracks the LOS angle error in the image screen, it has an adaptability to the image signal delay. More precisely, a longer image signal delay results in a lead compensation instead of a lag compensation from the pilot. Due to these two reasons, it is difficult for the lag time constant ${\tau }_{\mathrm{l}}$, representing the low-frequency characteristics, to accurately describe the real operation behavior of pilots. Furthermore, the authors in [37] show that it is hard to obtain through an identification experiment. Hence, the structure of the PEC model expressed in Equation (2) is optimized as follows:

$${G_{\mathrm{p}}}^{\prime }(s)=K_{\mathrm{p}}\left({\tau }_{\mathrm{L}}s+1\right)\frac{{\omega }_{\mathrm{nm}}^2}{s^2+2{\xi }_{\mathrm{nm}}{\omega }_{\mathrm{nm}}s+{\omega }_{\mathrm{nm}}^2}{\mathrm{e}}^{-{\tau }_{\mathrm{d}}s}$$

(3)

The following study is based on the pilot model shown in Equation (3), while the dynamic parameter vector to be identified is given by $\mathit{\theta }={\left[K_{\mathrm{p}},{\tau }_{\mathrm{L}},{\tau }_{\mathrm{d}},{\xi }_{\mathrm{nm}},{\omega }_{\mathrm{nm}}\right]}^{\mathrm{T}}$.

3.2. Dynamic Parameter Identification

Identifying the pilot model dynamic parameters from experimental measurement data is a highly nonlinear optimization problem, which results in an objective function containing many local minima besides the global minima. Compared to the frequency-domain parameter identification methods, time-domain parameter identification methods can identify the model parameters by directly fitting the time-domain data, while the excitation signal has fewer design constraints and the frequency information of the excitation signal has less influence on the parameter’s identifiability. MLE is a time-domain parameter estimation method that has been widely used in several fields [41]. It belongs to probabilistic Bayesian estimation methods. It consists of obtaining estimates of the model parameters by introducing the conditional probability density of the observations, constructing a likelihood function with the observations and the unknown parameters as independent variables, and taking the maximum probability of the observed value as the estimation criterion. In [37], the performance of the MLE algorithm is higher than that of the Fourier coefficient method in identifying the parameters of a pilot dynamic model.

The MLE method based on output error structure (OES) is used to improve identification accuracy. The hybrid gradient algorithm is used for optimization in order to deal with the problem that the gradient optimization algorithm of the parameter estimation process is sensitive to the initial values and prone to scattering. In other words, a GA is used to determine the initial value of the gradient algorithm. A Gauss–Newton algorithm is then used to speed up the operation. Finally, a simplex method is used to improve the robustness, which increases the probability of finding the global optima of a nonlinear optimization problem.

3.2.1. OES-Based MLE

Before applying the MLE method, the PEC model should be expressed in the form of state space. The pure delay link in the PEC model is a nonlinear link, which can be approximated by a high-order transfer function, such as the Padé approximation [39,42], expressed as follows:

$${\mathrm{e}}^{-{\tau }_{\mathrm{d}}s}=\frac{P(-{\tau }_{\mathrm{d}}s)}{P({\tau }_{\mathrm{d}}s)},P({\tau }_{\mathrm{d}}s)={\displaystyle \sum _{\widehat{k}=0}^{\lambda }\frac{(\lambda +\widehat{k})!}{\widehat{k}!(\lambda -\widehat{k})!}}{\tau }_{\mathrm{d}}{}^{-\widehat{k}}{s}^{\lambda -\widehat{k}}$$

(4)

where $\widehat{k}=0,\cdots ,\lambda$, $\lambda$ is the order of the approximation.

In order to ensure an accurate description of the pure delay link at high frequency, and avoid a high approximation of the PEC model, this paper uses a fifth-order Padé approximation [25]:

$${\mathrm{e}}^{-{\tau }_{\mathrm{d}}s}\approx \frac{-s^5+30\frac{1}{{\tau }_{\mathrm{d}}}{s}^4-420\frac{1}{{\tau }_{\mathrm{d}}^2}{s}^3+3360\frac{1}{{\tau }_{\mathrm{d}}^3}{s}^2-15120\frac{1}{{\tau }_{\mathrm{d}}^4}s+30240\frac{1}{{\tau }_{\mathrm{d}}^5}}{s^5+30\frac{1}{{\tau }_{\mathrm{d}}}{s}^4+420\frac{1}{{\tau }_{\mathrm{d}}^2}{s}^3+3360\frac{1}{{\tau }_{\mathrm{d}}^3}{s}^2+15120\frac{1}{{\tau }_{\mathrm{d}}^4}s+30240\frac{1}{{\tau }_{\mathrm{d}}^5}}$$

(5)

Substituting the expression of ${\mathrm{e}}^{-{\tau }_{\mathrm{d}}s}$ into Equation (3) results in the following:

$${G_{\mathrm{p}}}^{\prime }(s)\approx \frac{b_6{s}^6+b_5{s}^5+b_4{s}^4+b_3{s}^3+b_2{s}^2+b_{1}s+b_0}{s^7+a_6{s}^6+a_5{s}^5+a_4{s}^4+a_3{s}^3+a_2{s}^2+a_{1}s+a_0}$$

(6)

Due to the introduced Padé approximation, the coefficients in Equation (6) are nonlinear functions of the unknown identification parameters. For instance, the coefficients of the first and last terms of the denominator can be expressed as follows:

$$a_6=\frac{2{\xi }_{\mathrm{nm}}{\omega }_{\mathrm{nm}}{\tau }_{\mathrm{d}}^5+30{\tau }_{\mathrm{d}}^4}{{\tau }_{\mathrm{d}}^5},a_0=\frac{30240{\omega }_{\mathrm{nm}}^2}{{\tau }_{\mathrm{d}}^5}$$

(7)

The canonical form of the state space of Equation (6) is given by:

$$\dot{\mathit{x}}=\mathit{A}\left(\mathit{\theta }\right)\mathit{x}+\mathit{B}\left(\mathit{\theta }\right)\Delta q,u=\mathit{C}\left(\mathit{\theta }\right)\mathit{x}+n$$

(8)

where $\mathit{x}$ is the state vector of the system, $\mathit{A}\left(\mathit{\theta }\right)$ is the state matrix, $\mathit{B}\left(\mathit{\theta }\right)$ is the control matrix and $\mathit{C}\left(\mathit{\theta }\right)$ is the output matrix expressed as follows:

$$\mathit{A}(\mathit{\theta })=\left[\begin{array}{ccccc}0& 1& 0& \cdots & 0\\ 0& 0& 1& \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & 0\\ 0& 0& \cdots & 0& 1\\ -a_0& -a_1& -a_2& \cdots & -a_6\end{array}\right],\mathit{B}(\mathit{\theta })={\left[\begin{array}{ccccc}0& 0& 0& \cdots & 1\end{array}\right]}^T,\mathit{C}(\mathit{\theta })=\left[\begin{array}{ccccc}-b_0& -b_1& -b_2& \cdots & -b_6\end{array}\right]$$

(9)

Assuming that there is no influence of process noise, Equation (8) only contains an additional noise term n in the output equation. Therefore, the MLE of the OES method can be used to simplify the identification procedure. In the study of pilot model identification in [43], it is assumed that the residual signal of the pilot response is a zero-mean Gaussian noise, and this hypothesis was verified by experiments. Similarly, the residual signal n in Equation (8) is considered as a zero-mean Gaussian noise, and its mean and variance are defined as follows:

$$E\left\{\mathit{n}(k)\right\}=0,\hspace{1em}{\sigma }_n^2=E\left\{\mathit{n}(k){\mathit{n}}^T(k)\right\}$$

(10)

where $\mathit{n}(k)$ is the discrete sampling time series of the residual signal n, and $k=1,\cdots ,m$ represents the serial number, such that $m$ is the number of sampling points.

The MLE method requires finding an estimate of the model parameter vector which maximizes the likelihood function. The latter defines the conditional probability density function of the prediction error:

$$L(\mathit{\theta })=f\left(v(1),v(2),\dots ,v(k),\dots ,v(m)|\mathit{\theta }\right)$$

(11)

where the prediction error $v(k)$ is defined as the difference between the measured pilot output signal $u(k)$ and the PEC model output signal $\widehat{u}(k)$:

$$v(k)=u(k)-\widehat{u}(k)$$

(12)

According to the residual property defined in Equation (10), Equation (11) can be rewritten as follows:

$$L(\mathit{\theta })={\displaystyle \prod _{k=1}^m\frac{1}{\sqrt{2\mathsf{\pi }}{\sigma }_n}}{\mathrm{e}}^{-\frac{v^2(k)}{2{\sigma }_n^2}}={\left(2\mathsf{\pi }{\sigma }_n^2\right)}^{-\frac{m}{2}}{\mathrm{e}}^{-\frac{1}{2{\sigma }_n^2}{\displaystyle \sum _{k=1}^m{v}^2(k)}}$$

(13)

It can be seen from Figure 5 that the OES requires minimizing a quadratic penalty function related to the prediction error $v(k)$. Hence, the penalty function is defined as the match degree between the actual and simulated data, i.e., the difference between the actual estimates and the parameter estimates. The OES uses an iterative search algorithm in order to find the optimal parameter in the parameter space, and therefore minimize the penalty function.

When applying the MLE method, the model parameter vector estimation $\widehat{\mathit{\theta }}$ maximizes the likelihood function. On the other hand, using the MLE of OES requires obtaining the global minimum of the logarithmic likelihood function, thus simplifying it to a more direct optimization problem. The resulting parameter vector is the maximum likelihood estimate ${\widehat{\mathit{\theta }}}_{\mathbf{ML}}$, expressed as follows:

$${\widehat{\mathit{\theta }}}_{\mathbf{ML}}=\underset{\mathit{\theta }}{\arg }\min \left\{-\ln (L(\mathit{\theta }))\right\}=\underset{\mathit{\theta }}{\arg }\min \left\{\frac{m}{2}\ln (2\mathsf{\pi }{\sigma }_n^2)+\frac{1}{2{\sigma }_n^2}{\displaystyle \sum _{k=1}^m{v}^2(k)}\right\}$$

(14)

Equation (14) summarizes the parameter estimation problem for the PEC model, which defines a strongly nonlinear optimization problem. Similar to [39], the hybrid gradient optimization algorithm is then used to obtain the global optimal solution of this optimization problem.

3.2.2. Hybrid Gradient Optimization Algorithm

The traditional methods of solving the maximum likelihood minimum include the Levenberg–Marquardt, Gauss–Newton and Newton–Raphson methods. This paper proposes a hybrid gradient optimization algorithm (HGOA), which combines a genetic optimization algorithm and Gauss–Newton gradient optimization algorithm in order to estimate the parameters of the PEC model, as illustrated in Figure 6.

In the absence of a priori knowledge, the HGOA firstly applies the GA for the initial iteration, and then optimizes Equation (14) using the Gauss–Newton and simplex method algorithms. Because a large initial error may lead to scattered parameter estimates, the GA is used at the initial moment of the calculation in order to reduce the sensitivity of the algorithm to the initial value. In addition, Gauss–Newton has a fast convergence speed. However, it will produce inaccurate gradient information and lead to divergence when it is far from the real solution. Hence, in order to enhance the robustness of the algorithm, when the Gauss–Newton algorithm diverges, it is switched to the Nelder–Mead simplex method in the next iteration until the Gauss–Newton algorithm converges.

Due to the inherent randomness of a GA and the search for the optimal solution in the initial parameter vector set, the algorithm has a high probability of finding the global optimal value of the problem [44]. The GA first requires creating an initial population, i.e., the initial parameter vector set. The latter is randomly selected on the upper and lower bounds of the parameters, and it is represented by the following matrix:

$$\mathrm{Initial} \mathrm{population}=\left[\begin{array}{cccc}{\theta }_1^1& {\theta }_2^1& \cdots & {\theta }_{p_i}^1\\ {\theta }_1^2& {\theta }_2^2& \cdots & {\theta }_{p_i}^2\\ ⋮& ⋮& \cdots & ⋮\\ {\theta }_1^{p_j}& {\theta }_2^{p_j}& \cdots & {\theta }_{p_i}^{p_j}\end{array}\right]$$

(15)

where $p_i$ represents the number of model parameters, $p_j$ denotes the number of the parameter vector, element ${\theta }_{p_i}^{p_j}$ in each row represents the $p_i$ model parameter in the $p_j$ parameter vector, and each row of the population matrix represents the gene sequence of an individual.

The GA selects the individuals with higher fitness using gene manipulation. The population size should be large enough for each generation in order to have a high probability of searching for the global optimum.

The fitness of each individual in the population is assessed by the objective function presented in Equation (14). Figure 7 shows the process of generating a new population by applying three gene manipulations to each individual of the current population. The detailed steps are given by:

• Select: the fittest individuals are selected from the current population and copied to the next generation.
• Crossover: two individuals are randomly selected from the current population. The position of the crossover point in the gene is randomly selected to simulate the reproduction phenomenon in the evolution process, and therefore, the individuals with brand-new genes are obtained.
• Mutate: the current individual gene randomly selects the mutation point to change with a small probability, and then produces the new individual.

These steps are repeated until the termination condition is met. Because a GA uses probability search technology, the transfer mode and relationship from one search point to another are uncertain, which leads to inconsistent results in each iteration. Hence, the GA cannot be used only to estimate the PEC model parameters based on test data. The PEC model parameters estimated by the GA are very close to the global optimal solution. Furthermore, in order to obtain an accurate parameter estimation, the solution of the GA is used as the initial parameter estimation of the Gauss–Newton optimization algorithm. This gradient-based optimization algorithm is a classical method for solving MLE optimization problems [39,41]. The iterative equation of the Gauss–Newton optimization algorithm is given by:

$$\widehat{\mathit{\theta }}(i+1)=\widehat{\mathit{\theta }}(i)+\Delta \widehat{\mathit{\theta }},\hspace{1em}\Delta \widehat{\mathit{\theta }}=-{\mathit{M}}_{\mathit{\theta }=\mathit{\theta }(i)}^{-1}{\mathit{g}}_{\mathit{\theta }=\mathit{\theta }(i)}$$

(16)

where $\widehat{\mathit{\theta }}(i)$ is the iterative value of the parameter vector $\mathit{\theta }$ at step $i$, and ${\mathit{g}}_{\mathit{\theta }=\mathit{\theta }(i)}$ is the gradient of the likelihood function relative to the model parameters which can be obtained by calculating the Jacobian matrix of the model output relative to $\mathit{\theta }$.

The derivation process is presented in Appendix A.

In addition, to prevent the Gauss–Newton optimization algorithm from diverging and terminating the iteration, the algorithm is switched to the Nelder–Mead simplex method. The Nelder–Mead simplex method is a multi-dimensional direct search approach for local optimization. It adapts itself to the local landscape and contracts on to the final minimum. It is efficient and computationally compact. Its detailed process is provided in [45].

In this study, the residual variance of the PEC model is unknown. Its estimated value is provided in [46]:

$${\widehat{\sigma }}_n^2=\frac{1}{m}{\displaystyle \sum _{k=1}^m{v}^2(k)}$$

(17)

It can be seen that the variance estimate ${\widehat{\sigma }}_n^2$ is related to the model parameter vector $\mathit{\theta }$. For actual optimization, the relaxation technique is often used to optimize the objective function shown in Equation (14). It consists in alternately estimating the variance ${\sigma }_n^2$ and the model parameter vector $\mathit{\theta }$. More precisely, a variance estimate ${\widehat{\sigma }}_n^2$ is obtained based on a fixed model parameter vector $\mathit{\theta }$, by fixing ${\sigma }_n^2={\widehat{\sigma }}_n^2$. The variable to be optimized then becomes $\mathit{\theta }$. When the iteration step is very small, the variance estimate ${\widehat{\sigma }}_n^2$ is calculated based on the current model parameter vector $\mathit{\theta }$. Finally, these steps are repeated until the convergence criteria are met.

4. HIL Simulation Experiment

4.1. Experimental Setup

In HIL guidance, the pilot controls the inner loop of seeker, and its dynamic characteristics slightly change. In the parameter identification experiment of the pilot model, the pilot tracks and compensates for the error angle of the target in the field of view. In order to make the excitation signal have a wide frequency range and high signal-to-noise ratio, and show a certain randomness, the excitation signal was designed as the superposition of some sinusoidal signals [47]:

$$q_{\mathrm{t}}(t)={\displaystyle \sum _{j=1}^N{A}_j\sin ({\omega }_{j}t+{\phi }_j)}$$

(18)

In order to fully excite the pilot’s response, the spectrum of the excitation signal should be distributed halfway between high frequency and low frequency. Therefore, the frequency of the sinusoidal signal was designed to be evenly distributed in ${\omega }_j=0.1~20\mathrm{rad}/\mathrm{s}$, according to the logarithmic scale. The amplitude distribution was similar to a first-order low-pass filter.

The variance of the excitation signal was defined as the integral of the square of the sinusoidal signal amplitude. The relationship is shown in Equation (19), and the derivation process is provided in Appendix B.

$${\sigma }^2\left\{q_{\mathrm{t}}(t)\right\}=\frac{1}{2}{\displaystyle \sum _{j=1}^N{A}_j^2}$$

(19)

In addition, Damveld [40] mentioned that the phase has a small effect on the experimental results and was randomly selected. The excitation signal with a number of sinusoidal signals is shown in Figure 8. It can be seen that the maximum field angle is less than 15 degrees.

Figure 9 shows the hardware-in-the-loop simulation system and HIL experimental flow. The dynamic simulation module was developed on the basis of embedded system RTX, which collects the output of the handle and simulates the dynamic of the seeker system in real-time (cf. Figure 1), with a calculation frequency of 1000 Hz. The display module was developed using vista2D+MFC to present the LOS error angle through the UDP communication protocol. According to the error signal, the participants manipulate the handle to point the optical axis of the seeker to the target in order to reduce tracking error. It can be seen from the test time axis of Figure 9 that a preparation time of 10 s was set as the initial transient response time before the test in order to ensure that the participants could stably track the error signal. The test was divided into two stages:

• Training stage: the participants were trained for five rounds. The duration of each test round was 90 s. A fade-in time of 6 s and a fade-out time of 4 s was introduced to make the participants adapt to the task and turn into a linear control behavior.
• Test stage: it was performed under three conditions; no delay condition (phase I), small delay condition (phase II) and large delay condition (phase III). Phase II and phase III displayed the field error angle with an image signal delay of $\tau =250\mathrm{ms}$ and $\tau =500\mathrm{ms}$, respectively. As in the training stage, the participants also conducted five rounds of tests.

Six participants were selected to perform the HIL experiment. Participant A as a “new-hand”, did not carry out the preliminary training stage, but only carried out the test stage of five rounds. As an “old-hand”, participants B~F had five rounds of early training, and then five rounds of testing. During the test phase, a short rest period was provided to prevent the participants from becoming fatigued. After each test round, the participants were informed of their tracking score, so as to motivate them to improve their tracking performance during the initial familiarity period and maintain a stable performance after reaching a certain proficiency level. Mcruer [22] believes that participants need a transition process in order to adapt to changes in the dynamics of the controlled object. During the transition process, participants still maintain the original control mode, and then gradually adapt to the dynamic characteristics of the controlled object, and finally complete the compensation operation. In this experiment, the fade-in time provided sufficient time for the participants to make their control strategy adapt to the dynamics of the changed controlled object. That is, a delay link was added to the control loop.

4.2. Dynamic Parameter Selection

The signals of all the participants were recorded in the experiment. In order to reduce the signal noise influence and improve the accuracy of the pilot model parameter estimation, repeated experimental signals of all the participants were averaged. Taking “old-hand” participant B as an example, Figure 10 shows the LOS error angle signal, output signal, excitation signal (LOS angle) and tracking signal (optical axis space angle) for 25 s before the test. It can be seen from Figure 10c that participant B successfully aligned the optical axis with the target.

Figure 11 presents the amplitude spectrum of all the signals of Figure 10, where the frequency of the excitation signal shown in Figure 8a is indicated by a circle. It can be seen that, in the whole control loop, all the signals had a high signal-to-noise ratio at the excitation frequency, which met the requirements of the identification experiment.

The LOS error angle signal $\Delta q$ and pilot output signal $u$ in Figure 10 were used to identify the pilot model parameter vector $\mathit{\theta }={\left[K_{\mathrm{p}},{\tau }_{\mathrm{L}},{\tau }_{\mathrm{d}},{\xi }_{\mathrm{nm}},{\omega }_{\mathrm{nm}}\right]}^{\mathrm{T}}$ in Equation (3). The HGOA (cf. Figure 6) was used. The GA was iterated 100 times before the GN optimization algorithm was used. The upper and lower bounds of the GA parameter field were set as shown in

Table 3. Upper and lower bounds of the GA parameter field.

Parameter	$K_{\mathrm{p}}$	${\tau }_{\mathrm{L}}$	${\omega }_{\mathrm{nm}}$	${\xi }_{\mathrm{nm}}$	${\tau }_{\mathrm{d}}$
Lower bound	0	0	0	0	0
Upper bound	5	5	10	1	1

. The lower bound of the parameter field was set to 0, while the upper bound was set large enough to make the initial parameter field include all the possible solutions. In the process of GA optimization, the crossover probability and mutation probability of genes were set to 0.7 and 0.1, respectively. In order to simplify the calculation, the initial residual variance identification value ${\sigma }_n^2$ was set to 1. After the GA iteration, the result was used as the initial value for the GN optimization. If the iteration step of the parameter vector, objective function and residual variance identification value were smaller than the convergence index, the GN iterations were stopped.

GA is a search algorithm based on a probabilistic heuristic. The results of each iteration are not exactly the same under a certain iteration step constraint. The GA identification operation was repeated several times for the experimental results of each participant, and the eight solutions with the smallest iteration results were substituted into the GN algorithm for further optimization. Figure 12 shows all the iterative processes of these eight solutions for “old-hand” participant B. It can be seen that the final values given by the GA are different after eight repetitions. After performing the optimization using the GN algorithm, the same results were obtained. In other words, the value of the global minimum logarithmic likelihood function was obtained, and the corresponding model parameters were the identification results. Because the initial identification values of the residual variance in the two optimization algorithms were both set to 1, the obtained objective function values were large. Finally, the identification value of each algorithm was obtained through the previously mentioned alternate iteration process.

One hundred initial parameter vectors were randomly selected from the parameter domain in Table 3. The HGOA was used for optimization. The obtained results are presented in Figure 13, in ascending order. It can be seen that the HGOA can highly reduce the initial value selection, and accurately converge to the global minimum.

The test data of all the participants were collected for 80 s in the middle of the five rounds of the test stage. The parameter estimation value was first obtained using the MLE method. The global optimal solution was then determined using the HGOA, so as to obtain the dynamic parameter vector $\mathit{\theta }$. Figure 14 shows the difference in the mean value of the frequency response obtained by participant B, after five rounds of test stage. Under the same excitation signal and control object conditions, participant B had a very similar frequency response in the training stage and phase I stage. This shows that participant B had a good training effect and reached the control level of the “old-hand”. In addition, the frequency response of participant B had clear adaptive changes to the time delay. Its amplitude at low frequency decreased to ensure the stability margin of the seeker loop, and it provided the lead phase to compensate for the impact of the time delay. The identification results of all the participants are shown in

Table 4. Results of parameter identification for different subjects.

Dynamic Parameters	New-Hand	Old-Hand
	A	B	C	D	E	F
$K_{\mathrm{p}}$	0.18	0.32	0.32	0.27	0.29	0.31
${\tau }_{\mathrm{L}}/\mathrm{s}$	0.36	0.34	0.36	0.46	0.36	0.38
${\omega }_{\mathrm{nm}}/(\mathrm{rad}·{\mathrm{s}}^{-1})$	4.33	4.78	5.09	4.78	4.82	4.79
${\xi }_{\mathrm{nm}}$	0.26	0.15	0.27	0.31	0.27	0.26
${\tau }_{\mathrm{d}}/\mathrm{s}$	0.19	0.16	0.15	0.15	0.14	0.16

. Due to the individual differences of the participants and different control strategies, the experimental results were slightly different.

Figure 15 shows the Bode diagram of the identification model for all the participants. It can be seen that the amplitude frequency characteristics of participant A are significantly lower than those of other participants, which resulted in a lower tracking performance. In order to more truly and generally reflect the dynamic characteristics of the pilots, the identification results of participants B~F were selected for further analysis.

The correlation coefficient ${\rho }_{u\widehat{u}}$ is generally used to evaluate the accuracy of an identification model. It can be used to represent the ability of the identification model output to reproduce the measured output. If ${\rho }_{u\widehat{u}}$ is close to 1, the output of the identification model can well fit the measurement output, and the identification model can truly reflect the dynamic characteristics of the actual object. On the other hand, if ${\rho }_{u\widehat{u}}$ is close to 0, the identification result is poor and the dynamic characteristics of the actual object are not captured. ${\rho }_{u\widehat{u}}$ is expressed as follows [48]:

$${\rho }_{u\widehat{u}}=\frac{{\displaystyle \sum _{k=1}^m{\epsilon }_1{\epsilon }_2}}{\sqrt{{\displaystyle \sum _{k=1}^m{\epsilon }_1^2}}\sqrt{{\displaystyle \sum _{k=1}^m{\epsilon }_2^2}}},\hspace{1em}{\epsilon }_1=u(k)-\left(1/m\right){\displaystyle \sum _{k=1}^{m}u(k)},\hspace{1em}{\epsilon }_2=\widehat{u}(k)-\left(1/m\right){\displaystyle \sum _{k=1}^m\widehat{u}(k)}.$$

(20)

Figure 16 shows a comparison between the measured output of participant B and the PEC model output. It can be seen that the output of the PEC model can better fit the measured output.

Table 5. The ${\rho }_{u\widehat{u}}$ values of the PEC model.

Participant	B	C	D	E	F
${\rho }_{u\widehat{u}}/\%$	88.16	76.73	82.08	86.84	92.94

shows the ${\rho }_{u\widehat{u}}$ values of all the participants. It can be seen that the ${\rho }_{u\widehat{u}}$ values of the PEC model for participants B~F are greater than 70%. Therefore, the PEC identification model can better fit and reflect the dynamic characteristics of the pilots.

It can be seen from Figure 14 that trained pilots will adjust their control behavior in order to adapt to the dynamic changes of the controlled object caused by time delay. Simultaneously, it can be observed from Table 4 that differences in the identification results of trained pilots exist. The difference of pilot dynamics in [33] is mainly reflected in the visual perception lead time constant ${\tau }_{\mathrm{L}}$ and visual perception response delay ${\tau }_{\mathrm{d}}$. Moreover, the ${\tau }_{\mathrm{d}}$ of trained pilots is usually stable at almost 0.15 s. Thus, based on the identification results of participants B~F (cf. Table 4), the average value of $K_{\mathrm{p}},{\xi }_{\mathrm{nm}},{\omega }_{\mathrm{nm}}$ was considered and the value range of ${\tau }_{\mathrm{L}},{\tau }_{\mathrm{d}}$ was given, so as to represent the different control strategies and tracking compensation behaviors used by different pilots when facing the dynamic changes of the controlled object caused by delay. The values of the dynamic parameter vector $\mathit{\theta }$ are presented in

Table 6. Dynamic parameter settings.

Parameter	$K_{\mathrm{p}}$	${\tau }_{\mathrm{L}}/\mathrm{s}$	${\omega }_{\mathrm{nm}}/(\mathrm{rad}·{\mathrm{s}}^{-1})$	${\xi }_{\mathrm{nm}}$	${\tau }_{\mathrm{d}}/\mathrm{s}$
Value	0.3	0.34~0.46	4.85	0.25	0.1~0.2

.

It can be seen from the Bode diagram shown in Figure 14 and Figure 15 that the dynamic response of the new-hand and old-hand are quite different. The dynamic response of the new-hand trained by the simulator can be expressed by the PEC model shown in Equation (3). When the dynamic parameters of the new-hand meet Table 6, they can be judged to have a good operation performance and meet the conditions for becoming a professional pilot, so as to improve the efficiency of the pilot training and screening.

Figure 17 shows the open-loop frequency response (Bode diagram) of the seeker tracking loop, including the image signal delay, pilot link and stabilization loop. Figure 17 validates the crossing theorem proposed by McRuer. That is, in the amplitude frequency characteristic curve, the slope at the amplitude crossover frequency is almost $-20\mathrm{dB}/\mathrm{dec}$, while the slope at the high-frequency input area is less than $-20\mathrm{dB}/\mathrm{dec}$, because the participants have an adaptive advance adjustment behavior when dealing with image signal delay.

5. HIL Image Seeker Stability Analysis

5.1. Stability Domain Analysis

It can be seen from Figure 14 that participants will adjust their control behavior according to the existence of delay and the delay length, when the LOS angle error $\Delta q$ is displayed. For example, they will provide lead phase ${\tau }_{\mathrm{L}}$ to compensate for the delay impact, and adaptively adjust their own response delay ${\tau }_{\mathrm{d}}$, so as to improve the tracking performance.

Therefore, it is important to clarify the influence of the variation of the typical pilot model parameters on the stability of the image seeker under different image signal delays, for the design of the guidance control system and improvement of guidance accuracy. The influence of the pilot model parameters ${\tau }_{\mathrm{L}}$ and ${\tau }_{\mathrm{d}}$ on the stability of the seeker control system was studied according to the image seeker control loop shown in Figure 1. The closed-loop transfer function of the seeker stabilization loop is given by:

$${\Phi }_1(s)=\frac{K_2{K}_{\mathrm{T}}{G}_{\mathrm{c}}}{K_2{K}_{\mathrm{T}}{G}_{\mathrm{c}}{G}_{\mathrm{g}}+K_{\mathrm{T}}{K}_{\mathrm{E}}+Js(Ls+R)}$$

(21)

The dominant pole plays an important role in the transient response of the system. Therefore, by combining the lag correction network shown in Equation (1) and the platform seeker parameters shown in Table 1, Equation (21) is approximately equivalent to the following:

$${\Phi }_1(s)=\frac{213(s^2+427s+2.85\times {10}^5)}{(s+313)(s^2+166s+1.95\times {10}^5)}\approx \frac{1}{0.0032s+1}$$

(22)

The open-loop transfer function of the seeker tracking loop is given by:

$$G(s)=\frac{K_1{\mathrm{e}}^{-\tau s}{G_{\mathrm{p}}}^{\prime }(s){\Phi }_1(s)}{s}$$

(23)

The stability margin represents the relative stability of the control system, including the amplitude margin and phase margin. The phase margin $\gamma$ is expressed as follows:

$$\gamma =\mathsf{\pi }+\Phi ({\omega }_{\mathrm{c}})$$

(24)

The physical meaning of $\gamma$ is that, when $\Phi ({\omega }_{\mathrm{c}})$ lags behind the $\gamma$ angle, the system is in a critical stable state, and $\gamma =30°\sim 60°$ is generally considered.

The Nyquist stability criterion plays an important role in the frequency domain control theory, which can obtain the stability degree of the system. Equations (25) and (26), respectively, present the relationship between the image signal delay $\tau$ and the seeker tracking loop open-loop gain $K$ under the change of the visual perception lead time constant ${\tau }_{\mathrm{L}}$, and under the change of the visual perception response delay ${\tau }_{\mathrm{d}}$, when the phase margin is $\gamma$ and according to the Nyquist stability criterion. The derivation process is provided in Appendix C.

$$\omega =\frac{\pi /2-\gamma }{(\tau +1.75-{\tau }_{\mathrm{L}})}\hspace{1em}K\le \frac{\omega \sqrt{1-0.07{\omega }^2}}{\sqrt{1+{({\tau }_{\mathrm{L}}\omega )}^2}}$$

(25)

$$\omega =\frac{\pi /2-\gamma }{(\tau +1.2+{\tau }_{\mathrm{d}})}\hspace{1em}K\le \frac{\omega \sqrt{1-0.07{\omega }^2}}{\sqrt{1+{(0.4\omega )}^2}}$$

(26)

The stability region of the seeker under different parameter changes of the pilot model is shown in Figure 18. The stability region at special values is shown in

Table 7. Stability region of the seeker open-loop gain K.

$\tau$	$\gamma =40°,{\tau }_{\mathrm{d}}=0.15\mathrm{s}$			$\gamma =50°,{\tau }_{\mathrm{d}}=0.15\mathrm{s}$			$\gamma =40°,{\tau }_{\mathrm{L}}=0.4$			$\gamma =50°,{\tau }_{\mathrm{L}}=0.4$
	Visual perception lead time constant ${\tau }_{\mathrm{L}}$						Visual perception response delay ${\tau }_{\mathrm{d}}$
	0.34	0.4	0.46	0.34	0.4	0.46	0.1	0.15	0.2	0.1	0.15	0.2
0.1	0.56	0.58	0.59	0.45	0.47	0.49	0.60	0.58	0.56	0.48	0.47	0.45
0.2	0.53	0.54	0.56	0.43	0.44	0.45	0.56	0.54	0.53	0.45	0.44	0.43
0.3	0.50	0.51	0.53	0.40	0.41	0.43	0.53	0.51	0.50	0.43	0.41	0.40
0.4	0.47	0.48	0.50	0.38	0.39	0.40	0.50	0.48	0.47	0.40	0.39	0.38
0.5	0.45	0.46	0.47	0.36	0.37	0.38	0.47	0.46	0.45	0.38	0.37	0.36

. It can be seen that the larger the phase margin $\gamma$ (or the longer the image signal delay $\tau$), the smaller the stability region of the seeker, and therefore, the smaller the designed open-loop gain $K$. This reduces the seeker speed of target tracking. In addition, the smaller the visual perception lead time constant ${\tau }_{\mathrm{L}}$ (or the greater the visual perception response delay ${\tau }_{\mathrm{d}}$), the smaller the stability region of the seeker, and therefore, the smaller the designed open-loop gain $K$.

5.2. Disturbance Rejection Rate Analysis

The DRR of the seeker is an important index that characterizes its isolation ability to projectile disturbance. The greater the DRR, the lower the seeker’s ability to isolate the projectile disturbance. In order to achieve the stability of the seeker optical axis in the inertial space, the DRR should be reduced as much as possible in the seeker design. At present, the existing research on DRR [49,50,51,52,53,54] is limited to its amplitude percentage at the typical frequency, without considering the influence of the pilot dynamics in a real application. However, the DRR amplitude characteristics at different frequencies are generally concerned with engineering applications. In addition, the amplitude characteristics of the seeker DRR vary with the frequency change. Hence, the influence of the pilot dynamics on the DRR amplitude characteristics at different frequencies should be clarified in engineering applications.

The interference torque is generated by the relative rotation between the seeker frame and the base, which is the main factor leading to the DRR of platform seeker. It mainly includes spring torque and damping torque. The simplified model shown in Figure 19 is used to analyze the influence of the interference torque on the DRR and output accuracy of the seeker [49]. In the process of aircraft flight, the aircraft attitude disturbance is partially coupled to the seeker, so that the output LOS angular velocity has an error. The seeker DRR represents the decoupling ability of seeker to aircraft disturbance. It is defined as follows:

$$G_{GS}=\frac{\Delta \dot{q}}{\dot{\vartheta }}$$

(27)

where $\Delta \dot{q}$ is an additional component of the seeker output caused by projectile disturbance.

In an aircraft guidance system, due to the existence of DRR, the seeker outputs the wrong LOS angular velocity, which is introduced into the guidance system. It then generates an overload instruction according to the guidance law and passes it to the autopilot. The actuator on the aircraft deflects according to the control command, changes the aircraft attitude, and generates the aircraft motion. The motion passes through the seeker DRR model and causes the seeker to produce an error output. The block diagram of the platform seeker DRR under interference torque is shown in Figure 20. Considering the pilot model, the DRR transfer function caused by interference torque is given by:

$$G_{GS}(s)=\frac{(Ls+R)G_D{G}_{\mathrm{g}}s+K_{\mathrm{E}}{K}_{\mathrm{T}}{G}_{\mathrm{g}}s}{(Ls+R)J{s}^2+(Ls+R)G_{D}s+K_{\mathrm{E}}{K}_{\mathrm{T}}s+K_2{K}_{\mathrm{T}}{G}_{\mathrm{c}}{G}_{\mathrm{g}}s+K_1{K}_2{K}_{\mathrm{T}}{G}_{\mathrm{c}}{G_{\mathrm{p}}}^{\prime }{\mathrm{e}}^{-\tau s}}$$

(28)

The parameter settings are presented in

Table 8. Parameter settings.

Parameter	$K_{\omega }$	$K_{\mathrm{N}}$	${\tau }_{\mathrm{L}}$	${\tau }_{\mathrm{d}}$
Value	0.04	0.8	0.4	0.15

. Figure 21 shows the variation of the DRR amplitude near projectile frequency with different signal delays. It can be seen from Figure 21 that:

• When the pilot dynamics are not introduced into the system, the DRR is highly affected by the signal delay. In the vicinity of the projectile frequency, with the increase in signal delay, the DRR amplitude increases, resulting in the reduction of the seeker’s ability to isolate the projectile disturbance. In addition, the influence of signal delay on the DRR decreases with the increase in the projectile frequency. Therefore, when there is no pilot to track and compensate for the LOS error angle in the field of view, the large signal delay will reduce the stability of the seeker’s optical axis in the inertial space.
• When pilot dynamics are introduced into the system, the pilot will track and compensate for the LOS error angle in the field of view. The change of signal delay has a small influence on the DRR amplitude of the seeker.

6. Conclusions

The proposed optimized pilot model describes more accurately the pilot’s lead compensation characteristics for the image signal delay and high-frequency distribution of operation behavior. Multi-round HIL experiments were performed based on the simulation platform, and a wide frequency random excitation signal was designed to fully stimulate the pilot’s response. MLE was used to identify the dynamic parameters of the PEC model. In addition, the GA and HGOA were used to find the global optimal solution, while the simplex method was used to improve the robustness of the algorithm, so as to solve the nonlinear optimization problems encountered in the identification process. Moreover, the impacts of different image signal delays on the stability and DRR of image seeker control system were studied.

The simulation experiment shows that in the HIL tracking mode, the longer the image signal delay, the smaller the stability region of the seeker. This results in a slower seeker tracking speed and a larger error angle during stable tracking, which can easily cause the target to escape from the probe field. The adaptive control behavior of the pilot compensates for the phase lag caused by the signal delay to a certain extent. The smaller the visual perception lead time constant ${\tau }_{\mathrm{L}}$ (or the larger the visual perception response delay ${\tau }_{\mathrm{d}}$), the smaller the stability region of the seeker, and therefore, the slower the seeker tracking speed. When the pilot dynamics are introduced into the system, the pilot will track and compensate for the error angle of the probe field, and the change in image signal delay has a small impact on the DRR amplitude of the seeker. Furthermore, the optimized PEC model and identified dynamic parameters improve the efficiency of the pilot training and screening.