A Study on a Radio Source Location Estimation System Using High Altitude Platform Stations (HAPS)

1. Introduction

In recent years, mobile communication technology has made great progress, and is used not only in public places but also in our everyday life. For example, mobile communication technology is indispensable in the field of IoT, which is used in smart phones and electrical appliances that everyone has today, and in the field of automatic driving. In fact, the penetration rate of mobile devices is increasing every year [1]. The number of wireless stations has been increasing along with this growth, and as of the end of September 2022, a total of 298.27 million wireless stations were in operation nationwide in Japan as shown in Figure 1 [2].

As time progresses, technologies such as cell phones and IoT have been developed, and mobile communication devices have become more familiar to us, the general public. However, perhaps due to these developments, the number of illegal radio stations, which use power in excess of the limit or illegal frequencies without a license, has been increasing in recent years as shown in Figure 2 [3].

Today, these illegal radio stations are searched in various ways, and one of them is called a radio wave monitoring system. This system monitors radio waves received at sensor stations and estimates the direction of the radio wave source, based on which the location of illegal radio stations is identified by remotely operating sensor stations installed at various locations and sensor stations mounted on vehicles from a center station installed at each general communications station [3]. However, in urban areas where many buildings are densely built or where there is a lot of vegetation, these systems may not work well because there may be no line of sight from the sensor station to the illegal radio station or many reflected waves may be generated. In recent years, research has been conducted on using unmanned aerial vehicles (UAVs) such as drones to gain some altitude in order to improve the location estimation accuracy of the source of radio waves by making it easier to see the source of radio transmission [4]. However, the service area of this system is limited by the UAV’s heights and actual implementation has not yet been achieved.

To address this issue, this paper focuses on the use of the High Altitude Platform System (HAPS) which is a generic term for a system that provides communication services over a wide area by operating an unmanned aircraft, such as an aircraft flown in the stratosphere, like a communication base station [5]. As explained here, HAPS mainly operates unmanned aircraft as base stations, and papers such as "Research and Development of HAPS Mobile Communication System for Rapid Disaster Recovery" [6] have been published HAPS has the advantage of being able to deploy unmanned base stations regardless of ground traffic conditions or environments such as the sea. So far, research has been conducted with the idea of operating unmanned aircraft as the main base station. However, if HAPS is operated as an illegal radio surveillance system, it will be possible to monitor a wider area with less influence from reflected waves and at a lower cost than conventional surveillance systems. If HAPS can also be operated as a base station in times of disaster, it will be effective in estimating the location of mobile terminals in times of disaster.

This paper examines the operation of HAPS as an illegal radio monitoring system and a location estimation system for mobile terminals during disasters and demonstrates its effectiveness mainly through numerical simulations using MATLAB [7]. This paper aims to evaluate the performance of a location estimation system using HAPS for illegal wireless stations and mobile terminals during disasters that occur in various environments and distances through numerical simulations.

2. Methods for Estimating the Location of Radio Transmission Sources

There are two types of systems for estimating the location of radio wave sources such as illegal radio waves and mobile terminals: geometric methods that use various numerical data observed from actual radio wave sources to estimate the location, and statistical methods that estimate the location by comparing the actual observation data with a database obtained using a virtual transmitter in advance. This section describes the geometrical methods currently used to estimate the location of a radio transmitter. In Japan, the DEURAS system [4] is currently used to monitor illegal radio stations. The DEURAS system uses phase differences of the received signals at different antennas to estimate the direction of incoming radio waves. This section briefly explains the methods based on estimation of angle of arrival.

2.1. Estimation Method Using Phase Difference of Received Signals

This section revisits a well-known method for estimating the location of a radio wave source by determining the direction of arrival of a radio wave using the phase difference between received signals at multiple elements based on the same principle as that of interferometry. For simplicity, we will explain the estimation of the direction of arrival when a plane wave is incident on two array antennas as shown in Figure 3. When considering such a situation, the phase difference $∆\varphi$ between the received signals of arrays 1 and 2 can be expressed as in Equation (1).

$$\mathsf{\Delta }\varphi ={\displaystyle \frac{2\pi }{\lambda }}{d}_{\mathrm{a}\mathrm{r}\mathrm{r}\mathrm{a}\mathrm{y}}\sin \theta$$

(1)

If the complex received signal series of each antenna are ${\mathit{x}}_1$ and ${\mathit{x}}_2$, respectively, the phase difference $∆\varphi$ can be expressed as in Equation (2).

$$∆\varphi ={\tan }^{-1}{\displaystyle \frac{\mathrm{I}\mathrm{m}\left\{E\right[{\mathit{x}}_1{\mathit{x}}_2^{*}\left]\right\}}{\mathrm{R}\mathrm{e}\left\{E{[\mathit{x}}_1{\mathit{x}}_2^{*}\right]\}}}$$

(2)

where E[∙] denotes an ensemble average operation, $\lambda$ is the wavelength of the carrier wave, and $d_{\mathrm{a}\mathrm{r}\mathrm{r}\mathrm{a}\mathrm{y}}$ is the antenna spacing. When the antenna spacing is 1/2 of the carrier wave wavelength, $\theta = {\sin }^{-1}\left({\displaystyle \frac{∆\varphi }{\pi }}\right)$ from Equation (1), and the direction of arrival $\theta$ can be obtained in the range $-{\displaystyle \frac{\pi }{2}}<\theta <{\displaystyle \frac{\pi }{2}}$ by determining the phase difference $∆\varphi$. In a previous study [8], in order to improve the estimation accuracy, the direction of arrival is obtained by using several pairs of two-element arrays to obtain in advance the complex reception response $a_i\left(\theta \right)$ to the arriving wave from the angle $\theta$ of element i and introducing the evaluation function shown in Equation (3). $\theta$ when $P\left(\theta \right)$ becomes maximum is consistent with the direction of arrival.

$$P\left(\theta \right)={\displaystyle \frac{1}{{\left|\frac{E\left[{\mathit{x}}_1{\mathit{x}}_2^{*}\right]}{\sqrt{E\left[{|\mathit{x}}_1{\left.\right|}^2\right]E\left[|{\mathit{x}}_2{\left.\right|}^2\right]}}-\frac{a_1\left(\theta \right)a_2{\left(\theta \right)}^{*}}{\left|a_1\left(\theta \right)\right|\left|a_2\left(\theta \right)\right|}\right|}^2}}$$

(3)

$$[a_1\left(\theta \right),a_2\left(\theta \right)]=\left[1,exp\left(-{\displaystyle \frac{2\pi d_{\mathrm{a}\mathrm{r}\mathrm{r}\mathrm{a}\mathrm{y}}}{\lambda }}\mathit{sin}\theta \right)\right]$$

(4)

The estimated direction of arrival is then used to estimate the location of the radio source by triangulation, as shown in Figure 4.

This method provides highly accurate estimation when there is a line of sight (LoS) from the antenna to the source of the radio wave. However, when estimating the direction of arrival in urban areas, multiple waves are received due to reflected waves. The problem with this method is that the direction of arrival cannot be estimated properly at that time.

2.2. Estimation Method Using the Beamformer Method

The Beamformer method is a basic method for estimating the direction of arrival. It estimates the direction of a radio source by scanning the main lobe of a directional antenna with high gain in all directions as shown in Figure 5, outputting the output voltage change of the array at each angle, and obtaining the maximum output value.

For simplicity, the direction-of-arrival estimation method using the Beamformer method is described in detail using a K-element linear array antenna as shown in Figure 6 [9].

In this case, the output of the array is $y\left(t\right)$ can be expressed as follows.

$$y\left(t\right)={\mathit{W}}^H\mathit{X}\left(t\right)$$

(5)

$$\mathit{X}\left(t\right)={\left[x_1\left(t\right),x_2\left(t\right),\dots ,x_K\left(t\right)\right]}^T$$

(6)

$$\mathit{W}={\left[w_1,w_2,\dots ,w_K\right]}^T$$

(7)

And the output voltage $P_{out}$ is provided as follows.

$$P_{out}={\displaystyle \frac{1}{2}}E\left[{\left|y\left(t\right)\right|}^2\right]={\displaystyle \frac{1}{2}}{\mathit{W}}^H{\mathit{R}}_{xx}\mathit{W}$$

(8)

$${\mathit{R}}_{xx}=E\left[\mathit{X}\left(t\right){\mathit{X}}^H\left(t\right)\right]$$

(9)

Assuming that $L$ waves are incident on the incoming wave and that the respective signal waveforms and angles of arrival are $F_l\left(t\right)$ and ${\theta }_l\left(l=\mathrm{1,2},\dots ,L\right)$ the respective direction vectors ${\mathit{V}}_l$ can be expressed as follows.

$${\mathit{V}}_l={\left[\mathrm{e}\mathrm{x}\mathrm{p}\left(-j{\displaystyle \frac{2\pi }{\lambda }}{d}_1\sin {\theta }_l\right),...,\mathrm{e}\mathrm{x}\mathrm{p}\left(-j{\displaystyle \frac{2\pi }{\lambda }}{d}_K\sin {\theta }_l\right)\right]}^T\equiv \mathit{a}\left({\theta }_l\right)$$

(10)

Also $d_k(k=\mathrm{1,2},...,K)$ is the distance from the reference array to the position of the k-th element. At this time, the input vector is expressed as follows.

$$\mathit{X}\left(t\right)=\sum _{l=1}^L{F}_l\left(t\right)\mathit{a}\left({\theta }_l\right)+\mathit{N}\left(t\right)$$

(11)

$$=\mathit{A}\mathit{F}\left(t\right)+\mathit{N}\left(t\right)$$

(12)

$$\mathit{A}=\left[\mathit{a}\left({\theta }_1\right),\mathit{a}\left({\theta }_2\right),\dots ,\mathit{a}\left({\theta }_L\right)\right]$$

(13)

$$\mathit{F}\left(t\right)={\left[F_1\left(t\right),F_2\left(t\right),\dots ,F_L\left(t\right)\right]}^T$$

(14)

In the above equation, $\mathit{N}\left(t\right)$ is the thermal noise vector, whose components are independent complex Gaussian processes with zero mean and variance (power) ${\sigma }^2\left(=P_n\right)$. Also, the correlation matrix ${\mathit{R}}_{xx}$ is expressed as follows.

$${\mathit{R}}_{xx}=E\left[\mathit{X}\left(t\right){\mathit{X}}^H\left(t\right)\right]$$

$$=\mathit{A}E\left[\mathit{F}\left(t\right){\mathit{F}}^H\left(t\right)\right]{\mathit{A}}^H+E\left[\mathit{N}\left(t\right){\mathit{N}}^H\left(t\right)\right]$$

$$=\mathit{A}\mathit{S}{\mathit{A}}^H+{\sigma }^2\mathit{I}$$

(15)

$$\mathit{S}\triangleq E\left[\mathit{F}\left(t\right){\mathit{F}}^H\left(t\right)\right]$$

(16)

The matrix $\mathit{S}$ in Eq. (16), which represents the correlation between arriving waves, is called the signal correlation matrix and is a diagonal matrix with $P_l\left(l=\mathrm{1,2},\dots ,L\right)$ (each input power) in the diagonal components if all arriving waves are uncorrelated with each other. In order to direct the main lobe of the array antenna to angle $\theta$ using the above equation and others, the weight $w_k$ expressed in Equation (17), which is multiplied to each input signal based on the common-phase condition (the condition to align the phases to be in-phase), can be expressed as follows.

$$w_k=\mathrm{e}\mathrm{x}\mathrm{p}\left(-j{\displaystyle \frac{2\pi }{\lambda }}{d}_k\sin \theta \right) \left(k=\mathrm{1,2},\dots ,K\right)$$

(17)

The angle $\theta$ is varied from $-\pi /2$ to $\pi /2$ to find the peak of the output voltage of the array. The output voltage of the array is expressed by Equation (18), and the angular distribution by the Beamformer method is obtained by normalizing Equation (18) as follows.

$$P_{BH}\left(\theta \right)={\displaystyle \frac{P_{out}}{{\mathit{a}}^H\left(\theta \right)\mathit{a}\left(\theta \right)/2}}={\displaystyle \frac{{\mathit{a}}^H\left(\theta \right){\mathit{R}}_{xx}\mathit{a}\left(\theta \right)}{{\mathit{a}}^H\left(\theta \right)\mathit{a}\left(\theta \right)}}$$

(18)

The Beamformer method estimates the direction of arrival and the input power of the arriving wave from the position of the peak of the evaluation function as shown in Equation (18) using the input correlation matrix ${\mathit{R}}_{xx}$ and the weight expressed in Equation (17). Using the estimated direction of arrival, the radio wave source is estimated by triangulation as shown in Figure 4.However, in the case of multiple incoming radio waves, if the main lobe of the antenna pattern shown in Figure 5 faces a certain radio wave while the side lobes (in the direction of relatively high gain) coincide with the direction of other radio waves, the output voltage will include not only one wave but also other waves, it is not guaranteed that the evaluation function may always have peak in the desired direction of arrival. In addition, when radio waves with almost similar arrival angles are incident, it may not be possible to distinguish the two waves when scanning with the main lobe having a wide angle of arrival as shown in Figure 5. In such cases, the angular resolution of location estimation by the Beamformer method becomes low, and the peak value may not accurately represent the power of the signal.

In fact, as shown in [10] which studied the direction-of-arrival estimation using the Beamformer method in a multipath environment, the estimation accuracy is degraded depending on the number of arriving waves. There is also the Capon method, which directs the main lobe in one direction while at the same time attempting to reduce the output contribution from other directions. However, all of these methods use the main lobe of the antenna pattern, and the angular resolution of the estimation is limited. Therefore, it is difficult to apply the Beamformer method to the scenarios of ground emitters in urban areas where many radio waves arrive.

2.3. Estimation Method Using the MUSIC Method

This section describes in detail the MUSIC method, which uses eigenvalues and eigenvectors to estimate the direction of arrival at a higher angle resolution than the Beamformer method [9]. For the sake of simplicity, let assume a K-element linear array antenna as shown in Figure 6 is used. In this case, the input vector $\mathit{X}\left(t\right)$ of each antenna is the same as in Equation (12).

$$\mathit{X}\left(t\right)=\mathit{A}\mathit{F}\left(t\right)+\mathit{N}\left(t\right)$$

(19)

The correlation matrix of the input vectors is then as in Equation (15),

$${\mathit{R}}_{xx}\triangleq E\left[\mathit{X}\left(t\right){\mathit{X}}^H\left(t\right)\right]=\mathit{A}\mathit{S}{\mathit{A}}^H+{\sigma }^2\mathit{I}$$

(20)

where the signal correlation matrix $\mathit{S}$ can be expressed as follows.

$$\mathit{S}=\left[\begin{array}{cccc}E\left[{\left|F_1\left(t\right)\right|}^2\right]& E\left[F_1\left(t\right)F_2^{*}\left(t\right)\right]& \cdots & E\left[F_1\left(t\right)F_L^{*}\left(t\right)\right]\\ E\left[F_2\left(t\right)F_1^{*}\left(t\right)\right]& E\left[{\left|F_2\left(t\right)\right|}^2\right]& \cdots & E\left[F_2\left(t\right)F_L^{*}\left(t\right)\right]\\ ⋮& ⋮& \ddots & ⋮\\ E\left[F_L\left(t\right)F_1^{*}\left(t\right)\right]& E\left[F_L\left(t\right)F_2^{*}\left(t\right)\right]& \cdots & E\left[{\left|F_L\left(t\right)\right|}^2\right]\end{array}\right]$$

(21)

In the absence of thermal noise, if the arriving waves are uncorrelated, $\mathit{S}$ becomes a diagonal matrix and its rank is $L$. If the direction matrix $\mathit{A}$ also differs from the direction of arrival, its column vectors become independent, and its rank becomes $L$. Therefore, in this case, the correlation matrix ${\mathit{R}}_{xx}$ is a non-negative definite Hermitian matrix of rank $L$. Denoting the eigenvalues of this matrix by ${\mu }_i\left(i=\mathrm{1,2},\dots ,K\right)$ and the corresponding eigenvectors by ${\mathit{e}}_i\left(\mathrm{1,2},\dots ,K\right)$

$$\mathit{A}\mathit{S}{\mathit{A}}^H{\mathit{e}}_i={\mu }_i{\mathit{e}}_i\left(i=\mathrm{1,2},\dots K\right)$$

(22)

and each eigenvalue can be expressed as a real number.

$${\mu }_1\ge {\mu }_2\ge \cdots {\mu }_L>{\mu }_L+1=\cdots ={\mu }_K=0$$

(23)

And the eigenvectors are orthogonal to each other as follows.

$${\mathit{e}}_i^H{\mathit{e}}_k={\delta }_{ik} \left(i,k=\mathrm{1,2},\dots ,K\right)$$

(24)

In the presence of thermal noise, the thermal noise power is added to the eigenvalues and the eigenvectors remain unchanged.

$${\mathit{R}}_{xx}{\mathit{e}}_i=\left(\mathit{A}\mathit{S}{\mathit{A}}^H+{\sigma }^2\mathit{I}\right){\mathit{e}}_i$$

$$=\mathit{A}\mathit{S}{\mathit{A}}^H{\mathit{e}}_i+{\sigma }^2{\mathit{e}}_i$$

$$={\mu }_i{\mathit{e}}_i+{\sigma }^2{\mathit{e}}_i$$

$$=\left({\mu }_i+{\sigma }^2\right){\mathit{e}}_i (i=\mathrm{1,2},...K)$$

(25)

Then we have

$${\lambda }_i\triangleq {\mu }_i+{\sigma }^2 \left(i=\mathrm{1,2},\dots ,K\right)$$

(26)

Denoting the eigenvalues of the correlation matrix by

$${\lambda }_1\ge {\lambda }_2\ge \cdots \ge {\lambda }_L>{\lambda }_{L+1}=\cdots ={\lambda }_K={\sigma }^2$$

(27)

the number of arriving waves L can be obtained by finding the number of eigenvalues of the correlation matrix that have values larger than the noise power. Also, if we focus on the eigenvectors corresponding to the eigenvalues equal to the thermal noise power,

$${\mathit{R}}_{xx}{\mathit{e}}_i=\left(\mathit{A}\mathit{S}{\mathit{A}}^H+{\sigma }^2\mathit{I}\right){\mathit{e}}_i={\lambda }_i{\mathit{e}}_i={\sigma }^2{\mathit{e}}_i \left(i=L+1,\dots ,K\right)$$

(28)

based on Equation (20), we can derive the following property.

$$\mathit{A}\mathit{S}{\mathit{A}}^H{\mathit{e}}_i=0 \left(i=L+1,\dots ,K\right)$$

(29)

Also, from the fact that both $\mathit{A}$ and $\mathit{S}$ are full ranks.

$${\mathit{A}}^H{\mathit{e}}_i=0 \left(i=L+1,\dots ,K\right)$$

(30)

In other words

$${\mathit{a}}^H\left({\theta }_l\right){\mathit{e}}_i=0 \left(l=\mathrm{1,2},\dots ,L;i=L+1,\dots ,K\right)$$

(31)

From this equation, it can be seen that the eigenvectors corresponding to eigenvalues equal to the thermal noise power are orthogonal to the direction vectors of all arriving waves. Using this property, the evaluation function of the MUSIC method can be expressed as follows with normalization.

$$P_{M}U\left(\theta \right)\triangleq {\displaystyle \frac{1}{{\sum }_{i=L+1}^K{\left|{\mathit{e}}_i^H\mathit{a}\left(\theta \right)\right|}^2}}\times {\mathit{a}}^H\left(\theta \right)\mathit{a}\left(\theta \right)$$

(32)

$$={\displaystyle \frac{{\mathit{a}}^H\left(\theta \right)\mathit{a}\left(\theta \right)}{{\mathit{a}}^H\left(\theta \right){\mathit{E}}_N{\mathit{E}}_N^H\mathit{a}\left(\theta \right)}}$$

(33)

$${\mathit{E}}_N\triangleq \left[{\mathit{e}}_{L+1},\dots ,{\mathit{e}}_K\right]$$

(34)

The denominator part of Eq. (32) shows a minimum value at $\theta$, which is the direction vector (direction vector of the arriving wave) orthogonal to the eigenvector corresponding to the noise, as mentioned earlier. The MUSIC method, on the other hand, uses the noise portion (the portion where the power is small) to achieve high resolution estimation. However, when multiple radio waves are incident from a single source due to reflected waves, etc., there is a correlation between the incoming waves. In this case, the rank of the signal correlation matrix expressed in Equation (21) is not equal to the number of arriving waves, and the direction of arrival cannot be accurately estimated.

It is known that all of the three location estimation methods using direction-of-arrival estimation introduced in this section have a problem that the estimation accuracy is degraded by many reflected waves when a ground antenna is used.

3. Numerical Analysis and Results

Currently, location estimation using HAPS is not yet available and is considered to be operated as a base station. Because of this situation, no demonstration experiment of location estimation using HAPS has been realized. In this section, methods of direction-of-arrival estimation will be applied for HAPS for location estimation of ground emitter. Since the phase difference method is just a primitive version of the Beamformer one, only the Beamformer and MUSIC methods are introduced in this section for evaluation of the localization accuracy.

3.1. Comparison of LoS rate between HAPS and UAV

In Section 2, we discussed various location estimation methods. However, the accuracy of such methods varies greatly depending on whether the environment between the receiving antenna and the radio source is a LoS environment or an NLoS environment. In this section, we perform a ray tracing simulation using the SBR (Shooting and Bouncing Rays method) in MATLAB to compare the LoS rate between the case using HAPS and the case using a UAV. In this analysis, the case in which a direct wave (a radio wave that is not reflected by any obstacle) reaches the receiving antenna was considered as the LoS environment and vice versa is considered an NLoS environment.

3.1.1. Simulation Model

Simulations were performed using a 3D model environment that mimics the town of Chicago that pre-exists in MATLAB, shown in Figure 7.

An overall view of the simulation is shown in Figure 8. The simulation was conducted assuming that the HAPS was 2,000 m above the ground and the UAV was 200 m above the ground.

For simplicity, we will use the simulation of the LoS rate over the town to illustrate the simulation. First, random radio sources are generated in the town as shown in Fig. 9. Here, the height of each antenna was set to 10 meters.

Figure 9. Randomly generated radio wave sources.

Figure 9. Randomly generated radio wave sources.

Next, as shown in Figure 10, a receiving antenna (in this case, a UAV) is generated over the town, the number of direct waves is calculated using a ray tracing simulation, and the LoS rate is calculated by calculating the ratio of the number of direct waves to the number of radio wave sources. In this simulation, the number of radio wave sources is set to 1000. The same simulation is performed for HAPS.

When the simulation is performed away from the center, receivers are generated in several directions as shown in Figure 11, and the average value of the LoS ratio obtained from each receiver is used as the LoS ratio.

3.1.2. Results

The results of the simulation conducted above are shown in Figure 12. The LoS ratio of the HAPS was about 100% when the center of the urban area where the radio source was generated was directly below the HAPS, and about 60% at a distance of 20 km (elevation angle 45°). On the other hand, when a UAV was used, the LoS ratio was about 30% even when the center of the urban area was directly below the UAV, indicating that the majority of the radio sources were out of sight. Therefore, an altitude of about HAPS is considered necessary for position estimation by angle estimation.

3.2. Location Estimation Simulation

In order to realize a simulation of radio source location estimation using HAPS, we used MATLAB to estimate the location of radio sources. Several location estimation methods were discussed in Section 2 but considering the wide search area and required estimation accuracy, the Beamformer method and the MUSIC method were used in our simulation. As a comparison, a case in which location estimation was performed at the altitude of the UAV was also simulated. In our analysis, the location estimation simulation was performed for a non-mobile illegal radio source located outdoors in an urban area, as shown in Figure 7.

The simulation process is shown in Figure 13. (1) A radio source is randomly generated, (2) A ray traccing simulation is performed to obtain information about arriving radio waves (direction of arrival, number of reflections, propagation distance, etc.), and (3) Received signal is generated from the aforementioned information. The direction of arrival is estimated by using each estimation method applied on the received data generated in (3). Then, location estimation is performed using triangularization method. By comparing the randomly generated coordinates of the radio source with the estimated coordinates, the location estimation of the radio source using HAPS and other methods is evaluated. Figure 13 shows an overview of the simulation.

3.2.1. Receiving and Transmitting Antenna Configuration

Table 1 shows the antenna specifications for estimating the location of a radio wave source using ray tracing simulation. Hereafter, the transmitting antenna is denoted as Tx and the receiving antenna as Rx.

The receiving antenna was modeled after the cylinder antenna used in HAPS as a base station for the simulation. The actual prototype cylinder antenna can be found in [12].

In this paper, the cylinder antenna consists of $N=$ 36 elements as depicted in Figure 14, where each element was assumed to be omni-directional and the maximum gain was set to 0dBi in the simulation.

3.2.2. Radio propagation model

This section describes the radio propagation model used in the simulation of radio source location estimation using HAPS. In this simulation, statistical data were not used in the propagation model because a ray tracing simulation was conducted to determine specific radio waves. The free propagation loss was modeled using Friis formula as follows,

$${PL}_{fsp}=10lo{g}_{10}{\left({\displaystyle \frac{4\pi d}{\lambda }}\right)}^2\left[dB\right]$$

(37)

where $d$ and $\lambda$ are the distance between the transceiver pairs and the wavelength of the considering radio wave, respectively.

The overall path loss can be expressed by Equation (38), where the additional attenuation $L_{ref}$ due to earth and building reflections was calculated using the Fresnel equation. The permittivity and conductivity of the surface materials used in these calculations were taken from the ITU-R study [13]. Assuming that both the ground surface and the building are concrete, the relative permittivity and conductivity were calculated as 5.31 and 0.0548 [S/m], respectively.

$$PL=P{L}_{fspl}+P{L}_{ref} \left[dB\right]$$

(38)

3.2.3. Signal processing system

This section describes the signal processing performed in the location estimation simulation. First, the analytical signal input to the cylinder antenna is obtained by generating a single radio source in a 3D model and considering the case where L radio waves are incident from the source. The propagation loss $P{L}_l\left[dB\right] (l=\mathrm{1,2},\dots ,L)$ is obtained from Equation (39) using the propagation distance and other data obtained from the ray tracing simulation, and the input power is calculated by following formula.

$$P_l=10\times {10}^{-\left({\displaystyle \frac{P{L}_l}{10}}\right)} \left[W\right] (l=\mathrm{1,2},\dots ,L)$$

(39)

Using this value, the propagation time delay against the direct wave ${\tau }_l(l=\mathrm{1,2},\dots ,L)$, $({\tau }_1=0)$and the times of reflection $m_l(l=\mathrm{1,2},\dots ,L)$, the complex waveform of the l-th wave can be written by the following formula.

$$F_l\left(t\right)=\sqrt{P_l}\mathit{exp}\left\{j2\pi f\left(t-{\tau }_l\right)\right\}·{\left(-1\right)}^{m_l} \left(l=\mathrm{1,2},\dots L\right)$$

(40)

The array matrix for the cylinder antennas is then expressed by the following equation [9].

$$\mathit{A}=\left[\mathit{a}\left({\theta }_1,{\varphi }_1,f\right),\dots ,\mathit{a}\left({\theta }_L,{\varphi }_L,f\right)\right]$$

(41)

$$\mathit{a}\left({\theta }_l,{\varphi }_l,f\right)={\left[\exp \left\{j2\pi {\displaystyle \frac{f}{c}}{\mathit{r}}_1^T\mathit{L}\left({\theta }_l,{\varphi }_l\right)\right\} ,\dots ,\exp \left\{j2\pi {\displaystyle \frac{f}{c}}{\mathit{r}}_N^T\mathit{L}\left({\theta }_l,{\varphi }_l\right)\right\} \right]}^T(l=\mathrm{1,2},\dots ,L)$$

(42)

$$\mathit{L}\left(\theta ,\varphi \right)\triangleq \left[\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{c}\mathrm{o}\mathrm{s}\varphi ,\mathrm{s}\mathrm{i}\mathrm{n}\theta \mathrm{s}\mathrm{i}\mathrm{n}\varphi ,\mathrm{c}\mathrm{o}\mathrm{s}\theta \right]$$

(43)

where ${\mathit{r}}_k(k=\mathrm{1,2},\dots ,N)$ is the position vector of each element, and $\theta$ and $\varphi$ are the elevation and azimuth angles to the reference point, respectively.

The thermal noise power $P_N$ generated at the receiving antenna is given by,

$$P_N=10{\mathit{log}}_{10}\left(kTB\right)+NF \left[dBW\right]$$

(44)

where $k, T, B, NF$ are Boltzmann's constant, absolute temperature of the noise source, bandwidth, and noise figure, respectively. The absolute temperature in the HAPS environment was assumed to be 216.5 [K] and the bandwidth was 10 [MHz]. The noise figure was calculated as 3 dB [14]. As a result, the thermal noise power $P_N$was about $-102$ [dBm] in HAPSs environment while in the UAV environment, the thermal noise power was slightly higher since the absolute temperature at the UAV altitude is about 288 [K]. Using Equations (40), (41), and (44), the input vector $\mathit{X}\left(t\right)$ to each antenna is expressed by the following equation.

$$\mathit{X}\left(t\right)=\mathit{A}\mathit{F}\left(t\right)+\mathit{N}\left(t\right)$$

(45)

$$\mathit{F}\left(t\right)={\left[F_1\left(t\right),F_2\left(t\right),\dots ,F_L\left(t\right)\right]}^T$$

(46)

In the above equation,$\mathit{N}\left(t\right)$is the thermal noise vector, whose components are independent complex Gaussian processes with mean zero and variance (power) of $P_N$. The input vector obtained by substituting appropriate values for t in Equation (45) for different time samples, was used as the input complex signal for the cylinder array. The number of samples was set to 100 in this paper. By applying the aforementioned Beamformer and MUSIC methods to the obtained input vector $\mathit{X}\left(t\right)$, the direction of arrival from the radio source was estimated and used for the location estimation method described in the next section. The direction of arrival with the highest power obtained by the Beamformer and MUSIC methods was used as the estimated direction of arrival.

3.2.4. Location Estimation Method

This section describes a method for obtaining the estimated coordinates of a radio source from the direction of arrival obtained by the MUSIC and Beamformer methods. In this simulation, the direction of arrival of radio waves is estimated from three points as shown in Figure 15.

The latitude and longitude of each HAPS observation point are converted to the Cartesian coordinate system, and a direction line is created in the x-y coordinate system from the direction of arrival. Figure 16 shows two localization methods from the estimated angular information. In the left-hand side case that all angular information is available from the three sensors, the center of gravity of the intersections of each directional line is considered as the estimated coordinates of the emitter. In the right-hand side case that the angular information from one of three sensors is missing, the estimated coordinates are the intersection points of the direction lines of the remaining sensors where angular information is available. For all other cases, the location of the emitter is considered to be unidentified (undetected).

3.2.5. Results

The 3D city model shown in Figure 7 was used to randomly generate 1,000 radio wave sources, and the location of each source was estimated using the method described in the previous section. The horizontal distance between the center of the urban area where there exist illegal radio sources and the hovering sensors was varied as different scenarios in our below evaluation results. Three scenarios of sensor altitudes i.e. HAPS, UAV and DEURAS are also separately evaluated in this section.

▪

HAPS altitude

Table 3 shows the results of the localization using HAPS at three points 5 km away from the urban area. The Beamformer and MUSIC methods were used to estimate the direction of arrival. The estimation probability denotes the percentage that the emitter can be identified (with localization errors) via our proposed system as mentioned in Figure 16. The stochastics of the estimation errors were derived only from emitters which were identified by our proposed method. From Table 3, it is observed that there was no significant difference in the results between the Beamformer and MUSIC methods, and 90% of the estimated locations were within 30 meters of the true coordinates.

Table 4 shows the results of the estimation using HAPS at three locations 10 km away from the urban area. Since the sensors are more distanced from the sources compared to the 5-km case, we can observe a slight degradation in terms of localization accuracies as compared to these in Table 3.

Furthermore, Table 5 shows the results of the estimation using HAPS at three locations 20 km away from the urban area. Similar to the two aforementioned scenarios, neither the Beamformer nor the MUSIC method showed a significant difference in results. Also, the accuracy of the estimation was poorer at greater distances compared to the cases at 5km and 10 km. Furthermore, only about 83.6% of the radio transmitting sources in urban areas could be obtained.

The CDF of localization errors for these scenarios is summarized in Figure 17.

▪

UAV altitude

Similarly, a simulation was conducted assuming a UAV with the height of the receiving antenna changed to 200 meters. In addition, because the effect of reflections is likely to increase at the altitude of the UAV, the simulation was performed considering radio waves that can be reflected up to four times. Table 6, Table 7 and Table 8 show the results of the estimation using UAVs at three locations 5 km, 10km and 20km away from the urban area, respectively.

In Table 6, only 20% of the total number of radio sources could be identified due to the low LoS ratio as described in Section 3-1. The localization estimation error of the radio sources that could be identified was also larger than that of the HAPS case.

Table 7 shows the results of the estimation using UAVs at three locations 10 km away from the urban area. As in the 5-km case, the number of radio sources that could be identified was small. However, the estimation error and the CDF90% value was slightly improved against the 5 km case. This is thought to be due to the fact that the 10 km environment has smaller multipath components compared to the 5km case, such that angular estimation has higher accuracy.

Table 8 shows the results of the estimation using UAVs at three locations 20 km away from the urban area. Both the estimated probability and localization accuracy are degraded in this case. It is thought to be due to the fact that the sensors are located too far away such that the arrival waves after multiple reflections are too weak compared to the noise level.

The CDF of localization errors for these different scenarios is summarized in Fig. 18.

Figure 18. CDF of errors using UAV at each distance.

Figure 18. CDF of errors using UAV at each distance.

▪

DEURAS

Finally, a model simulating the DEURAS system was created by installing the receiving antennas on the rooftop of buildings with good visibility in an urban area, as shown in Figure 19. Other system parameters were the same as those for UAVs.

Table 9 shows the evaluation results of position estimation simulations using the Beamformer and MUSIC methods where DEURAS sensors are located as shown in Figure 19. Compared to location estimation using UAVs, the estimation accuracy is still poorer, except for the fact that DEURAS can provide a higher estimated probability since the sensors are close to the sources.

For a fair comparison, we conducted new simulations that HAPS and UAVs were deployed at three points at a distance of 500 m from the center of the urban area (just right above the area similar to the case of DEURAS system).

Table 10 shows the evaluation results of the estimation using HAPS at three locations 500 m away from the urban area. Similarly, Table 11 shows the results of the evaluation when estimating using UAVs at three locations 500 m away from the urban area.

Figure 20 shows summarized the CDF of localization errors of all scenarios. These results show that HAPS shows the best performance in terms of both estimated probability and localization accuracy. It is owing to the fact that HAPS has higher LoS ratio to guarantee higher estimated probability, and also direct waves are more likely to arrive at HAPS rather than multiple reflection paths that increase the angular estimation accuracy.

4. Conclusions

In this study, we compared HAPS, UAV and DEURAS systems for location estimation of radio sources, and found that the HAPS estimation showed excellent accuracy at distances up to 20 km due to its high LoS ratio, while the UAV estimation showed poor performance due to interference from reflected waves, especially in urban environments. The performance of the UAV estimation was low, especially in urban environments, due to interference from reflected waves. The position estimation using an antenna configuration that mimics the DEURAS system also showed lower position estimation accuracy compared to HAPS due to the antenna height and reflected waves. Overall, HAPS is excellent for monitoring radio sources over a wide area, and we believe that integrating HAPS into the existing DEURAS setup has the potential to enhance the current system.

In this paper, we used MATLAB to simulate the location estimation of a radio source, it is possible that diffracted and scattered waves, in addition to other arriving and reflected waves, may exist in a real environment, and these may affect the estimation results. Therefore, we would like to conduct actual experiments in the future to confirm the accuracy of the estimation. In addition, we would like to consider the location estimation of low-power radio sources, moving radio sources, and indoor radio sources, which were not covered in this study.