*Article* **IN-ME Position Error Compensation Algorithm for the Near-Field Beamforming of UAVs**

**Yinan Zhang, Guangxue Wang, Yi Leng \*, Guowen Yu and Shirui Peng \***

Department of Information Countermeasure, Air Force Early Warning Academy, Wuhan 430010, China **\*** Correspondence: lengyi200209@163.com (Y.L.); psr99@21.cn (S.P.)

**Abstract:** The target of an unmanned aerial vehicle swarm will present near-field characteristics when it is integrated as an array, and the existence of the unmanned aerial vehicle swarm motion error will greatly deteriorate the beam pattern formed by the array. To solve these problems, a near-field array beamforming model with array element position error is constructed, and the Taylor expansion of the phase difference function is used to approximately simplify the model. The improved Newton maximum entropy algorithm is proposed to estimate and compensate for the phase errors. The maximum entropy objective function is established, and the Newton iterative algorithm is used to estimate the phase error iteratively. To select the proper Newton iteration initial value, based on a single reference source signal, the initial value of the phase error is estimated through the phase gradient information of the received array signal. Beamforming is carried out after phase error compensation regarding the array. In order to assess the mismatch of the phase error compensation function based on the proposed method, when the beam is scanning, the effectively compensated spatial area of the array beamforming is divided, which lays a foundation for subsequent spatial region division and unmanned aerial vehicle swarm path planning. The simulation results show that the beam formed by the method proposed in this paper has a lower sidelobe level, and as the signal-to-noise ratio changes, the robustness of the proposed method is better validated. The proposed algorithm can effectively suppress the adverse influence of array element position error on array beamforming, and when the beam is scanning, the effectively compensated area of the phase error compensation function is divided, based on the proposed method.

**Keywords:** unmanned aerial vehicle swarm; antenna array; near-field beamforming; array element position error compensation; spatial area division

**MSC:** 94-08

#### **1. Introduction**

With the development of beamforming technology and unmanned aerial vehicle swarm (UAVs) technology, the antenna array elements equipping the UAVs are combined into a distributed beamforming system through collaboration among the swarm to achieve directional high gain signal transmission. Therefore, using beamforming technology to improve the combat capability of UAVs in a fierce confrontation environment is a growing trend [1]. It is obvious that taking UAVs as an array has many advantages, such as using low-cost UAVs to achieve high gain signal transmission, which is also known as the beam pattern, but the targets of the array are often located in the near-field, and the motion error [2–5] of a single UAV as an array element will deteriorate the beam pattern formed by the array. At present, the array position error calibration is mainly founded on the calibration method of three or more known radiant sources, based on far-field signals, or the rotating array calibration method, based on a single known as a radiant source [6–8], but the latter is, in essence, also a multi-radiant source calibration method exchanging time for space. Reference [2] proposes using three simultaneously existing auxiliary sources

**Citation:** Zhang, Y.; Wang, G.; Leng, Y.; Yu, G.; Peng, S. IN-ME Position Error Compensation Algorithm for the Near-Field Beamforming of UAVs. *Mathematics* **2022**, *10*, 3256. https://doi.org/10.3390/ math10183256

Academic Editor: Paolo Crippa

Received: 19 July 2022 Accepted: 30 August 2022 Published: 7 September 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

to calibrate the position error of the array, but this method can only use the Taylor firstorder expansion of the array steering vector with error when the position error has a very small disturbance, which means the position error cannot be corrected by this method when it is too large. In array beamforming, the beam cannot be formed when an array element position error exists, but as long as the phase error caused by the position error is compensated for, the beam can be formed effectively. It is easier to solve a one-dimensional phase error than to solve two- or three- dimensional position errors, and the phase error can be obtained by using a single known radiant source as a reference source [9–12]. As for the solution to phase error, based on the information of the reference source, the phase error can be estimated using the phase gradient (PG) value of the received signal between the array elements [13–17], but this method requires a high signal-to-noise ratio (SNR), and under the condition of low SNR, the deviation of the phase error estimation value from the real value is large. Phase error can also be solved by iteration. The authors of [18–21] take the entropy of the image as the objective function and estimate the phase error by using the Newton iterative method, but this method presents problems, such as the selection of the initial value affecting the algorithm time, causing the algorithm to fall into local optimization. When the expected main-lobe direction of beamforming is the same as the reference source direction, the phase error in main-lobe direction can be compensated for by the estimated phase error compensation function, but when it is not the same as the reference source direction, the mismatch of the phase error compensation function occurs, which will affect the beamforming pattern. Therefore, it is necessary to study the influence of phase error compensation function mismatch, which is the same as the division of the area effectively compensated by the estimated phase error compensation function. The authors of [22] propose that when the near-field radial compensation filter is used to filter out interference sources, the beamforming pattern will deteriorate when the sources are far away from the array, but the study does not specifically analyze the pattern deterioration level. The work in [23] proposes a method to form a multi-focus antenna array in the near-field, with specified amplitude and phase conditions for the multi-focus problem of the near-field array, but this study does not involve the analysis of the mismatch of the single-focus filter.

Based on the above analysis, firstly, a near-field array beamforming model with UAV position error is constructed in Section 2, and the model is approximately simplified by the Taylor expansion of the phase difference function of the near-field signal. Secondly, the improved Newton maximum entropy (IN-ME) algorithm is used to estimate and correct the phase error in Section 3. Under the condition of low SNR, taking the single known radiant source as a reference source, the initial value of the phase error is estimated through the PG information of the received array signal, and then, taking the maximum entropy function of the array synthetic power at the reference source as the objective function, the Newton iterative solution to the phase error is carried out. After the phase error is compensated, beamforming is carried out. Thirdly, when phase error estimation based on proposed IN-ME algorithm is used to compensate for beamforming, if the array forms beam scanning, the phase compensation function will be mismatched. Based on the phase error compensation function of the IN-ME algorithm, the effectively compensated area which can be compensated for to effectively form the beam pattern is further divided in Section 4, laying a foundation for the subsequent spatial region division and UAV path planning. The simulation results in Section 5 show the effectiveness of the algorithm, and the conclusions are given in Section 6. The adverse influence of UAV position errors on near-field array beamforming is effectively suppressed, and the effectively compensated area, with phase compensation function mismatch, is divided.

#### **2. Near-Field Beamforming Model**

#### *2.1. Array Signal Model*

Under the condition of targets being in the near-field, the assumption that the electromagnetic wave is the plane wave is not tenable, and it should be modeled as a spherical

wave. As shown in Figure 1, it is assumed that the near-field antenna array of *N* UAVs which can be seen as *N* elements is linearly distributed along the axis *x* at equal intervals, the total length of the array is *L*, the coordinate of the *n*th element is (*xn*, 0), where, *<sup>n</sup>* <sup>=</sup> 1, 2, ··· , *<sup>N</sup>*, and <sup>−</sup> *<sup>L</sup>* <sup>2</sup> <sup>≤</sup> *xn* <sup>≤</sup> *<sup>L</sup>* <sup>2</sup> , *rn* refers to the distance from the *n*th antenna array element to the source *P*, *R* refers to the distance from the source *P* to the array center *O*, *θ* is the included angle between *OP* and the axis *y*, and when the direction of *OP* projection on the *x* axis is the same as positive direction of the *x* axis, *θ* is positive, and vice versa. The coordinate of source *P* is *P*(*xP*, *yP*).

**Figure 1.** Near-field array model.

For any radiant source *P* in the near-field, according to the spherical wave theory, the radiant signal received by the array element can be expressed as

$$s(\mathbf{x}\_n) = \frac{A}{r\_n} \exp\left[\mathbf{j}\left(2\pi ft - 2\pi \frac{\Delta r\_n}{\lambda}\right)\right] = \frac{\tilde{A}}{r\_n} \exp\left(-\mathbf{j}k\Delta r\_n\right) \tag{1}$$

where, *<sup>A</sup>*<sup>7</sup> <sup>=</sup> *<sup>A</sup>*exp(*j*2*<sup>π</sup> f t*) is the radiant signal of the source *<sup>P</sup>*, j is the imaginary unit, *<sup>f</sup>* is the radiant signal frequency, *λ* is the wavelength, *k* = <sup>2</sup>*<sup>π</sup> <sup>λ</sup>* is the wavenumber, and Δ*rn* is the wave path difference from the source *P* to the *n*th antenna array element, with the array center *O* as the reference point.

**Theorem 1.** ([24]) *Suppose that <sup>r</sup>*(*x*) *is a differentiable function for which <sup>r</sup>*(*a*), ··· ,*r*(*n*)(*a*) *all exist, and there is neighborhood U*(*a*) *of* a *, in which the function value of r*(*x*) *at x* = a *can be written as r*(*x*) = *<sup>r</sup>*(*a*) + *<sup>r</sup>*(*a*)(*<sup>x</sup>* <sup>−</sup> *<sup>a</sup>*) <sup>+</sup> *<sup>r</sup>*(*a*)(*x*−*a*) 2 2! + *o* (*x* − *a*) 2 .

Theorem 1 is also known as Taylor's theorem. For the antenna array composed of UAVs, although the targets of interest in wireless communication, electronic reconnaissance, and jamming are generally located in the near-field of the array, the distances from the targets of interest to the antenna center are usually much greater than the aperture of the antenna. Therefore, it can be assumed that *R L* ≥ 2|*xn*|. Then, the Taylor series expansion of Δ*rn* can be obtained by

$$
\Delta r\_n = r\_n - R = \sqrt{(\mathbf{x}\_n - \mathbf{x}\_P)^2 + (y\_P)^2} - R \approx -\sin\theta \cdot \mathbf{x}\_n + \frac{\cos^2\theta}{2R} \mathbf{x}\_n^2 \tag{2}
$$

Through Equation (2), the phase difference between the received signal of the *n*th array element and the received signal at the array center *O* is

$$
\Delta \boldsymbol{\varphi}\_n = -\frac{2\pi}{\lambda} \Delta \boldsymbol{r}\_n \approx -k \left( -\sin \theta \cdot \mathbf{x}\_n + \frac{\cos^2 \theta}{2R} \mathbf{x}\_n^2 \right) \tag{3}
$$

By substituting Equation (3) into Equation (1) and ignoring the influence of distance on signal amplitude, the received signal model of the *n*th array element in the near-field can be approximately written as

$$s(\mathbf{x}\_{\rm tr}) \approx \tilde{A} \exp\left(-\mathrm{jk}\left(-\sin\theta \cdot \mathbf{x}\_{\rm tr} + \frac{\cos^2\theta}{2R} \mathbf{x}\_{\rm tr}^2\right)\right) \tag{4}$$

For *M* radiant sources in the near-field, the radiant signal amplitude of the *m*th radiant source is *Am*, the frequency is *fm*, the polar coordinate is (*θm*, *Rm*), and the distance from the *n*th array element to the *m*th radiant source is *rnm*, where, *m* = 1, ··· , *M*. Then, the received signal of the array is

$$\mathbf{X} = \sum\_{m=1}^{M} s\_m(\boldsymbol{x}\_n) + \boldsymbol{V} = \boldsymbol{\tilde{A}}\boldsymbol{A}(\boldsymbol{\theta}, \boldsymbol{R}) + \boldsymbol{V} \tag{5}$$

where, *sm*(*xn*) <sup>=</sup> *<sup>A</sup>*7*m*exp −*jk* <sup>−</sup> sin *<sup>θ</sup><sup>m</sup>* · *xn* <sup>+</sup> cos<sup>2</sup> *<sup>θ</sup><sup>m</sup>* <sup>2</sup>*Rm <sup>x</sup>*<sup>2</sup> *n* ,*A*7*<sup>m</sup>* <sup>=</sup> *Am*exp(*j*2*<sup>π</sup> fmt*), *<sup>X</sup>* <sup>=</sup>

[*x*<sup>1</sup> *x*<sup>2</sup> ··· *xN*] *<sup>T</sup>* is a *<sup>N</sup>* <sup>×</sup> 1 dimensional receiving data vector, *<sup>A</sup>*(*θ*, *<sup>R</sup>*) is a *<sup>N</sup>* <sup>×</sup> *<sup>M</sup>* dimensional array steering vector, *<sup>S</sup>* is a *<sup>M</sup>* <sup>×</sup> 1 dimensional signal vector, and *<sup>V</sup>* is a *<sup>N</sup>* <sup>×</sup> <sup>1</sup> dimensional noise vector. There is

$$A(\theta, R) = \begin{bmatrix} a(\theta\_1, R\_1) & a(\theta\_2, R\_2) & \cdots & a(\theta\_M, R\_M) \end{bmatrix} \tag{6}$$

$$\begin{aligned} \mathbf{a}(\theta\_{\mathsf{m}}, \mathsf{R}\_{\mathsf{m}}) &= \begin{bmatrix} a\_{1}(\theta\_{\mathsf{m}}, \mathsf{R}\_{\mathsf{m}}) & a\_{2}(\theta\_{\mathsf{m}}, \mathsf{R}\_{\mathsf{m}}) & \cdots & a\_{N}(\theta\_{\mathsf{m}}, \mathsf{R}\_{\mathsf{m}}) \end{bmatrix}^{\mathsf{T}} \\\\ \mathbf{b} &= \begin{bmatrix} \begin{pmatrix} -\operatorname{jk}\left(-\sin\theta\_{\mathsf{m}} \,\boldsymbol{x}\_{1} + \frac{\cos^{2}\theta\_{\mathsf{m}}}{2\mathsf{M}\_{\mathsf{m}}} \boldsymbol{x}\_{1}^{2}\right) \end{bmatrix}^{\mathsf{T}} & \begin{pmatrix} -\operatorname{jk}\left(-\sin\theta\_{\mathsf{m}} \,\boldsymbol{x}\_{2} + \frac{\cos^{2}\theta\_{\mathsf{m}}}{2\mathsf{M}\_{\mathsf{m}}} \boldsymbol{x}\_{2}^{2}\right) \end{bmatrix} & \begin{pmatrix} \begin{pmatrix} \cos\theta\_{\mathsf{m}} \,\boldsymbol{x}\_{1} + \frac{\cos^{2}\theta\_{\mathsf{m}}}{2\mathsf{M}\_{\mathsf{m}}} \boldsymbol{x}\_{2}^{2} \end{pmatrix} & \begin{pmatrix} \cos\theta\_{\mathsf{m}} \,\boldsymbol{x}\_{3} + \frac{\cos^{2}\theta\_{\mathsf{m}}}{2\mathsf{M}\_{\mathsf{m}}} \boldsymbol{x}\_{3}^{2} \end{pmatrix} \end{bmatrix} \end{aligned} \tag{7}$$
 
$$V = \begin{bmatrix} v\_{1} & v\_{2} & \cdots & v\_{N} \end{bmatrix}^{\mathsf{T}} \tag{8}$$

where, [·] <sup>T</sup> represents matrix transpose, and the *<sup>N</sup>* <sup>×</sup> 1 dimensional vector *<sup>a</sup>*(*θm*, *Rm*) is the steering vector of the array under the *m*th radiant source.

The received signals of each array element are weighted and summed to obtain the array output as, which is beamforming

$$\mathbf{Y} = \mathbf{W}^{\mathrm{H}} \mathbf{X} = \mathbf{W}^{\mathrm{H}} \left( \tilde{A} \mathbf{A} (\theta, \mathbf{R}) + \mathbf{V} \right) \tag{9}$$

where, [·] <sup>H</sup> represents the conjugate transpose of the matrix, and *W* = *w*<sup>1</sup> *w*<sup>2</sup> ··· *wN* T is the *N* × 1 dimensional weighted vector of the array. When the polar coordinate of the expected signal is (*θP*, *RP*), and the distance from the *n*th array element is *rnP*, there is *wn* <sup>=</sup> exp(−j*k*Δ*rnP*), where, <sup>Δ</sup>*rnP* <sup>=</sup> <sup>−</sup> sin *<sup>θ</sup><sup>P</sup>* · *xn* <sup>+</sup> cos2 *<sup>θ</sup><sup>P</sup>* <sup>2</sup>*RP <sup>x</sup>*<sup>2</sup> *n*.

#### *2.2. Actual Near-Field Beamforming*

When there are errors in the positions of the array elements, the actual coordinate of the *n*th array element is (*x*ˆ*n*, *y*ˆ*n*), where, *x*ˆ*<sup>n</sup>* = *xn* + Δ*xn*, and *y*ˆ*<sup>n</sup>* = *yn* + Δ*yn*. Under the model of this paper, *yn* = 0. Similar to *2.1.*, the modified actual wave path difference Δ*r*ˆ*<sup>n</sup>* in the Taylor series expansion is

$$\begin{split} \Delta \hat{r}\_n &\approx -\sin\theta \hat{x}\_n - \cos\theta \hat{y}\_n + \frac{\cos^2\theta}{2R} \hat{x}\_n^2 \\ &+ \frac{\sin^2\theta}{2R} \hat{y}\_n^2 + \frac{(\cos\theta + \sin\theta)}{2R^2} \hat{x}\_n \hat{y}\_n \end{split} \tag{10}$$

By using Equation (10) to approximately simplify *s*(*x*ˆ*n*), the actual *n*th array signal model can be expressed as

$$s(\pounds\_{\mathfrak{n}}) \approx \tilde{A} \exp \left( -\mathrm{jk} \begin{pmatrix} -\sin \theta \pounds\_{\mathfrak{n}} - \cos \theta \pounds\_{\mathfrak{n}} + \frac{\cos^{2} \theta}{2R} \pounds\_{\mathfrak{n}}^{2} \\ + \frac{\sin^{2} \theta}{2R} \pounds\_{\mathfrak{n}}^{2} + \frac{(\cos \theta + \sin \theta)}{2R^{2}} \pounds\_{\mathfrak{n}} \pounds\_{\mathfrak{n}} \end{pmatrix} \right) \tag{11}$$

Then the actual received signal of the array is

$$\hat{\mathbf{X}} = \sum\_{m=1}^{M} s\_m(\hat{\mathbf{x}}\_n) + \mathbf{V} = \tilde{A}\hat{\mathbf{A}}(\theta, \mathbf{R}) + \mathbf{V} \tag{12}$$

where, *X*ˆ = *<sup>X</sup>*<sup>ˆ</sup> <sup>1</sup> *<sup>X</sup>*<sup>ˆ</sup> <sup>2</sup> ··· *<sup>X</sup>*<sup>ˆ</sup> *<sup>N</sup> <sup>T</sup>* is an actual *<sup>N</sup>* <sup>×</sup> 1 dimensional received data vector, *sm*(*x*ˆ*n*) is the actual *<sup>m</sup>*th, *<sup>m</sup>* <sup>=</sup> 1, ··· , *<sup>M</sup>* is the received radiant source signal, *<sup>A</sup>*<sup>ˆ</sup> (*θ*, *<sup>R</sup>*) is an actual *N* × *M* dimensional array steering vector, and

$$s\_{\mathfrak{m}}(\mathbf{\hat{x}}\_{\mathfrak{n}}) = \mathbf{s}\_{\mathfrak{m}}(\mathbf{x}\_{\mathfrak{n}}) \cdot w\_{\mathfrak{m}} = \mathbf{s}\_{\mathfrak{m}}(\mathbf{x}\_{\mathfrak{n}}) \cdot \exp\left(-\mathrm{jk}\begin{pmatrix} -\sin\theta\_{\mathfrak{m}} \cdot \Delta \mathbf{x}\_{\mathfrak{n}} - \cos\theta\_{\mathfrak{m}} \cdot \Delta y\_{\mathfrak{n}} + \frac{\cos^{2}\theta\_{\mathfrak{m}}}{R\_{\mathfrak{m}}} \mathbf{x}\_{\mathfrak{n}} \Delta \mathbf{x}\_{\mathfrak{n}} + \frac{\cos^{2}\theta\_{\mathfrak{m}}}{2R\_{\mathfrak{m}}} \Delta \mathbf{x}\_{\mathfrak{n}} \\\\ + \frac{\sin^{2}\theta\_{\mathfrak{m}}}{2R\_{\mathfrak{m}}} \Delta y\_{\mathfrak{n}}^{2} + \frac{(\cos\theta\_{\mathfrak{m}} + \sin\theta\_{\mathfrak{m}})}{2R\_{\mathfrak{m}}^{2}} (\mathbf{x}\_{\mathfrak{n}} + \Delta \mathbf{x}\_{\mathfrak{n}}) \Delta y\_{\mathfrak{n}} \end{pmatrix}\right) \tag{13}$$

$$\hat{A}(\theta, R) = \begin{bmatrix} \hat{\mathfrak{a}}(\theta\_1, R\_1) & \hat{\mathfrak{a}}(\theta\_2, R\_2) & \dots & \hat{\mathfrak{a}}(\theta\_M, R\_M) \end{bmatrix} \tag{14}$$

$$\operatorname{Ad}(\theta\_{m\prime}, \mathbb{R}\_m) = [\mathbb{A}\_1(\theta\_{m\prime}, \mathbb{R}\_m), \mathbb{A}\_2(\theta\_{m\prime}, \mathbb{R}\_m) \quad \cdots \quad \mathbb{A}\_N(\theta\_{m\prime}, \mathbb{R}\_m)]^T$$

$$=\left[\exp\left(-\mathrm{jk}\left(\begin{array}{c} -\sin\theta\_{\mathrm{m}}\hat{\mathbf{x}}\_{1} - \cos\theta\_{\mathrm{m}}\hat{\mathbf{y}}\_{1} + \frac{\cos^{2}\theta\_{\mathrm{m}}}{2R\_{\mathrm{m}}}\hat{\mathbf{x}}\_{1}^{2} \\ + \frac{\sin^{2}\theta\_{\mathrm{m}}}{2R\_{\mathrm{m}}}\hat{\mathbf{y}}\_{1}^{2} \frac{(\cos\theta\_{\mathrm{m}}\sin\theta\_{\mathrm{m}})}{2R\_{\mathrm{m}}^{2}}\hat{\mathbf{x}}\_{1}\boldsymbol{\varrho}\_{1} \end{array}\right)\right)\dots\exp\left(-\mathrm{jk}\left(\begin{array}{c} -\sin\theta\_{\mathrm{m}}\boldsymbol{\varrho}\_{N} - \cos\theta\_{\mathrm{m}}\boldsymbol{\varrho}\_{N}\frac{\cos^{2}\theta\_{\mathrm{m}}}{2R\_{\mathrm{m}}}\boldsymbol{\varrho}\_{N}^{2} \\ \frac{\sin^{2}\theta\_{\mathrm{m}}}{2R\_{\mathrm{m}}}\hat{\boldsymbol{\varrho}}\_{N}^{2} \frac{(\cos\theta\_{\mathrm{m}}\sin\theta\_{\mathrm{m}})}{2R\_{\mathrm{m}}^{2}}\boldsymbol{\varrho}\_{N}\boldsymbol{\varrho}\_{N} \end{array}\right)\right)\right)^{\mathrm{T}}\tag{15}$$

where, *wen* is the error value of the *n*th array element. When *M* = 1, the actual beamforming formula can be written as

$$\hat{Y} = \mathbf{W}^{\text{H}} \hat{\mathbf{X}} = \sum\_{n=1}^{N} w\_n^\* \cdot (w\_{\text{en}} s(\mathbf{x}\_n) + v\_{\text{\text{\textell}}}) \tag{16}$$

where, *w*∗ *<sup>n</sup>* is the conjugate complex of *wn*.

Obviously, due to the presence of array position errors and noises, it is impossible to effectively carry out near-field beamforming by using the signal weight that can be obtained when the array is at the nominal value to compensate for the actual signal.

#### **3. IN-ME Phase Error Compensation Algorithm Based on Reference Source**

Only by compensating for phase error can the array effectively form a beam pattern. We propose the IN-ME algorithm, which, compared to Newton maximum entropy (N-ME) algorithm, can obtain higher convergent value within a shorter period, as well as form a lower sidelobe-level beam pattern with a higher SNR.

#### *3.1. Phase Error Estimation Based on Newton Maximum Entropy Algorithm*

Entropy is a measurement of information uncertainty. The maximum entropy principle is a criterion of probabilistic model learning, and when the probability distribution of random variables obeys uniform distribution, the entropy reaches its highest level. As for beamforming, we expect that after phase error compensation, the synthetic powers in the main-lobe direction of the array's received signals in all snapshots increase, which becomes a multi-objective optimization problem. If we regard the ratio of synthetic power in the main-lobe direction of the array's received signals in one snapshot to the sum of the synthetic powers in all snapshots as the probability distribution, compared to the maximum entropy principle, it can be seen as that when this ratio obeys uniform distribution, the synthetic power in every snapshot can increase, which means that the phase error compensation is effective for the signals in every snapshot. Therefore, the objective function should be the maximum entropy of that ratio.

Assuming that there is a reference source in the target scene, if the array beamforming is expected to form a main-lobe at the source, the weighted value of the array can be obtained as *W* = *w*<sup>1</sup> *w*<sup>2</sup> ··· *wN* T . When the snapshots are certain, *Qg*, representing the total ideal synthetic power of all snapshots at the reference point, is a constant, and can be obtained by

$$Q\_{\mathcal{S}} = \sum\_{ss=1}^{SS} \left| G(ss) \right|^2 \tag{17}$$

where, *SS* is the total number of snapshots, and *G*(*ss*) is the synthesized signal at the reference source after ideally weighting and summing the received signals in the *ss*th snapshot. The synthetic power in the main-lobe direction of the array's received signals in one snapshot can be obtained by

$$\log(ss) = \sum\_{n=1}^{N} w\_n^\* \cdot (s(\pounds\_n) + v\_n)\_{ss} \cdot e^{i\rho\_n} \tag{18}$$

where, (*s*(*x*ˆ*n*) + *vn*)*ss* is the received signal in the *ss*th snapshot of the array element, and *ϕ<sup>n</sup>* is the phase error of the *n*th array element to be estimated. Because *Qg* is a constant, the entropy objective function of *g*(*ss*) can be written as

$$E\_{\mathcal{S}} = -\frac{1}{Q\_{\mathcal{S}}} \sum\_{ss=1}^{SS} \left| \mathcal{g}(ss) \right|^2 \ln \frac{\left| \mathcal{g}(ss) \right|^2}{Q\_{\mathcal{S}}} \tag{19}$$

The entropy objective function *Eg* is the function of the phase *ϕ<sup>n</sup>* to be estimated. Thus, the phase estimation based on the maximum entropy can be expressed as

$$\boldsymbol{\phi}\_{\boldsymbol{n}} = \mathbf{argmax}\_{\boldsymbol{\varphi}\_{\boldsymbol{n}}} E\_{\boldsymbol{\mathcal{S}}}(\boldsymbol{\varphi}\_{\boldsymbol{n}}) \tag{20}$$

**Theorem 2.** ([25]) *Suppose that E*(*x*), *which is the* second derivation of *E*, *is continuous in an open neighborhood of* a *, and that E*(a) = 0 *and E*(*x*) *is negative definite. Then* a *is a strict local maximizer of E*(*x*).

**Proof.** Because *E*(*x*), which is also known as the Hessian, is a continuous and negative definite at a, we can choose a radius *r* > 0 so that *E*(*x*) remains negative definite for all *x* in the open ball *D* = {*z*|*z* − a < *r* }. Taking any nonzero vector *p* with *p* ≤ *r*, we have <sup>a</sup> <sup>+</sup> *<sup>p</sup>* <sup>∈</sup> *<sup>D</sup>*, and so *<sup>E</sup>*(<sup>a</sup> <sup>+</sup> *<sup>p</sup>*) <sup>=</sup> *<sup>E</sup>*(a) <sup>+</sup> *<sup>p</sup>*T*E* (a) + <sup>1</sup> <sup>2</sup> *<sup>p</sup>*T*E*(*z*)*<sup>p</sup>* <sup>=</sup> *<sup>E</sup>*(a) <sup>+</sup> <sup>1</sup> <sup>2</sup> *<sup>p</sup>*T*E*(*z*)*p*, where, *<sup>z</sup>* <sup>=</sup> <sup>a</sup> <sup>+</sup> *tp* for some *<sup>t</sup>* <sup>∈</sup> (0, 1). Since *<sup>z</sup>* <sup>∈</sup> *<sup>D</sup>*, we have *<sup>p</sup>*T*E*(*z*)*<sup>p</sup>* <sup>&</sup>lt; 0, and therefore *E*(a + *p*) < *E*(a), giving the results.

By combining Theorem 1 and Theorem 2, we can use the Newton method for Formula (20). There is an iterative solution of *ϕ<sup>n</sup>* as

$$\left. \boldsymbol{\varrho}\_{n}^{(l+1)} \right| = \left. \boldsymbol{\varrho}\_{n}^{(l)} - \left( \frac{\partial^{2} E\_{\boldsymbol{\xi}}}{\partial \boldsymbol{\varrho}\_{n}^{2}} \right)^{-1} \frac{\partial E\_{\boldsymbol{\xi}}}{\partial \boldsymbol{\varrho}\_{n}} \right|\_{\boldsymbol{\varrho}\_{n} = \boldsymbol{\varrho}\_{n}^{(l)}} \tag{21}$$

where, superscript (*l*) indicates the *l*th iteration. To solve Equation (21), the first and second derivation of *Eg* against *ϕ<sup>n</sup>* need to be calculated. The first derivative expression can be expressed as

$$\frac{\partial E\_{\%}}{\partial \varphi\_{n}} = -\frac{1}{Q\_{\%}} \sum\_{ss=1}^{SS} \left( 1 + \ln \frac{\left| \mathcal{g}(ss) \right|^{2}}{Q\_{\%}} \right) \cdot \frac{\left\| \mathcal{g}(ss) \right\|^{2}}{\partial \varphi\_{n}} \tag{22}$$

The second derivative expression can be expressed as

$$\frac{\partial^2 E\_{\mathcal{S}}}{\partial \rho\_n^2} = -\frac{1}{Q\_{\mathcal{S}}} \sum\_{ss=1}^{SS} \left( \frac{1}{|\mathcal{g}(ss)|^2} \cdot \left( \frac{\partial |\mathcal{g}(ss)|^2}{\partial \rho\_n} \right)^2 + \left( 1 + \ln \frac{|\mathcal{g}(ss)|^2}{Q\_{\mathcal{S}}} \right) \cdot \frac{\partial^2 |\mathcal{g}(ss)|^2}{\partial \rho\_n^2} \right) \tag{23}$$

where,

$$\begin{split} \frac{\partial \left[ \mathcal{g}(ss) \right]^2}{\partial \boldsymbol{\varrho}\_n} &= \mathcal{g}(ss) \frac{\partial \mathcal{g}^\*(ss)}{\partial \boldsymbol{\varrho}\_n} + \mathcal{g}^\*(ss) \frac{\partial \mathcal{g}(ss)}{\partial \boldsymbol{\varrho}\_n} \\ &= 2 \text{Re} \left( \mathcal{g}^\*(ss) \frac{\partial \mathcal{g}(ss)}{\partial \boldsymbol{\varrho}\_n} \right) \end{split} \tag{24}$$

$$\frac{\partial g(ss)}{\partial \varphi\_n} = jw\_n^\* \cdot (w\_{en}s(\chi\_n) + v\_{ll})\_{ss} \cdot e^{\mathrm{i}\varphi\_n} \tag{25}$$

$$\begin{split} \left(\frac{\partial \left|\mathcal{g}(s\mathbf{s})\right|^{2}}{\partial \boldsymbol{\varrho}\_{n}}\right)^{2} &= \left(\mathcal{g}(s\mathbf{s})\frac{\partial \boldsymbol{\chi}^{\*}(s\mathbf{s})}{\partial \boldsymbol{\varrho}\_{n}} + \mathcal{g}^{\*}(s\mathbf{s})\frac{\partial \boldsymbol{\chi}(s\mathbf{s})}{\partial \boldsymbol{\varrho}\_{n}}\right)^{2} \\ &= 2\text{Re}\left(\left(\mathcal{g}^{\*}(s\mathbf{s})\frac{\partial \boldsymbol{\chi}(s\mathbf{s})}{\partial \boldsymbol{\varrho}\_{n}}\right)^{2}\right) + 2|\mathcal{g}(s\mathbf{s})|^{2} \left|\frac{\partial \boldsymbol{\chi}(s\mathbf{s})}{\partial \boldsymbol{\varrho}\_{n}}\right|^{2} \end{split} \tag{26}$$

$$\frac{\left|\partial^2 \vert \mathcal{g}(\text{ss})\right|^2}{\partial \rho\_n^2} = 2 \left|\frac{\partial \mathcal{g}(\text{ss})}{\partial \rho\_n}\right|^2 + 2 \text{Re}\left(\mathcal{g}^\*(\text{ss}) \frac{\partial^2 \mathcal{g}(\text{ss})}{\partial \rho\_n^2}\right) \tag{27}$$

$$\frac{\partial^2 g(\text{ss})}{\partial \rho\_n^2} = -w\_n^\* \cdot \left( w\_{\text{en}} s(\text{x}\_n) + v\_n \right)\_{\text{ss}} \cdot e^{i\rho\_n} \tag{28}$$

where, Re(·) represents the real part operation.

In practical cases, the phase errors in signals received by different array elements are relatively independent from each other [26]. Therefore, the phase error of each array element is searched separately, and in Formula (23), only the diagonal elements of the Hessian, which is also the second derivative *<sup>∂</sup>*2*Eg ∂ϕ*<sup>2</sup> *n* , of the entropy with respect to phase error are derived, while the off-diagonal elements *<sup>∂</sup>*2*Eg ∂ϕn∂ϕ<sup>m</sup>* , *<sup>n</sup>* = *<sup>m</sup>* are ignored.

By taking Equations (24) and (25) into Equation (22), the analytical formula of the first derivative is obtained

$$\frac{\partial E\_{\mathcal{S}}}{\partial \rho\_n} = -\frac{2}{Q\_{\mathcal{S}}} \sum\_{ss=1}^{\text{SS}} \left( 1 + \ln \frac{|\mathbf{g}(\mathbf{s}s)|^2}{Q\_{\mathcal{S}}} \right) \cdot \text{Re} \left( \dot{\mathbf{y}} \cdot \mathbf{g}^\*(\mathbf{s}s) \cdot w\_n^\* \cdot \left( w\_{\sigma n} \mathbf{s}(\mathbf{x}\_n) + v\_n \right)\_{\text{ss}} \cdot \mathbf{e}^{i\rho\_n} \right) \tag{29}$$

By taking Equation (24) to Equation (28) into Equation (23), the analytical formula of the first derivative is obtained

$$\frac{\partial^{2}E\_{\varepsilon}}{\partial\varphi\_{n}^{2}} = -\frac{2}{Q\_{\varepsilon}}\sum\_{n=1}^{\infty} \left| \frac{\text{Re}\left( \left( \text{jg}^{\circ}\left(\text{ss}\right)\text{w}\_{n}^{\circ}\left(\text{w}\_{\text{as}}\text{s}\left(\text{x}\_{\text{s}}\right)+\text{v}\_{\text{u}}\right)\_{\text{ue}}e^{i\theta\_{n}} \right)^{2} \right) + \left| \text{g}\left(\text{ss}\right) \right|^{2} \left| \text{jw}\_{n}^{\circ}\left(\text{w}\_{\text{as}}\text{s}\left(\text{x}\_{\text{s}}\right)+\text{v}\_{\text{u}}\right)\_{\text{ue}}e^{i\theta\_{n}} \right|^{2} \right) \right|$$

$$+ \left( 1 + \ln\frac{\left| \text{g}\left(\text{ss}\right) \right|^{2}}{Q\_{\varepsilon}} \right) \left( \left| \text{jw}\_{\text{u}}^{\circ}\left(\text{w}\_{\text{as}}\text{s}\left(\text{x}\_{\text{s}}\right)+\text{v}\_{\text{u}}\right)\_{\text{ue}}e^{i\theta\_{\varepsilon}} \right|^{2} - \text{Re}\left( \text{g}^{\circ}\left(\text{ss}\right)\text{w}\_{\text{u}}^{\circ}\left(\text{w}\_{\text{as}}\text{s}\left(\text{x}\_{\text{s}}\right)+\text{v}\_{\text{u}}\right)\_{\text{ue}}e^{i\theta\_{\varepsilon}} \right) \right| \right) \tag{30}$$

After performing phase error calibration for Equation (16), the beamforming with error calibration is

$$\tilde{\mathbf{Y}} = \mathbf{W}\_c^\mathrm{H} \odot \mathbf{W}^\mathrm{H} \mathbf{X} = \sum\_{n=1}^{N} w\_{cn}^\* \cdot w\_n^\* \cdot (w\_{cn} s(\mathbf{x}\_n) + v\_n) \tag{31}$$

ρ

ϕ=

˖

˖

where, represents the Hadamard product of the matrix, and *<sup>W</sup><sup>c</sup>* <sup>=</sup> *wc*<sup>1</sup> *wc*<sup>2</sup> ... *wcN*<sup>T</sup> = *e*−*jϕ*ˆ1 *e*−*jϕ*ˆ2 ... *e*−*jϕ*ˆ*<sup>N</sup>* T .

#### *3.2. Initial Value Estimation of Phase Error Based on PG*

According to the principle of Newton's maximum entropy algorithm, the selection of the initial value of the phase error will directly affect the iteration efficiency. If the initial value of the phase error is randomly selected, such as directly setting the initial value of phase error to zero, which is the N-ME algorithm, it may not only greatly slow down the convergence speed of the algorithm, but also make the algorithm fall into the local optimal solution. Therefore, a method based on a reference source to select the initial value of the phase error is proposed to improve the N-ME algorithm.

When there are noises and array position errors, after weighting the received signals, the initial received signals of the array without phase error calibration can be written as

$$\begin{aligned} \mathbf{S} &= \mathbf{W}^{\prime} \odot \hat{\mathbf{X}} = \begin{bmatrix} S\_1 & S\_2 & \dots & S\_N \end{bmatrix}^T \\ &= \begin{bmatrix} \check{A}\text{w}\_1^{\prime} e^{-\bar{\boldsymbol{\beta}}\left(-\sin\theta \,\boldsymbol{z}\_1 + \frac{\cos^2\theta}{2R} \boldsymbol{z}\_1^{\prime}\right)} \cdot e^{\bar{\boldsymbol{\beta}}\boldsymbol{\theta}^{(0)}} & \check{A}\text{w}\_2^{\prime} e^{\bar{\boldsymbol{\beta}}\left(-\sin\theta \,\boldsymbol{z}\_2 + \frac{\cos^2\theta}{2R} \boldsymbol{z}\_2^{\prime}\right)} \cdot \boldsymbol{e}^{\bar{\boldsymbol{\beta}}\boldsymbol{\theta}^{(0)}} & \cdots & \check{A}\text{w}\_N^{\prime} e^{\bar{\boldsymbol{\beta}}\left(-\sin\theta \,\boldsymbol{z}\_N + \frac{\cos^2\theta}{2R} \boldsymbol{z}\_N^{\prime}\right)} \cdot e^{\bar{\boldsymbol{\beta}}\boldsymbol{\theta}^{(0)}} \end{bmatrix}^{\prime} \end{aligned} \tag{32}$$

where, [·] <sup>∗</sup> represents conjugation, *Sn*, *n* = 1, ··· , *N* represents the initial output signal of the *n*th array element, and *<sup>ϕ</sup>*(0) *<sup>n</sup>* , *n* = 1, ··· , *N* represents the initial phase error of the *n*th array element caused by noise and array position error.

Ideally, after weighting the received signals, the phase error caused by different distances is fully compensated, and the PG between the output signals of the array elements is constant. For discrete sequences, the phase difference between the output signals of adjacent array elements is the PG. We define the correlation sequence of the array as

$$\boldsymbol{\rho} = \begin{bmatrix} \rho\_1 & \rho\_2 & \cdots & \rho\_{N-1} \end{bmatrix}^T = \begin{bmatrix} \mathbb{S}\_1 \cdot \mathbb{S}\_2 & \mathbb{S}\_2 \cdot \mathbb{S}\_3 & \cdots & \mathbb{S}\_{N-1} \cdot \mathbb{S}\_N \end{bmatrix}^T \tag{33}$$

The PG between adjacent elements is estimated as

$$
\Delta \mathfrak{g}^{(0)} = \begin{bmatrix} \Delta \mathfrak{q}\_1^{(0)} & \Delta \mathfrak{q}\_2^{(0)} & \cdots & \Delta \mathfrak{q}\_{N-1}^{(0)} \end{bmatrix}^\mathrm{T} = \mathrm{arg}(\mathfrak{p})\tag{34}$$

By taking the phase of the first array element as a reference, the initial value of phase error to be compensated for each array element is estimated by calculating the cumulative sum between adjacent array elements. There is

$$\boldsymbol{\mathfrak{q}}^{(0)} = \begin{bmatrix} \boldsymbol{\varrho}\_1^{(0)} & \boldsymbol{\varrho}\_2^{(0)} & \cdots & \boldsymbol{\varrho}\_N^{(0)} \end{bmatrix}^T = \begin{bmatrix} 0 & \text{cusum} \left(\boldsymbol{\Delta}\boldsymbol{\mathfrak{q}}^{(0)}\right) \end{bmatrix}^T \tag{35}$$

 ϕ

> =

> > ϕ

 =

 =

=where, cusum(·) represents the cumulative sum.

The IN-ME phase error compensation algorithm, based on the reference source, is shown in Algorithm 1.

 =

ρ

=

ϕ=

**Algorithm 1** The IN-ME phase error compensation algorithm based on reference source.

**Input:** received signals of array *X*ˆ , reference source *P*(*xP*, *yP*) **Output:** beamform after phase error calibration *<sup>Y</sup>*<sup>7</sup>

Step 1: Calculate weighted signal, weights *W* = *w*<sup>1</sup> *w*<sup>2</sup> ... *wN* T , and weighted signal *<sup>S</sup>* <sup>=</sup> *<sup>W</sup>*<sup>∗</sup> *<sup>X</sup>*<sup>ˆ</sup> ;

Step 2: Calculate the correlation sequence of the array, *ρ* = *S*<sup>1</sup> · *S*<sup>2</sup> *S*<sup>2</sup> · *S*<sup>3</sup> ... *SN*−<sup>1</sup> · *SN* T ; Step 3: Estimate initial value of phase error based on PG: 1. Estimate the PG between adjacent elements,

$$\Delta \boldsymbol{\varrho}^{(0)} = \begin{bmatrix} \Delta \boldsymbol{\varrho}\_1^{(0)} & \Delta \boldsymbol{\varrho}\_2^{(0)} & \cdots & \Delta \boldsymbol{\varrho}\_{N-1}^{(0)} \end{bmatrix}^T = \arg(\boldsymbol{\rho});$$

$$\begin{aligned} &2. \text{ Estimate initial value of phase error,} \\ &\varphi^{(0)} = \begin{bmatrix} \varphi\_1^{(0)} & \varphi\_2^{(0)} & \dots & \varphi\_N^{(0)} \end{bmatrix}^T = \begin{bmatrix} 0 & \text{cusum}\left(\Delta\varphi^{(0)}\right) \end{bmatrix}^T; \end{aligned}$$

Step 4: Improved Newton Maximum Entropy Algorithm:

1. Calculate the array's synthetic signal at the reference source,

$$\mathcal{g}(\mathbf{s}s) = \sum\_{n=1}^{N} \mathbb{H}\_{\mathrm{H}}^{\*} \cdot \left(\mathbf{s}(\mathbb{X}\_{\mathrm{H}}) + \boldsymbol{\nu}\_{\mathrm{H}}\right)\_{\mathrm{ss}} \cdot \mathbf{c}^{\|\boldsymbol{\varphi}\_{\mathrm{g}}\|\_{\mathrm{g}}}$$

2. Form the entropy objective function, *Eg* <sup>=</sup> <sup>−</sup> <sup>1</sup> *Qg SS* ∑ *ss*=1 |*g*(*ss*)| <sup>2</sup> ln <sup>|</sup>*g*(*ss*)<sup>|</sup> 2 *Qg* ; 3. Calculate the first and second derivatives of *Eg* to *ϕn*,

*∂Eg ∂ϕ<sup>n</sup>* <sup>=</sup> <sup>−</sup> <sup>2</sup> *Qg SS* ∑ *ss*=1 1 + ln <sup>|</sup>*g*(*ss*)<sup>|</sup> 2 *Qg* · Re *j* · *g*∗(*ss*) · *w*<sup>∗</sup> *<sup>n</sup>* · (*wens*(*xn*) <sup>+</sup> *vn*)*ss* · *<sup>e</sup>*j*ϕ<sup>n</sup>* , *∂*<sup>2</sup>*Eg ∂ϕ*<sup>2</sup> *n* <sup>=</sup> <sup>−</sup> <sup>2</sup> *Qg SS* ∑ *ss*=1 ⎛ ⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎝ 1 |*g*(*ss*)| 2 · ⎛ ⎜⎜⎜⎝ Re *jg*∗(*ss*)*w*∗ *<sup>n</sup>*(*wens*(*xn*) + *vn*)*ssejϕ<sup>n</sup>* 2 + |*g*(*ss*)| 2 *jw*<sup>∗</sup> *<sup>n</sup>*(*wens*(*xn*) + *vn*)*ssejϕ<sup>n</sup>* 2 ⎞ ⎟⎟⎟⎠ + 1 + ln <sup>|</sup>*g*(*ss*)<sup>|</sup> 2 *Qg* · ⎛ ⎜⎜⎝ *jw*<sup>∗</sup> *<sup>n</sup>*(*wens*(*xn*) + *vn*)*ssejϕ<sup>n</sup>* 2 − Re *g*∗(*ss*)*w*∗ *<sup>n</sup>*(*wens*(*xn*) + *vn*)*ssejϕ<sup>n</sup>* ⎞ ⎟⎟⎠ ⎞ ⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎠ ; *<sup>∂</sup>*<sup>2</sup>*Eg* <sup>−</sup><sup>1</sup> *<sup>∂</sup>Eg* ;

4. Use Newton iterative method for optimization, *<sup>ϕ</sup>*(*l*+1) *<sup>n</sup>* <sup>=</sup> *<sup>ϕ</sup>*(*l*) *n* − *∂ϕ*<sup>2</sup> *n ∂ϕ<sup>n</sup> <sup>ϕ</sup>n*=*ϕ*(*l*) *n*

Step 5: Beamform after error calibration, *<sup>Y</sup>*<sup>7</sup> <sup>=</sup> *<sup>W</sup>*<sup>H</sup> *<sup>c</sup> <sup>W</sup>*H*X*<sup>ˆ</sup> <sup>=</sup> *<sup>N</sup>* ∑ *n*=1 *w*∗ *cn* · *w*<sup>∗</sup> *<sup>n</sup>* · (*wens*(*xn*) + *vn*).

#### **4. Effectively Compensated Area**

The phase error compensation based on phase error compensation function of the IN-ME algorithm is aimed at the situation in which the expected main-lobe of beamforming is just at the reference source *P*(*xP*, *yP*). When the expected main-lobe of beamforming is not at the reference source, there is a mismatch between the phase error compensation function and the actual phase error. Thus, it is necessary to discuss the effectively compensated area based on the phase error compensation function of the IN-ME algorithm.

Assuming that the position errors of the array elements are independent of each other and meet the normal distribution with the mean value of 0 and the variance of *σ*2, which is, for the *n*th UAV element, Δ*xn* ∼ *N* 0, *σ*<sup>2</sup> and Δ*yn* ∼ *N* 0, *σ*<sup>2</sup> , and the two-dimensional joint distribution of the two also meets the normal distribution, according to the "3*σ* rule", the elements can be regarded approximately as that they are located in a circle with the radius centered on the ideal array element position, as shown in Figure 2. When the expected main-lobe of beamforming is at the source *B*(*xB*, *yB*), but the reference source of the IN-ME algorithm is *P*(*xP*, *yP*), the effectively compensated area based on the phase error compensation function of the IN-ME algorithm can be approximately solved by geometric methods.

**Figure 2.** Geometric relation between the array elements, with position error and radiant source.

Obviously, when the actual position of the array element is on the circle, the phase error is larger than when it is in the circle. Point *A* is set as a moving point on the circle to represent the actual position of the *n*th array element. When the two sources, *P*(*xP*, *yP*) and *<sup>B</sup>*(*xB*, *yB*), have relative phase errors <sup>|</sup>Δ*ϕBP*<sup>|</sup> <sup>≤</sup> <sup>π</sup> <sup>4</sup> , it is considered that the compensation is effective [27,28]. There is

$$\begin{cases} \max fit = |(PA - PO) - (BA - BO)| \le \frac{\lambda}{8} \\ \text{s.t.} \\ \mathbf{x}\_A^2 + \mathbf{y}\_A^2 = r\_\varepsilon^2 \end{cases} \tag{36}$$

where, *xA* and *yA* are respectively the *x* axis and *y* axis coordinates of point *A.*

As shown in Figure 2, crossing point *A*, the plumb line *AF*, which is perpendicular to *PO*, is made, and the point of intersection is point *F*. The plumb line *AE*, which is perpendicular to *BO*, is made, and the point of intersection is point *E*. We let ∠*AOB* = *αB*, ∠*AOP* = *αP*, and ∠*BOP* = Δ*α*. When *PO re* and *BO re*, the electromagnetic wave from *P*(*xP*, *yP*) or *B*(*xB*, *yB*) to any point on the circle can be approximately simplified to the plane wave, and we have |*PA* − *PO*| ≈ *FO* and |*BA* − *BO*| ≈ *EO*. Then Equation (36) is approximately simplified as

$$\begin{cases} \max fit = |r\_\varepsilon \cos(\alpha\_P) - r\_\varepsilon \cos(\alpha\_B)| \le \frac{\lambda}{8} \\ \text{s.t.} \\ PO \gg r\_\varepsilon \\ BO \gg r\_\varepsilon \\ \alpha\_P = \alpha\_B + \Delta \alpha \end{cases} \tag{37}$$

By using sum-to-product addition formulas, there is *fit* <sup>=</sup> <sup>2</sup>*re* sin <sup>Δ</sup>*<sup>α</sup>* 2 · sin *α<sup>B</sup>* + <sup>Δ</sup>*<sup>α</sup>* 2 . Then there is

$$\max(fit) = 2r\_t \sin\left(\frac{\Delta a}{2}\right), \text{ when } a\_B + \frac{\Delta a}{2} = \frac{\pi}{2} + i\pi, \text{ } i = \pm 1, \pm 2, \dots \tag{38}$$

As shown in Figure 2, the perpendicular line of the Δ*α* angular bisector which is made through the center *O* of the circle intersects with the circle at A and C . Equation (38) shows that the expected point *A* which can get max(*fit*) approximately locates at A or C , according to

$$\max(fit) = 2r\_{\varepsilon}\sin\left(\frac{\Delta\alpha}{2}\right) \le \frac{\lambda}{8} \tag{39}$$

The area of *B*(*xB*, *yB*) that can be compensated based on *P*(*xP*, *yP*) for the *n*th array element can be obtained by

$$
\Delta a\_{\text{max}} = 2 \arcsin(\lambda / (16r\_{\varepsilon})) \tag{40}
$$

After finding out the upper and lower bounds of the area of *B*(*xB*, *yB*) the phase error for each array element can be effectively compensated, and the area can be divided by finding the intersection of the area of all array elements. As shown in Figure 3, the connecting line from the *n*th array element to *P*(*xP*, *yP*) is marked as line *an*, *n* = 1, 2, 3, ... , *N*, which is the angular bisector of 2Δ*α*, and the line towards the negative angle direction is line *bn*, while the line towards the positive angle direction is line *cn*. The intersection point of line *c*<sup>1</sup> and line *bN* is marked as point *D*. The intersection point of line *c*<sup>1</sup> and line *c*<sup>2</sup> is marked as point *G*. The line *g* parallel to line *cN* is started from point *G*, and the connecting line between *G* and *D* is marked as line *e*. The intersection point of *bN*−<sup>1</sup> and *bN* is marked as point *H*, and line *f* parallel to *b*<sup>1</sup> is started from *H*. The connecting line between *H* and *D* is marked as line *d*. Then the effectively compensated area of the array for *B*(*xB*, *yB*) based on *P*(*xP*, *yP*) is within the area encircled by four sides of *g*, *e*, *d*, and *f*.

**Figure 3.** The effectively compensated area of emitter *B*(*xB*, *yB*).

As shown in Figure 3, the problem of finding the effectively compensated area is transformed into the problem of finding the intersection point of the lines, and the equations of the black line cluster *an*, *n* = 1, 2, 3, . . . , *N* can be expressed as

$$y = \frac{y\_P}{\mathbf{x}\_P - \mathbf{x}\_n} (\mathbf{x} - \mathbf{x}\_n) \tag{41}$$

The equation of the green line cluster *bn*, *n* = 1, 2, 3, . . . , *N* can be expressed as

$$y = \tan\left(\arctan\left(\frac{y\_P}{\mathbf{x}\_P - \mathbf{x}\_\mathrm{n}}\right) + \Delta a\right)(\mathbf{x} - \mathbf{x}\_\mathrm{n})\tag{42}$$

The equation of the blue line cluster *cn*, *n* = 1, 2, 3, . . . , *N* can be expressed as

$$y = \tan\left(\arctan\left(\frac{y\_P}{\mathbf{x}\_P - \mathbf{x}\_n}\right) - \Delta\alpha\right)(\mathbf{x} - \mathbf{x}\_\mathbb{N})\tag{43}$$

Let *kan* = *yP xP*−*xn* , *kbn* <sup>=</sup> tan arctan *yP xP*−*xn* + Δ*α* , and *kcn* <sup>=</sup> tan arctan *yP xP*−*xn* − Δ*α* , and coordinates of *H*, *D* and *G* can be obtained by

$$\begin{cases} \mathbf{x}\_{H} = \frac{k\_{b(N-1)}\mathbf{x}\_{(N-1)} - k\_{bN}\mathbf{x}\_{N}}{k\_{b(N-1)} - k\_{bN}}, y\_{H} = k\_{b(N-1)}k\_{bN}\frac{\mathbf{x}\_{(N-1)} - \mathbf{x}\_{N}}{k\_{b(N-1)} - k\_{bN}} \\\\ \mathbf{x}\_{D} = \frac{k\_{c1}\mathbf{x}\_{1} - k\_{bN}\mathbf{x}\_{N}}{k\_{c1} - k\_{bN}}, y\_{D} = k\_{c1}k\_{bN}\frac{\mathbf{x}\_{1} - \mathbf{x}\_{N}}{k\_{c1} - k\_{bN}} \\\\ \mathbf{x}\_{G} = \frac{k\_{c1}\mathbf{x}\_{1} - k\_{c2}\mathbf{x}\_{2}}{k\_{c1} - k\_{c2}}, y\_{G} = k\_{c1}k\_{c2}\frac{\mathbf{x}\_{1} - \mathbf{x}\_{2}}{k\_{c1} - k\_{c2}} \end{cases} \tag{44}$$

The analytic expressions of the red four sides *g*, *e*, *d*, and *f* are

$$\begin{cases} f: y = k\_{b1}(\mathbf{x} - \mathbf{x}\_{H}) + y\_{H} \\ d: y = k\_{bN}(\mathbf{x} - \mathbf{x}\_{D}) + y\_{D} \\ e: y = k\_{c1}(\mathbf{x} - \mathbf{x}\_{D}) + y\_{D} \\ g: y = k\_{cN}(\mathbf{x} - \mathbf{x}\_{G}) + y\_{G} \end{cases} \tag{45}$$

The area encircled by the four edges shown in Equation (45) is the effectively compensated area of *B*(*xB*, *yB*) that can be compensated based on IN-ME algorithm.

By estimating and compensating the phase errors, the phase errors of the received signals caused by the position of the array elements in the array are calibrated so that the antenna array of UAVs can effectively carry out beamforming. The analysis of the effectively compensated area based on the reference source can provide a reference for the target area division and path planning of UAVs.

#### **5. Simulation Analysis**

#### *5.1. IN-ME Phase Error Compensation Algorithm Based on Reference Source*

We suppose that the signal frequency *f*<sup>0</sup> = 300 MHz, wavelength *˘* = 1 m, number of array elements *N* = 26, array aperture *L* = 500*˘*, main-lobe of beam is at known radiant source *P*(*xP*, *yP*), whose polar coordinates are *θ<sup>P</sup>* = 0 ◦ or *θ<sup>P</sup>* = 6 ◦ and *RP* = 30 km, the number of signal snapshots *SS* = 100, array element position error Δ*xn* and Δ*yn* are independent of each other, and they both meet the normal distribution with a mean value of 0 and a standard deviation of 3*˘*, and the low SNR is −3 dB. We use the N-ME algorithm, IN-ME algorithm, adaptive differential evolution algorithm (ADE) [29], and genetic algorithm (GA), respectively, to estimate the phase error; the changes in the maximum entropy of the array signal synthetic power, along with the number of iterations of the algorithms, are shown in Figure 4. The reason why we adopt ADE and GA is because they are both artificial intelligence algorithms, which are different from our algorithm and are widely used in solving multi-objective search problems and beamforming [29–33]. The objective function of both ADE and GA is the entropy objective function, as shown in Formula (19). The parameters of ADE are that the population number of the genetic algorithm is 100, the mutation rate is 0.5, and the crossover probability is 0.9. The parameters of GA are that the hybridization rate is 0.7, the mutation rate is 0.05, and the "roulette wheel" selection is used. The elitist retention strategy is adopted. It is worth noting that all of our simulation tests are respectively based on the average of 1000 Monte Carlo tests.

As shown in Figure 4, the convergence rates of ADE and GA are slower than those of N-ME and IN-ME, and within 50 iterations, the entropy value of ADE and GA cannot reach the same high rate as those reached by N-ME or IN-ME, which shows that ADE and GA are inferior to N-ME and IN-ME, and this is because the search direction of ADE and GA is aimless, while N-ME and IN-ME search along the gradient descent direction. Compared with IN-ME, the N-ME algorithm will require more iterations and a longer period to converge, and will easily fall into local optimization, eventually reaching a smaller entropy value. Compared with the N-ME algorithm, the IN-ME algorithm proposed in

this paper makes the entropy objective function of the algorithm converge faster and reach a larger entropy value, which verifies the effectiveness of IN-ME. Because of the adverse performance of ADE and GA, N-ME and IN-ME are used to form beams and make comparisons. After phase error compensation, when *θ<sup>P</sup>* = 0 ◦ , beamforming is shown in Figure 5a. When the beam is scanning and *θ<sup>P</sup>* = 6 ◦ , the compensation effect is shown in Figure 5b.

**Figure 4.** Comparison of different algorithms. (**a**) *θ<sup>P</sup>* = 0 ◦ ; (**b**) *θ<sup>P</sup>* = 6 ◦ .

**Figure 5.** Comparison of the beam forming effect after phase compensation. (**a**) *θ<sup>P</sup>* = 0 ◦ ; (**b**) *θ<sup>P</sup>* = 6 ◦ .

As shown in Figure 5, whether beam scans or not, beams cannot be formed because of the uncompensated phase error; the sidelobe level of the N-ME algorithm is approximately −9.5 dB, and the sidelobe level of IN-ME algorithm is approximately −12 dB, which means that the beamforming algorithm proposed in this paper can better compensate for the phase error of array elements and obtain a lower sidelobe level.

In order to compare the ideal beamforming (IB) with noise, phase gradient beamforming (PGB), the N-ME beamforming algorithm, and the IN-ME beamforming algorithm in different SNR backgrounds, by taking beamforming without scanning as an example, scenarios in which SNR changes from −12 dB to 9 dB are simulated, and the beamforming results are shown in Figure 6.

As shown in Figure 6, with the SNR changing from −12 dB to 9 dB, both the N-ME and IN-ME algorithms can effectively form a beam pattern. When the SNR is low, using the phase error estimated by PG information to compensate the array cannot effectively form a beam pattern, but with increasing SNR, the beam formed by PG gradually improves, compared to the beam formed by the N-ME. The IN-ME beamforming algorithm proposed in this paper can continuously remain nearly the same as ideal beamforming with noise; when the SNR is low, it performs better than PGB, and when SNR is high, it performs better than the N-ME beamforming algorithm, which shows the advantages of the IN-ME algorithm.

**Figure 6.** Beamforming contrast of two algorithms when SNR changes. (**a**) SNR = −12 dB; (**b**) SNR = −9 dB; (**c**) SNR = −6 dB; (**d**) SNR = −3 dB; (**e**) SNR = 0 dB; (**f**) SNR = 3 dB; (**g**) SNR = 6 dB; (**h**) SNR = 9 dB.

#### *5.2. Effectively Compensated Area*

We assume that the position errors corresponding to the phase errors Δ*xr* and Δ*yr* are independent of each other and meet the normal distribution of 0 mean value and *σ<sup>r</sup>* = *<sup>λ</sup>* <sup>2</sup> , which is Δ*xr* ∼ *N* 0, *<sup>λ</sup>* 2 2 and Δ*yr* ∼ *N* 0, *<sup>λ</sup>* 2 2 . According to the "3*σ* rule", we can get *re* <sup>=</sup> 1.5*λ*, and <sup>Δ</sup>*<sup>α</sup>* <sup>≈</sup> 4.77◦ from Equation (40). The area that can be effectively compensated by each array element is within the area whose angular bisector is the line from the known radiant source *P*(*xP*, *yP*) shown as red \* to the ideal position of the array element, and the angular area is between 2Δ*α* = 9.55◦ , as shown in Figure 3. As shown in Figure 2, we set *θ<sup>P</sup>* = 0 ◦ or *θ<sup>P</sup>* = 6 ◦ , *OP* <sup>≡</sup> <sup>30</sup> km, and <sup>∠</sup>*POB* <sup>≡</sup> <sup>Δ</sup>*<sup>α</sup>* <sup>=</sup> 4.77◦ or <sup>∠</sup>*POB* ≡ −Δ*<sup>α</sup>* <sup>=</sup> <sup>−</sup>4.77◦ . According to Equation (36) through Equation (45), the effectively compensated area and the corresponding dividing points *D*, *H*, and *G* are shown as red \* in Figure 7.

**Figure 7.** Effectively compensated area. (**a**) *θ<sup>P</sup>* = 0 ◦ ; (**b**) *θ<sup>P</sup>* = 6 ◦ .

According to Figure 7, we see that when the radiant source is located on the vertical line in the array, the effectively compensated area is greater than when it is located at a certain angle relative to the vertical line of the array. In order to verify the effectively compensated area shown in Figure 7, taking Figure 7a as an example, let the radiant source *M* be a moving source, whose coordinates are *yM* ≡ *yH* and *xM* = −5 km : 0.01 km : 5 km, the main-lobe is expected to point to the direction of the radiant source *M*, and the phase calibration function based on the IN-ME algorithm is used to compensate the beamforming. The movement track of the radiant source *M* is shown by the thick solid green line in Figure 8a. The power value change in the direction of the radiant source *M* is calculated, which is in the direction of the main-lobe, and the boundary points *H*, *G,* and *xP* of *P*(*xP*, *yP*) are marked, as shown in Figure 8b.

**Figure 8.** Change of main-lobe power within the effectively compensated area; (**a**) moving radiant source *M* trajectory; (**b**) variation of main-lobe power value with *xM*.

As shown in Figure 8, due to the influence of noise, when *yM* ≡ *yH*, the power value in the main-lobe direction, where the radiant source *M* is located, fluctuates, but changes regularly as a whole with the change of *xM*. When *xM* ≤ *xH* ∪ *xM* ≥ *xG*, the radiant source *M* is outside the effectively compensated area, and the power value in the main-lobe direction decreases by more than 3 dB, which shows that the compensation is invalid. When *xH* ≤ *xM* ≤ *xG*, the radiant source *M* is within the effectively compensated area, and the power value in the main-lobe direction drops within 3 dB, which shows that the compensation is effective. The closer *xM* is to *xP*, the greater the power value is in the main-lobe direction, which means the better the compensation is.

Similarly, the main-lobe power change within the effectively compensated area during beam scanning is simulated in Figure 9.

**Figure 9.** Change of main-lobe power within the effectively compensated area using beam scanning; (**a**) moving radiant source *M* trajectory; (**b**) variation of main-lobe power value with *xM*.

A conclusion similar to that in Figure 8 is obtained from Figure 9. When the radiant source *M* is outside the effectively compensated area, the power value in the main-lobe direction decreases by more than 1 dB, which is considered that the compensation is noneffective. When the radiant source *M* is within the effectively compensated area and the power value in the main-lobe direction drops within 1.5 dB, which is considered that the compensation is effective. The closer *xM* is to *xP*, the greater the power value is in the main-lobe direction, which means the better the compensation is.

All in all, as can be seen from Figures 8b and 9b, the method proposed in this paper divides the area into a subset of the area with the power value in the main-lobe direction decreasing by 3 dB as a boundary value, ensuring the effectiveness of the division area for compensation. This effectively compensated area is related to the selection of the effective standard for phase compensation in Equation (36).

#### **6. Conclusions**

In order to solve the problems of UAV near-field beamforming, UAV position error compensation, and the effectively compensated area of compensation, this paper firstly constructs a near-field array beamforming model based on element position error, and approximately simplifies the model using the Taylor expansion of the near-field signal phase difference function. Moreover, an IN-ME algorithm is proposed to estimate and compensate for the phase error. Based on the known reference source signal, the initial value of the phase error is estimated through the PG information of the array's received signal, and then the array signal entropy objective function is established to iteratively optimize the phase error value using the Newton iteration method, and the beam pattern is formed after the phase error compensation. Thirdly, when the beam scanning leads to the phase compensation function mismatch, based on the phase error compensation function of the IN-ME algorithm, the effectively compensated area is further divided, which lays a foundation for the subsequent area region division and UAV path planning. Finally, the validity of the conclusion is verified by simulation. By using the IN-ME algorithm proposed in this paper, the influence of a single UAV's position error on array beamforming is effectively suppressed, and the effectively compensated area is divided because of phase compensation function mismatch, which provides a theoretical basis for the UAV beamforming application in wireless communication, electronic reconnaissance, and jamming. As the signal wavelength becomes shorter, the adverse influence of the UAV's position error increases, and the beam pattern will deteriorate sharply, which is worthy of further research.

**Author Contributions:** Conceptualization, S.P., G.W., Y.Z., G.Y. and Y.L.; data curation, Y.Z. and G.W.; formal analysis, Y.Z. and G.W.; funding support, Y.L.; writing—original draft, Y.Z.; writing—review and editing, S.P., G.W., G.Y. and Y.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Military Scientific and Technological Commission, and the Airforce of the People's Liberation Army, Government of China, under grant number 21KJ3C1-0090R and PGGC-2021-003.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

#### **References**

