*Article* **Continuous Hidden Markov Model Based Spectrum Sensing with Estimated SNR for Cognitive UAV Networks**

**Yuqing Feng <sup>1</sup> , Wenjun Xu 2,\*, Zhi Zhang <sup>1</sup> and Fengyu Wang <sup>3</sup>**


**Abstract:** In this paper, to enhance the spectrum utilization in cognitive unmanned aerial vehicle networks (CUAVNs), we propose a cooperative spectrum sensing scheme based on a continuous hidden Markov model (CHMM) with a novel signal-to-noise ratio (SNR) estimation method. First, to exploit the Markov property in the spectrum state, we model the spectrum states and the corresponding fusion values as a hidden Markov model. A spectrum prediction is obtained by combining the parameters of CHMM and a preliminary sensing result (obtained from a clustered heterogeneous two-stage-fusion scheme), and this prediction can further guide the sensing detection procedure. Then, we analyze the detection performance of the proposed scheme by deriving its closed-formed expressions. Furthermore, considering imperfect SNR estimation in practical applications, we design a novel SNR estimation scheme which is inspired by the reconstruction of the signal on graphs to enhance the proposed CHMM-based sensing scheme with practical SNR estimation. Simulation results demonstrate the proposed CHMM-based cooperative spectrum sensing scheme outperforms the ones without CHMM, and the CHMM-based sensing scheme with the proposed SNR estimator can outperform the existing algorithm considerably.

**Keywords:** cognitive UAV networks; clustered two-stage-fusion cooperative spectrum sensing; continuous hidden Markov model; SNR estimation

## **1. Introduction**

With the advantage of high flexibility and low deployment cost, unmanned aerial vehicles (UAVs) have been widely used in military communications, weather monitoring, emergency rescue [1] and some other UAV-assisted Internet of Things (IoT) applications [2]. The large-scale deployment of UAVs has exacerbated the shortage of spectrum resources. However, the existing spectrum allocation strategies cannot effectively use the scarce spectrum resources, which becomes the bottleneck for enhancing the communication performance of UAVs [3]. Cognitive radio (CR) is proposed to solve the problem, which improves the spectrum efficiency by perceiving spectrum holes and providing secondary UAVs with opportunities to reuse idle spectrum. Thus, CR can further guide the spectrum utilization in cognitive unmanned aerial vehicle networks (CUAVNs), including the resource allocation for low-latency communications [4], high-quality services with limited resources [5], maximum achievable throughput [6], and optimal power allocation. To enable CUAVNs, accurate spectrum sensing attaches great importance.

With the development of UAVs, spectrum sensing in unmanned aerial vehicle networks (UAVNs) has attracted attention from both academia and industry. The detection performance is enhanced by using multiple secondary UAVs in [7,8]. Authors in [3,9]

**Citation:** Feng, Y.; Xu, W.; Zhang, Z.; Wang, F. Continuous Hidden Markov Model Based Spectrum Sensing with Estimated SNR for Cognitive UAV Networks. *Sensors* **2022**, *22*, 2620. https://doi.org/10.3390/s22072620

Academic Editor: Omprakash Kaiwartya

Received: 10 January 2022 Accepted: 22 March 2022 Published: 29 March 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/).

consider the combination of UAVs and terrestrial communication equipment for air-groundintegrated spectrum sensing. The heterogeneity of spatial information is further taken into account by using the 3D [10], so as to improve the spectrum sensing detection performance. The above studies introduce information from multiple users and spatial dimensions to improve the spectrum utilization of CUAVNs. However, they do not consider the temporal information of the spectrum, which affects the detection performance of UAVNs. Exploring the temporal correlation of spectrum states by focusing on the Markov property is an effective and novel idea to enhance the detection probability [11].

Another concern is that knowledge of the signal-to-noise ratio (SNR) is required before spectrum sensing [12], e.g., energy detection and cyclostationary feature detection. Therein, it is assumed that SNR is perfectly known. However, UAVs need to move across a large area in lots of applications, and the links between two UAVs that are far apart frequently break and reestablish [13], which degrades the sensing performance and the spectrum utilization. Besides, large flying areas, unstable links, and dynamic network topologies lead to variable SNRs in CUAVNs, which makes the assumption of pre-known SNRs no longer applicable, bringing new challenges for sensing in UAVNs.

To overcome these problems, we propose a continuous hidden Markov model (CHMM) based sensing scheme with a novel space smoothing second- and fourth-order moments (SS-M2M4) SNR estimator. Note that the true states of the unauthorized spectrum are not observable, however, the sensing results originating from the unobservable spectrum states can be easily obtained. Therefore, a hidden Markov chain model fits well with the spectrum sensing procedure. To fully exploit the temporal information and further enhance the detection performance, we combine the spectrum sensing scheme with CHMM. To the best of our knowledge, our work is the first to consider the continuous hidden Markov model in non-centralized CUANs. Moreover, we further provide theoretical analysis of the detection probability and the false alarm probability, while previous works mainly focus on numerical simulations. Besides, for generalized estimation and higher accuracy, based on a general SNR estimation method of a second- and fourth- moments (M2M4) estimator [14], we propose a SS-M2M4 SNR estimator. The SNR smoothness of the neighboring CUAVs has been taken into account, enlightened by the smoothness of the graph signal [15]. Compared to the widely-used M2M4 estimator, our proposed SNR estimator achieves more accurate estimation. With the SNR provided by the proposed estimator, our CHMM-based spectrum sensing scheme can achieve a high detection probability. Our contributions are summarized as follows:


The remainder of this paper is organized as follows. Section 2 discusses the related work. Section 3 introduces the system model. In Section 4, the CHMM based spectrum sensing scheme is proposed. In Section 5, the closed-form detection probability and false alarm probability expressions of the proposed method are derived. In Section 6, we design a novel SNR estimator to enhance the proposed CHMM-based spectrum sensing scheme with a more accurate SNR. The simulation results are presented to evaluate the proposed algorithms in Section 7. Section 8 concludes the work and discusses the possible future work.

#### **2. Related Work**

#### *2.1. Hidden Markov Model (HMM) Based Spectrum Sensing*

Hidden Markov model (HMM) based spectrum sensing means using the HMM to model the spectrum sensing procedure. As the HMM has a hidden layer and an observable layer, it fits well with the sensing procedure whose spectrum states are unknown but the receiver can be obtained. Exploring the Markov property of the spectrum states is an effective way to enhance the sensing performance [11]. Compared with deep learning based spectrum sensing [16,17], CHMM-based spectrum sensing has a stronger interpretability and smaller delay, and the initial probability distribution can be used to calculate the initial spectrum utilization. In addition, the obtained parameters can also be used for digital twinning of the communication system. It has been validated that the spectrum state can be modeled as a Markov chain by analyzing real-world measurements [18]. A hidden Markov model based scheme [19] is proposed to predict the arrival of the primary user (PU). Authors in [20] evaluate the reliability of HMM based cooperative spectrum sensing in cognitive radio networks, in the presence of random malfunctioning of secondary user nodes participating in the process. Occupancy prediction schemes based on a discrete hidden Markov model (DHMM) and a continuous hidden Markov model are investigated in [21–24], respectively. Authors in [21,22] adopt DHMM to model the spectrum sensing procedure, which do not make full use of the information obtained by the receiver. Authors in [23,24] use CHMM to model the sensing procedure of CUAVNs. However, they do not consider the dynamicity of UAVs, which is an important characteristic of the UAV networks [25]. Besides, centralized spectrum sensing methods in [23,24] do not work well in CUAVNs since the global fusion center is usually unreachable to secondary CUAVs. Thus, in this paper, we propose a CHMM-based spectrum sensing method to enhance the spectrum utilization in CUANs.

#### *2.2. SNR Estimation*

SNR estimation refers to the calculation of SNR by using signal information [14]. Various algorithms require SNR estimation for optimal performance if the SNR is not constant, such as linear diversity combining techniques and Viterbi algorithms with softdecision [14]. Note that the knowledge of SNR is also required for typical commonly-used spectrum sensing [12]. However, in practical UAV applications, it is difficult to obtain accurate SNR. Conventionally, SNR estimators require knowledge of the signal or the channel, such as the maximum likelihood (ML) SNR estimator [26] and particle swarm optimization (PSO) SNR estimator, based on parameters of hardwares or channels [27]. In addition, there exist estimators designed for specific signals or specific feature spectrum sensing methods, such as the SNR estimator [28] for M-ary amplitude phase shift keying (M-APSK) modulated signals, the SNR estimator for signal with Polar code [29], and the estimator for eigenvalue-based spectrum detectors [30]. Lacking the prior knowledge of the signal/channel and design for specific signals/sensing methods make it hard to generally adopt the above estimators to various CUAVNs. Besides, some scholars have paid attention to deep learning based SNR estimation methods, such as the convolutional neural networks (CNN)—long short term memory (LSTM) based SNR estimators [31] and the CNN-based SNR estimators designed for UAVNs [32]. However, in the spectrum sensing of CUAVNs, the deep learning based methods are too complicated and take more time, leading to the sensing term being missed. Taking all the above into consideration, we design a generalized and low-complexity SNR estimator named SS-M2M4. Based on the proposed SNR estimator, the CHMM-based spectrum sensing method can further enhance detection performance.

#### **3. System Model**

In order to detect whether the authorized spectrum of PU is occupied, a spectrum sensing method is adopted by multiple UAVs [3,33]. In this paper, we propose a CHMMbased spectrum sensing method with consideration of the imperfect SNR estimation to

enhance the detection probability. The overview of this unified scheme is shown in Figure 1. First, to obtain the single-time spectrum sensing fusion result, we adopt the max-min distance clustering algorithm [34] and heterogeneous two-stage-fusion spectrum sensing scheme [35], which are described in Sections 3.1 and 3.2, respectively. Then, to fully take advantage of the temporal correlation of the spectrum states, we propose a CHMMbased spectrum sensing method, which will be introduced in Section 4. Considering the imperfect SNR estimation in practical applications, we propose a novel SNR estimator shown in Section 6 to offer the sensing scheme more desirable SNRs. Finally, combining the CHMM-based spectrum sensing and the proposed SNR estimator, we propose the unified CHMM-based spectrum sensing scheme with advanced SNR estimator.

**Figure 1.** The unified framework of the CHMM-based spectrum sensing scheme with the SS-M2M4 estimator.

#### *3.1. Clustering Method*

Considering the fusion delay brought by distributed cooperative spectrum sensing (DCSS), and to further provide the proposed SS-M2M4 SNR estimator with graph topology, a max-min distance clustering algorithm [34] is adopted. It divides CUAVs with similar position and mobility into the same cluster, and selects the one with the highest trust value in each cluster as the cluster head. Thus, we can adopt intra-cluster centralized cooperative spectrum sensing. Due to the similar locations, the UAVs in the same cluster have similar SNRs, which facilitates the SNR estimation. The cluster heads are less than the number of the CUAVs, so the fusion delay will be reduced and, therefore, we adopt inter-cluster distributed cooperative spectrum sensing.

#### *3.2. Heterogeneous Two-Stage-Fusion Spectrum Sensing*

To further improve the sensing performance in CUAVNs, a heterogeneous cooperative spectrum sensing scheme is employed with the clustering outcome, as shown in Figure 2. A symbol table of the notations is shown in Table 1.

**Figure 2.** Heterogeneous Two-Stage-Fusion.

**Table 1.** Notations.


#### 3.2.1. Heterogeneous Spectrum Sensing Based on Clustering

To be specific, in the sensing state, a heterogeneous cooperative spectrum sensing scheme [36] is used according to the clustering result. Here, the "heterogeneous" means two different detection schemes: energy detection and cyclostationary detection [36]. Cluster heads adopt cyclostationary detection, while other secondary UAVs in the same cluster adopt energy detection. Cyclostationary feature detection is adopted to the cluster head, since it has great detection accuracy and can maintain good performance even in environments with low SNR. The cluster members (CMs) adopt energy detection, since it is simple to implement and does not require prior knowledge of the channels. Similarly to [36], we assume that the transmitted signal of the primary user is a sinusoidal signal, in which the n-th sample is expressed as

$$s(n) = ae^{j2\pi f\_c n + \varphi}.\tag{1}$$

where *a* is the amplitude, *f<sup>c</sup>* is the carrier frequency, and *ϕ* is the carrier phase offset. The received signal of the cluster heads is written as

$$y\_{\mathbb{C}}(n) = \begin{cases} \omega(n) & H\_0 \\ s(n) + \omega(n) & H\_1 \end{cases} \tag{2}$$

where *ω*(*n*) denotes the additive Gaussian noise with zero mean and unit variance. *H*<sup>0</sup> indicates that the spectrum is absent, while *H*<sup>1</sup> indicates that the spectrum is occupied. The test statistic for the first-order cyclostationary detection [36] is expressed as

$$T\_{\mathbb{C}} = |\frac{1}{M} \sum\_{n=1}^{M} y\_{\mathbb{C}(n)} e^{-j2\pi an}|,\tag{3}$$

where *α* is the cyclic frequency, M denotes the number of sampling point. The received signal *yE*(*n*) of the CUAV cluster member is modeled as the same as *yC*(*n*)

$$y\_E(n) = \begin{cases} \omega(n)\_\prime & H\_0 \\ s(n) + \omega(n)\_\prime & H\_1 \end{cases} \tag{4}$$

The test statistic for the energy detector [36] is expressed as

$$T\_E = \frac{1}{M} \sum\_{n=1}^{M} |y\_E(n)|^2,\tag{5}$$

where *M* denotes the number of sampling points. Let *Tic* = [*TC<sup>i</sup>* , *TEi*<sup>1</sup> , *TEi*<sup>2</sup> , . . . , *<sup>T</sup>EiKE* ] denotes the vector of the statistics from the *i*th cluster, *TC<sup>i</sup>* is the cyclostationary detection statistics of the *i*th cluster head, and its distribution [36] can be expressed as

$$T\_{\mathbb{C}\_i} \sim \begin{cases} \mathcal{FN}(\sqrt{\frac{2}{\pi \mathsf{M}}}, \frac{\pi - 2}{\pi \mathsf{M}}), & H\_0 \\ \mathcal{FN}(\zeta\_i, 2\gamma\_{\mathbb{C}\_i} + \frac{1}{\mathsf{M}} - \zeta\_i^2), & H\_1 \end{cases} \tag{6}$$

where *ζ<sup>i</sup>* = q 2 *<sup>π</sup><sup>M</sup> e* <sup>−</sup>*MγCi* + p 2*γC<sup>i</sup>* (1 − 2*Q*( p 2*MγC<sup>i</sup>* )), *γC<sup>i</sup>* denotes the SNR of the *i*th cluster head, and FN is the folded normal distribution [37]. *TEij* is the energy detection statistics of the *j*th the node in the *i*th cluster. When M is large enough, *TEij* approximatively follows normal distribution according to the central limit theorem [38], but slightly different to [38], our *TEij* is 1/M times the energy detection statistics in [38], thus the our SNR is *M* times the SNR in [38], and we assume the energy of the received signal is 1. Thus, the distribution of *TEij* can be expressed as

$$T\_{\mathbb{E}\_{ij}} \sim \begin{cases} \mathcal{N}(\gamma\_{\mathbb{E}\_{ij}}^{-1}, \frac{2\gamma\_{\mathbb{E}\_{ij}}^{-2}}{M}), & H\_0 \\ \mathcal{N}(1 + \gamma\_{\mathbb{E}\_{ij}}^{-1}, \frac{2(\gamma\_{\mathbb{E}\_{ij}}^{-2} + 2\gamma\_{\mathbb{E}\_{ij}})^{-1})}{M}), & H\_1 \end{cases} \tag{7}$$

where *γEij* is the SNR of the *j*th node in the *i*th cluster.

#### 3.2.2. Two-Stage-Fusion

In the fusion duration, we adopt a two-stage-fusion scheme [34], which includes the intra-cluster fusion stage and the inter-cluster fusion stage. The CUAV cluster members adopt energy detection, and we apply high-accuracy centralized soft fusion in the intracluster since the performance of energy detection is somewhat not precise. In the intracluster fusion stage, according to the assumption that the observations are independent, we can attain the likelihood ratio test (LRT) [36] of the *i*th cluster

$$\begin{split} T\_i &= \prod\_{j=1}^{K} \frac{P[T\_{ij}|H\_1]}{P[T\_{ij}|H\_0]} \\ &= \prod\_{j=1}^{K\_E} \frac{P[T\_{E\_{ij}}|H\_1]}{P[T\_{E\_{ij}}|H\_0]} \times \frac{P[T\_{Ci}|H\_1]}{P[T\_{Ci}|H\_0]} \end{split} \tag{8}$$

there are *K* CUAVs, *K<sup>E</sup>* CMs in the *i*th cluster, where *K<sup>E</sup>* = *K* − 1. *Tij* denotes cyclostationary detection statistics or energy detection statistics, *P*[*Tij*|*H*1] and *P*[*Tij*|*H*0] represent the probability density under hypotheses *H*<sup>1</sup> and *H*0, respectively. According to the *TEij* , *TC<sup>i</sup>* , the LRT in Equation (8) can be simplified [36] as

$$T\_i = \sum\_{j=1}^{K\_E} \omega\_{ij} T\_{E\_{ij}} + \rho\_i T\_{\mathbb{C}i}.\tag{9}$$

where *ωij* = *γEij* <sup>2</sup>(1+*γEij*) , and *ρ<sup>i</sup>* = q 2*M*2*γC<sup>i</sup>*

As for the inter-cluster fusion stage, considering the good sensing performance of cyclostationary detection and the large distance between the cluster heads, distributed consensus-based fusion [38] is performed. Each cluster head communicates with its neighbouring cluster heads to exchange information, and the exchange process is iteratively done. The initial information (which is the intra-cluster fusion result) of the *i*th cluster head is denoted as *Ti*(0). Then, according to the network topology, these cluster heads repeatedly iterate until *Ti*(*k*) covers a common value. The consensus-based scheme [38] is

.

$$T\_i(k+1) = \mathcal{W}T\_i(k). \tag{10}$$

where *W* = *I* − *α*∆ <sup>−</sup>1*L*, and *L* is the Laplacian matrix of the cluster heads topology. *α* is the step size, and it satisfies 0 < *α* < *d* −1 , *d* is maximum node degree of the graph. ∆ = *diag*{*δh*<sup>1</sup> , *δh*<sup>2</sup> , . . . , *δh<sup>n</sup>* }, where *δhi* [38] is the weight according to the channel condition of the *i*th cluster head, and it satisfies *δh<sup>i</sup>* ≥ 1. The cluster heads communicate with their own neighbors, then a final consensus [38] is reached as

$$T^\* = \lim\_{k \to \infty} T(k) = \frac{\sum\_{i=1}^n \delta\_{h\_i} T\_i(0)}{\sum\_{h\_{i-1}}^n \delta\_{h\_i}}.\tag{11}$$

#### **4. Continuous Hidden Markov Model Based Spectrum Sensing**

To fully take advantage of the temporal correlation of the spectrum states, we propose a CHMM-based spectrum sensing scheme as shown in Figure 3. Firstly, hidden HMM and CHMM are introduced, and the suitability regarding the spectrum states as a continuous hidden Markov model are analyzed. Then, the model is trained with the forward-back algorithm and Baum–Welch algorithm. The prediction obtained according to CHMM model is used to assist the sensing.

#### *4.1. Hidden Markov Model*

The hidden Markov model is a double stochastic process with a hidden layer and an observable layer. The hidden process is an unobservable Markov chain, which can be obtained through the observed states.

The true states of the spectrum are not observable but the sensing results can be easily obtained. Therefore, the hidden Markov chain model fits well with the PU spectrum state. For HMM, there are three basic problems that need to be solved, i.e., the evaluation problem that computes the probability of the observed fusion result sequence, the learning problem that adjusts the model parameters to maximize the probability of the observed sequence, and the predication problem that calculates the most likely hidden spectrum state sequence according to the observation sequence and model parameters.

In order to avoid the distortion caused by the discretization of continuous variables in the cluster heads, we consider the continuous HMM, which replaces the discrete observation states with continuous characteristics. With more specific spectrum information, we can obtain better detection performance.

#### *4.2. Continuous Hidden Markov Model of Spectrum States*

The PU spectrum state at time instant *t* is given by *x<sup>t</sup>* , and it can be 0 or 1, where 0 represents spectrum absence, 1 denotes spectrum occupancy. The sequence of the PU states *X* = (*x*1, *x*2, . . . , *xt*) can be seen as the hidden Markov chain. *o<sup>t</sup>* is the fusion value of the heterogeneous two-fusion-stage spectrum sensing at time instant *t*, and *O* = (*o*1, *o*2, . . . , *ot*) is the observable layer. The hidden Markov chain and the observable layer constitute a continuous hidden Markov model, which can be formulated as *λ* = (*π*, *A*, *µ*, Σ, *C*) [39], where *π* represents the initial probability vector of the hidden spectrum state, *A* is the transition matrix of the two states. The continuous hidden Markov model can be represented

in Figure 4. *µ*, Σ, *C* are the parameters of the observation probability distribution. The Gaussian mixture model (GMM) is used to model the probability, as the Gaussian process has good adaptability in dealing with complex regression problems and classification. Thus, the observation probability in state *i* according to the GMM can be written as

$$b\_{i}(o) = \sum\_{m=1}^{M} \mathbb{C}\_{im} N(o\_{\prime} \mu\_{im} \Sigma\_{im})\_{\prime} i = 0, 1,\tag{12}$$

which is composed of M Gaussian mixtures. *i* denotes the spectrum state. *Cim* is the proportion of the *m*th mixture coefficient in state *i*, *o* denotes the fusion result calculated by the clustered heterogeneous two-fusion-stage scheme, *µim* and Σ*im* represent the mean and the covariance of *m*th mixture in state *i*, respectively.

**Figure 4.** Continuous hidden Markov model.

To employ the continuous hidden Markov model to cognitive UAV networks, the forward-backward algorithm and Baum–Welch algorithm are utilized to solve the evaluation problem and the training problem, respectively. As for the prediction problem, the Viterbi algorithm is utilized.

#### *4.3. Evaluation Process and Learning Process of Continuous Hidden Markov Model*

For the evaluation problem, the forward-back algorithm is used to calculate the probability of the observed fusion value sequence, and it can be divided into two parts: the forward algorithm and backward algorithm. For a given *λ* and the spectrum state at time *t*, the forward quantity *αt*(*i*) is the joint probability of sequence *O* from the initial time to time *t*, and the state in *S<sup>i</sup>* at time *t*, *βt*(*i*) is the joint probability of sequence from time *t* + 1 to the final time and the state in *S<sup>i</sup>* at time *t*. The forward and backward quantities [39] are defined as

$$\mathfrak{a}\_{\mathfrak{t}}(i) = \mathbf{P}(o\_1, o\_2, \dots, o\_{\mathfrak{t}}, q\_{\mathfrak{t}} = \mathbf{S}\_{\bar{\mathbf{i}}} | \lambda), \tag{13}$$

$$\beta\_l(i) = \mathbb{P}(o\_{t+1}, o\_{t+2}, \dots, o\_T | q\_t = \mathbb{S}\_i, \lambda)\_r \tag{14}$$

with the initializations

$$
\pi\_1(i) = \pi\_i b\_i(o\_1), i = 0, 1,\tag{15}
$$

$$\beta\_T(i) = 1, i = 0, 1. \tag{16}$$

The forward and backward quantities can be calculated by

$$a\_{\mathbf{f}}(j) = [\sum\_{i=1}^{2} a\_{\mathbf{f}}(i) a\_{\mathbf{i},j}] b\_{\mathbf{j}}(o\_{\mathbf{l}}), \mathbf{t} = 1, 2, \dots, T, j = 0, 1,\tag{17}$$

$$\beta\_t(i) = \sum\_{j=1}^{2} a\_{i,j} b\_j(o\_{t+1}) \beta\_{t+1}(j), \\ t = T - 1, \\ T - 2, \dots, 1, \\ i = 0, 1. \tag{18}$$

Then, combining the forward and backward algorithm, we can get the forwardbackward algorithm, then the probability of the observed fusion result sequence *O* with the given model parameters *λ* is obtained [39] as

$$P(O|\lambda) = \sum\_{i=1}^{2} a\_t(i)\beta\_t(i). \tag{19}$$

To solve the learning problem, the Baum–Welch algorithm is adopted, which is one of the expectation–maximization algorithms and uses the forward-backward algorithm in each expectation process. Before using this algorithm, we need to define three parameters: *γt*(*i*), *ξt*(*i*, *j*) and *γt*(*j*, *m*). *γt*(*i*) denotes the probability of the *i*th spectrum state at time *t*, *ξt*(*i*, *j*) denotes the probability that at time *t* the spectrum state is *S<sup>i</sup>* and at time *t* + 1 the spectrum state is *S<sup>j</sup>* , *γt*(*j*, *m*) denotes the probability of the *m*th Gaussian mixture of state *S<sup>j</sup>* at time *t* . *γt*(*i*), *ξt*(*i*, *j*) and *γt*(*j*, *m*) can be calculated [39] as follows:

$$\gamma\_t(i) = \frac{\alpha\_t(i)\beta\_t(i)}{\sum\_{i=1}^2 \alpha\_t(i)\beta\_t(i)},\tag{20}$$

$$\mathfrak{F}\_t(i,j) = \frac{a\_t(i)a\_{i\bar{j}}b\_{\bar{j}}(o\_{t+1})\beta\_{t+1}(j)}{\sum\_{i=1}^2 \sum\_{j=1}^2 a\_t(i)a\_{i\bar{j}}b\_{\bar{j}}(o\_{t+1})\beta\_{t+1}(j)},\tag{21}$$

$$\gamma\_t(j,m) = \frac{\alpha\_t(j)\beta\_t(j)}{\sum\_{i=1}^2 \alpha\_t(i)\beta\_t(i)} \cdot \frac{c\_{j,m}N(o\_{t\prime}\mu\_{j,m\prime}\Sigma\_{j,m})}{\sum\_{n=1}^2 N(o\_{t\prime}\mu\_{j,m\prime}\Sigma\_{j,m})}.\tag{22}$$

When enough training data are provided, that is, the sequences of fusion results obtained from the heterogeneous two-stage-fusion sensing scheme and corresponding spectrum states, the Baum–Welch algorithm offers a way to train the model, outputting good CHMM parameters. In specific, initial model parameters are first selected according to the spectrum condition. Second, *αt*(*i*), *βt*+1(*j*), *γt*(*i*), *ξt*(*i*, *j*) and *γt*(*j*, *m*) are calculated. Third, the parameters are updated according to Appendix A. The forward-backward procedure and the updating procedure are repeated until the probability of the observation sequence *P*(*O*|*λ*) satisfies the convergence condition or the increments of parameters are less than threshold. Finally, we can obtain the trained model parameters *λ* = (*π*, *A*, *µ*, Σ, *C*).

#### *4.4. Predication of Spectrum State with CHMM*

In this section, we adopt the Viterbi algorithm to solve the predication problem. With the learned model parameters and the observed sequence of fusion results, we can calculate the joint probability of the observed sequence. The real spectrum state sequence is calculated [39] as

$$P(O, X|\lambda) = P(O|X, \lambda)P(X, \lambda)$$

$$= \pi\_{X\_1} b\_{X\_1}(o\_1) a\_{X\_1 X\_2} b\_{X\_2}(o\_2) a\_{X\_{T-1} X\_T} b\_{X\_T}(o\_T) \, \tag{23}$$

where *x<sup>t</sup>* denotes the spectrum state at time *t*, *axtxt*+<sup>1</sup> represents the transition probability from state *x<sup>t</sup>* to *xt*+1, *bx<sup>t</sup>* (*ot*) denotes the observation probability of *o<sup>t</sup>* , when the real spectrum state is *x<sup>t</sup>* . Then, select the sequence with the maximum probability as the prediction sequence

$$X = \underset{allX}{\text{argmax}} \, P(O, X|\lambda). \tag{24}$$

Next, we can get the prediction *x<sup>T</sup>* and the prediction sequence *X*, and compare the predication result X with the real state sequence to get the prediction accuracy *P<sup>r</sup>* . When the prediction is "busy", *P<sup>r</sup>* can denote the probability that the spectrum is really occupied at time *T*, 1 − *P<sup>r</sup>* can denote the probability that the spectrum is really absent at time *T*. Similarly, when the prediction is "idle", the probability of occupancy is 1 − *P<sup>r</sup>* , and the probability of absence is *P<sup>r</sup>* .

Combining the prediction accuracy of cluster heads with the fusion result of the detectors, we can get a new false alarm probability and detection probability. There are mainly two kinds of predictions: busy and idle. When the prediction is "busy", we multiply the fusion result *T* by *e* (*e* > 1). Similarly, when the prediction result is "idle", we multiply the fusion value by *η* (*η* < 1). Then we obtain the final fusion statistic adjusted by the prediction. After that, we compare the final fusion statistic with the threshold. When the fusion statistic is larger than the threshold, the decision is "busy" and vice versa. Thus, the final detection probability *P<sup>D</sup>* can be calculated as

$$P\_D = P\_r \cdot P\_d(\epsilon T) + (1 - P\_r) \cdot P\_d(\eta T). \tag{25}$$

#### **5. Analysis of Detection Performance for CHMM-Based Spectrum Sensing**

Under *H*<sup>1</sup> means the spectrum is occupied. When the sampling points are large enough, *TEij* approximately follows normal distribution, as does *TC<sup>j</sup>* (which will be explained in Section 7.1). Thus, when the spectrum is occupied, the fusion result can be calculated as

$$T = \sum\_{i=1}^{n} \delta\_{l\_i} (\sum\_{j=1}^{K\_{i\to}} \omega\_{ij} T\_{\mathcal{E}\_{ij}} + \rho\_i T\_{\mathcal{C}i})$$

$$= \sum\_{i=1}^{n} \sum\_{j=1}^{K\_{i\to}} \delta\_{l\_i} \omega\_{ij} T\_{\mathcal{E}\_{ij}} + \sum\_{i=1}^{n} \delta\_{l\_i l} \rho\_i T\_{\mathcal{C}i} \tag{26}$$

where *n* is the number of clusters, *δhi* is the weight according to the channel condition of the *i*th cluster head, *KiE* is the number of cluster members (CMs) in the *i*th cluster. According to [36], the weight can be simplified as *<sup>ω</sup>ij* <sup>=</sup> *<sup>γ</sup>Eij* 2(1 + *γEij*), and *ρ<sup>i</sup>* = q 2*M*2*γC<sup>i</sup>* . The fusion result *T* is an approximately normally distributed random variable with mean *µ<sup>T</sup>* and variance *σ<sup>T</sup>* 2 , the mean and the variance are expressed as

$$\mu\_T = \sum\_{i=1}^n \sum\_{j=1}^{K\_{\bar{l}\bar{l}}} \delta\_{l\bar{l}\_i} \omega\_{i\bar{j}} (1 + \gamma\_{E\_{\bar{i}\bar{j}}} -^1) + \sum\_{i=1}^n \delta\_{l\bar{l}\_i} \rho\_{\bar{i}} \sqrt{2\gamma\_{C\_{\bar{i}}}} $$

$$= \sum\_{i=1}^n \sum\_{j=1}^{K\_{\bar{l}\bar{l}}} \frac{\delta\_{l\bar{l}\_i}}{2} + \sum\_{i=1}^n 2M \delta\_{l\bar{l}\_i} \gamma\_{C\bar{i}\prime} \tag{27}$$

and

$$\left|\sigma\_{T}\right|^{2} = \sum\_{i=1}^{n} \sum\_{j=1}^{K\_{\rm i\rm E}} \frac{\delta\_{h\_{i}}}{2M(1+\gamma\_{\rm E,j})^{2}} + \sum\_{i=1}^{n} 2M\delta\_{h\_{i}}{}^{2}\gamma\_{\rm Ci}{}^{2}.\tag{28}$$

Then we can obtain the final detection probability *P<sup>D</sup>* of the CUAVNs,

$$P\_D = P\_r \cdot Q(\frac{\frac{\lambda\_T}{\epsilon} - \mu\_T}{\sigma\_T}) + (1 - P\_r) \cdot Q(\frac{\frac{\lambda\_T}{\eta} - \mu\_T}{\sigma\_T}) \tag{29}$$

where *λ<sup>T</sup>* is the threshold.

Under *H*0, when *M* is large enough, *TEij* approximately follows normal distribution, and the probability density *TCij* can be approximately represented as

$$f(T\_{\mathbb{C},i}) = \frac{2}{\sqrt{2\pi}\sigma\_{T\_{\mathbb{C},i}}} \exp\left\{-\frac{{T\_{\mathbb{C},i}}^2}{2\sigma\_{T\_{\mathbb{C},i}}}\right\}.\tag{30}$$

The probability density of the fusion result *T* can be represented as

$$f(T) = \frac{n}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(T-\gamma)^2}{2\sigma^2}\right\} \times \operatorname{erfc}\left(\frac{\sigma\_\mathbb{C}(T-\gamma)}{\sqrt{2\sigma\_\mathbb{E}^2\sigma^2}}\right),\tag{31}$$

where *γ* = ∑ *n <sup>i</sup>*=<sup>1</sup> ∑ *K<sup>E</sup> j*=1 *δhi* {2(1 + *γEij*)}, *σ* <sup>2</sup> = *σ<sup>E</sup>* <sup>2</sup> + *σ<sup>C</sup>* 2 , *σ<sup>C</sup>* <sup>2</sup> = ∑ *n i*=1 2*Mδh<sup>i</sup>* <sup>2</sup>*γC<sup>i</sup>* , *σE* <sup>2</sup> = ∑ *n <sup>i</sup>*=<sup>1</sup> ∑ *K<sup>E</sup> j*=1 *δhi* 2 {2*M*(1 + *γEij*) 2 }. The proof of *f*(*T*) is shown in Appendix B.

The false alarm probability can be calculated as follows, where *F*(*T*) = R *<sup>T</sup>* 0 *f*(*x*)*dx*,

$$\begin{split} P\_F &= P\_r \cdot \left( 1 - F(\eta T) \right) + \left( 1 - P\_r \right) \cdot \left( 1 - F(\varepsilon T) \right) \\ &= 1 - P\_r F(\eta T) - (1 - P\_r) F(\varepsilon T). \end{split} \tag{32}$$

#### **6. CHMM-Based Spectrum Sensing with Practical SNR Estimation**

The above work is based on the assumption that the SNR is perfectly known, i.e., we use a perfectly-known SNR when calculating the fusion weight, detection probability and false alarm probability. In this section, we further consider the scenario that the SNR is not perfectly known, and design a novel SNR estimation algorithm: space smoothing-based M2M4 (SS-M2M4), which modifies M2M4 with spatial smoothness, a technique that is used in the field of signal reconstruction in graphs. This algorithm can provide a more accurate SNR for the proposed CHMM-based spectrum sensing method.

#### *6.1. Typical SNR Estimator*

M2M4 is one of the most widely used blind estimators [14]. M2 and M4 represent the second and the fourth moment of *yn*, respectively, where *y<sup>n</sup>* refers to the samples of the received signal. M2 and M4 can be calculated as follows [14]

$$M\_2 = E\{y\_n y\_n^\*\},\tag{33}$$

$$M\_4 = E\{ (y\_n y\_n^\*)^2 \}.\tag{34}$$

Then with second-order moments and fourth-order moments, we can estimate SNR as follows,

$$\mathfrak{p} = \frac{\sqrt{2M\_2^2 - M\_4}}{M\_2 - \sqrt{2M\_2^2 - M\_4}} \cdot \tag{35}$$

In practice, the second and fourth moments are usually calculated by their own time averages:

$$M\_2 = \frac{1}{M} \sum\_{n=0}^{M} y\_n^2 \,\text{.}\tag{36}$$

$$M\_4 = \frac{1}{M} \sum\_{n=0}^{M} y\_n{}^4. \tag{37}$$

#### *6.2. SS-M2M4 SNR Estimation Algorithm*

In CUAVNs, the CUAVs within the same cluster are close to each other, and thus, their large-scale fading and shadow fading are generally similar [40]. Therefore, their SNRs are correlated. However, the M2M4 algorithm estimates these SNRs individually, ignoring the spatial correlation among the SNRs of those neighbouring users. CUAVNs have their own topology, and the SNR estimation problem can be naturally modeled as graph signal processing problems. Therefore, we propose a novel SNR estimation algorithm as shown in Algorithm 1, which considers the smoothness between neighbors [15]. Firstly, the M2M4 algorithm is applied to calculate the SNR of each CUAV, and the estimation result can be represented as **p**<sup>0</sup> = {*p*01, *p*02, . . . , *p*0*N*}. After that, to ensure that the secondary CUAVs that are close to each other are in the same cluster, the max–min distance clustering algorithm described in Section 3.1 is adopted. Based on the clustering result, we consider each cluster as a new graph, and then correct the original SNR to get the final estimated SNR **p**ˆ = {*p*ˆ1, *p*ˆ2, . . . , *p*ˆ*N*} by Equation (42).

The problem of estimating **pˆ** from the original SNR **p**<sup>0</sup> can be modeled as the following optimization problem,

$$\min\_{\hat{\mathbf{p}}} \quad \frac{1}{2} \| \| \mathbf{p}\_0 - \hat{\mathbf{p}} \|\_{2}^{2} + \frac{\rho}{2} \| \| \mathbf{H} \hat{\mathbf{p}} \|^{2} \tag{38}$$

where **H** is a high-pass graph filter, and *ρ* is the regularization parameter. The first term penalizes the error of the estimated graph signal, the second term encourages the smoothness of the estimated SNR. Similar to [15], we set **H** = **L** 1/2, where **L** is the Laplacian matrix, so that k **Hp**ˆ k 2 <sup>2</sup> = **p**ˆ *<sup>T</sup>***Lp**ˆ. The smoothness of the estimated graph signal can be characterized by the graph Laplacian quadratic form

$$S(\mathfrak{p}) = \mathfrak{p}^T \mathbf{L} \mathfrak{p}.\tag{39}$$

It can be written as

$$\mathcal{S}(\mathfrak{p}) = \sum\_{(i,j \in \varepsilon)} \mathcal{W}\_{ij} (\mathfrak{p}\_j - \mathfrak{p}\_i)^2 \,. \tag{40}$$

where *ε* is the set of edges. The smaller the function value *S*(**p**ˆ) is, the smoother the SNR difference of the cluster, especially when neighboring CUAVs connected by an edge with a large weight have similar values [15]. We can then represent the optimization problem as

$$\min\_{\boldsymbol{\Phi}} \quad \frac{1}{2} \left\| \left. \mathbf{p}\_0 - \boldsymbol{\hat{\mathbf{p}}} \right\| \right\|\_2^2 + \frac{\rho}{2} \boldsymbol{\hat{\mathbf{p}}}^T \mathbf{L} \boldsymbol{\hat{\mathbf{p}}}.\tag{41}$$

To get the best estimated SNRs, we take the derivation of the formula, and get the optimal solution as

$$
\mathfrak{p} = (\mathbf{1} + \rho \mathbf{L})^{-1} \mathbf{p}\_0. \tag{42}
$$

**Algorithm 1** Space smoothing-based M2M4 SNR estimation algorithm

1: **Input:** The received signal sample of each CUAV {*yi*<sup>1</sup> , *yi*<sup>2</sup> , . . . , *yin*}.

2: **Output:** The SNR of all the *N* CUAVs **p**ˆ = {*p*1, *p*2, . . . , *pN*}.

#### **Step 1: M2M4 SNR Estimation**

3: Calculate SNR of all the CUAVs with the traditional M2M4 algorithm in Section 6.1, obtain {*p*01, *p*02, . . . , *p*0*N*}.

## **Step 2: Space Smoothing**

4: Adopt the Max-Min distance clustering algorithm, obtain the clustering result.

5: With the clustering result, calculate the SNR of CUAVs in each cluster by Equation (34), then obtain the final estimated SNR **p**ˆ = {*p*ˆ1, *p*ˆ2, . . . , *p*ˆ*N*}.

#### *6.3. CHMM-Based Spectrum Sensing with SS-M2M4 SNR Estimation*

The estimated SNRs can provide better SNR information for the CHMM-based spectrum sensing procedure, and thus we can obtain more accurate fusion weight in a heterogenous two-stage-fusion stage. Therefore, the CHMM-based spectrum sensing scheme with estimated SNRs achieves a higher detection probability, and the utilization of spectrum sensing can be further enhanced.

The proposed continuous hidden Markov model based spectrum sensing with space smoothing M2M4 SNR estimator can achieve good detection performance in CUAVNs, since it has some advantages over the existing HMM-based spectrum sensing method and the existing SNR estimators. The proposed method can either achieve better performance or be more suitable for UAV applications. The characteristics of these existed methods compared to the proposed method are summarized in Table 2.


**Table 2.** The characteristics of the existing methods.

#### **7. Evaluation and Numerical Results**

In this section, the validity of the approximation of folded normal distribution is firstly presented. Then, the detection performance of the proposed CHMM-based spectrum sensing scheme is evaluated by comparing with the non-CHMM ones. Next, we verify the effectiveness of the proposed SS-M2M4 estimator and further demonstrate the performance of the unified CHMM-based spectrum sensing with the SS-M2M4 estimator. Finally, we further consider mutipath effects in CUAVNs, and verify the effectiveness of the proposed scheme under the Rice channel.

#### *7.1. Approximation of Folded Normal Distribution*

The cluster heads adopt cyclostationary detection to sense the primary spectrum. Under the hypotheses *H*1, the cyclostationary detection statistic *x* follows the folded normal distribution, which can be represented as

$$f(\mathbf{x}) = \frac{1}{\sqrt{2\pi \frac{1}{M}}} \exp\left\{ \frac{(\mathbf{x} - \sqrt{2\gamma})^2}{2\frac{1}{M}} \right\} + \frac{1}{\sqrt{2\pi \frac{1}{M}}} \exp\left\{ \frac{(\mathbf{x} + \sqrt{2\gamma})^2}{2\frac{1}{M}} \right\}, \mathbf{x} > 0. \tag{43}$$

Figure 5 shows the probability density function of the cyclostationary feature under hypotheses *H*<sup>1</sup> when *SNR* = −15 dB, *M* = 2048. As shown in Figure 5, the folded normal distribution mainly coincides with the normal distribution. When *x* is less than 0.0025, which is already 10.8*σ* away from the mean, the two distributions start diverging. In fact, the two terms of folded normal distribution *f*(*x*) can be seen as two normal distributions with opposite means and the same variance. In spectrum sensing, when the channel is occupied, the mean √ 2*γ* is away from 0, and it is much larger than the variance. Thus, when *x* > 0, the second term of *f*(*x*) contributes little to the folded normal distribution. Therefore, in CUAVNs, folded normal distribution can be seen as its first term, that is, a normal distribution in the positive axis.

**Figure 5.** Approximation of folded distribution.

#### *7.2. Performance of CHMM-Based Spectrum Sensing with Perfect SNR Estimation*

In the simulations, the number of secondary CUAVs is set to 20. The CUAVs move according to the random walk mobility model [34], in which the maximum velocity of the CUAVs is 36km/h, and the sensing time is 20 µs. The range of the SNR in our CUAVNs is set as [−15, −3] dB according to [32,33,41]. Locations of secondary nodes lead to different SNRs, we assume that the maximum distance from the transmitter to the secondary UAV is about 1.5 to 2 times the minimum [1]. According to the large-scale fading calculation formula [40], assuming that the fading coefficient is 2, we can obtain 5 dB as a SNR range in our simulation, in other words, the maximum SNR difference of secondary users is within the range of 5 dB. We assume that all secondary CUAVs experience additive white Gaussian noise.

According to [42], the spectrum utilization rate below 3G (The Federal Communications Commission (FCC) of the United States and the European Union (EU) have set 2.4 GHz and 5.8 GHz as the band of civil UAVs. The EU also allocates 433 MHz and 863–870 MHz to UAVs. Similarly, China has set 840.5–845 MHz, 1438–1444 MHz and 2408–2440 MHz as that for UAVs. Compared to 5.8 GHz, more UAVs work on lower bands, below 3 GHz, thus the dilemma of spectrum scarcity is more serious on the below-3G bands) in Berkeley is about 0.3, therefore, we set the initial distribution of the spectrum state distribution as *π* = (0.7, 0.3) *T* , that is to say, the probability of spectrum presence is 0.3, and the probability of spectrum absence is 0.7. We assume a 1st order Markov chain, and the spectrum state at time *t* is known, where the distribution is either (1, 0) *T* or (0, 1) *T* . As stated before, the spectrum utilization is around 0.3, the transition probability of absence

to presence is thus set as 0.25, and in the same way, the transition probability of presence to presence is set as 0.35. Therefore, the transition matrix is *A* = [0.75, 0.25; 0.65, 0.35]. Next, we use MATLAB to generate a spectrum state sequence with a length of 8000 under the parameters above. According to each single spectrum state (hidden state) and different clustering results, the simulated energy detection statistics *TEij* and cyclostationary statistics *TC*,*<sup>i</sup>* are obtained, respectively. Then corresponding observation values *o<sup>t</sup>* are calculated according to the two-step-fusion method.

Figures 6 and 7 show the received operating curve (ROC) of the CHMM-based sensing method and non-CHMM-based ones under the AWGN channel with 20 CUAVs. In Figure 6, the soft–soft represents the heterogeneous two-stage sensing scheme represented in Section 2.2, in which both the intra-cluster fusion stage and the inter-cluster fusion stage adopt a soft combining rule [43]. It can be observed from Figure 6 that the soft–soft heterogeneous sensing scheme with CHMM predication outperforms the non-predication one (soft–soft). Due to CHMM avoiding the distortion caused by the discretization of DHMM, the proposed CHMM-based sensing scheme can further improve the detection probability compared with the DHMM-based soft–soft scheme, which adopts DHMM to model the sensing procedure. What stands out in Figure 6 is the achieved high detection probability of 0.91 when false alarm probability is around 0.1. A higher *P<sup>D</sup>* with small *P<sup>F</sup>* indicates that the proposed algorithm can offer the secondary CUAVs more opportunities to access the spectrum and maintain tolerable interruption to the primary user.

**Figure 6.** ROC of CHMM soft–soft scheme (the proposed CHMM-based heterogeneous two-stage fusion sensing scheme), DHMM soft–soft scheme and non-CHMM soft–soft scheme.

**Figure 7.** ROC of soft–or, or–or scheme and CHMM soft–or, or–or scheme.

In addition, we also simulate the soft–or and or–or schemes [36] to verify the universality of the CHMM model for spectrum sensing in Figure 7, where soft–or means a soft combining rule at the intra-cluster fusion stage and/or a combining rule [43] at the intercluster fusion stage. Similarly, or–or means an or combining rule at both the intra-cluster fusion stage and the inter-cluster fusion stage. As shown in Figure 7, in addition to the good detection performance offered by the heterogeneous soft–soft scheme, the proposed CHMM can also achieve obvious improvement when it is implemented into the other two fusion schemes: or–or and soft–or. Then, with the help of a better sensing performance, the spectrum efficiency and throughput can be further improved.

#### *7.3. Performance of CHMM-Based Spectrum Sensing with SS-M2M4 SNR Estimator*

In this section, the simulation evaluation of our proposed SS-M2M4 in comparison with the existing M2M4 algorithm is presented. The simulation results with SNR in [−15, −3] are shown in Figure 8. Note that the SNR = −15 dB represents the SNR range from [−15, −10]. Similarly, −3 dB represents the range [−3, 2]. *MSE* is used to measure the effect of SNR estimation, and the calculation scheme of *MSE* is given as

$$MSE = \frac{1}{N} \sum\_{i=1}^{N} \left(p\_i - \mathfrak{p}\_i\right)^2,\tag{44}$$

where *p*ˆ*<sup>i</sup>* represents the estimated SNR of the *i*th CUAV, *p<sup>i</sup>* represents the actual SNR of the *i*th CUAV. Figure 8 shows that the proposed scheme offers a good SNR error decrease. It can be seen that Figure 8 also shows that the proposed scheme can improve the performance of SNR estimation up to 3dB (when a 2.3 <sup>×</sup> <sup>10</sup>−<sup>3</sup> *MSE* is required, the M2M4 estimator needs the actual SNR to be −6 dB, but SS-M2M4 only needs it to be −9 dB). An accurate SNR estimator can further help the CUAV get a better detection performance. Consequently, it is reasonable to expect good detection performance when the CHMM-based spectrum sensing adopts the novel SNR estimation algorithm.

**Figure 8.** MSE of SS-M2M4 estimator and M2M4 estimator.

Figure 9 shows the CHMM-based spectrum sensing with the SS-M2M4 estimator significantly outperforming the original ones, which either do not use the CHMM model (spectrum sensing + SS-M2M4 SNR estimator) or employ the M2M4 estimator (CHMM spectrum sensing + M2M4 SNR estimator), or neither (spectrum sensing + M2M4 SNR estimator). Specifically, when the false alarm is 0.1, the unified scheme achieves a detection probability of 0.95, while the CHMM-based spectrum sensing with the original SNR estimator (CHMM spectrum sensing + M2M4 SNR estimator) can only achieve around 0.82. Besides, the unified scheme can further enhance the detection performance compared with the DHMM-based spectrum sensing with the proposed SS-M2M4 (DHMM spectrum

sensing + SS-M2M4 SNR estimator), and the DHMM-based spectrum sensing with the proposed SNR estimator can also outperform the DHMM-based spectrum sensing with the original SNR estimator (DHMM spectrum sensing + M2M4 SNR estimator). In other words, the proposed SNR estimator has a universality to a different spectrum sensing method and the unified scheme can utilize the spectrum more efficiently and maintain a tiny interference to the primary user.

**Figure 9.** Detection probability versus false alarm probability.

#### *7.4. Performance of CHMM-Based Spectrum Sensing with SS-M2M4 SNR Estimator under Rice Channel*

In the above simulation, we assume the AWGN channel and verify the effectiveness of the proposed scheme under the AWGN channel. In this section, we have further considered the Rice channel when mutipath effects have been introduced in UAV applications [44]. The simulation of the proposed scheme under the Rice channel is shown in Figure 10, according to [44], we set the Rician factor as K = 10. Figure 10 shows that even under the Rice channel, the proposed method (Rice CHMM sensing+SS-M2M4 estimator) can obtain better performance compared with DHMM-based spectrum sensing with the proposed SS-M2M4 SNR estimator (Rice DHMM sensing + SS-M2M4 estimator), CHMM-based spectrum sensing with M2M4 SNR estimator (Rice CHMM sensing+M2M4 estimator), and other methods. Although the detection probability under the Rice channel is lower than that under the AWGN channel, the proposed scheme can also achieve better performance than other schemes under the Rice channel, in other words, our scheme is effective in CUANs.

**Figure 10.** Detection probability versus false alarm probability under the Rice channel.

#### **8. Conclusions**

In this paper, we consider modeling primary user states as a Markov chain, and propose a spectrum sensing scheme based on a continuous hidden Markov model with perfect SNR estimation. We derive the detection probability and false alarm probability of the heterogenous-fusion clustered CUAVNs. Taking the similarity of the SNRs of CUAVs in the same cluster into account, we propose a space smoothing based SNR estimator for sensing in CUAVNs to offer a more accurate SNR to the proposed sensing method. Simulation results show that the unified CHMM-based sensing scheme with the proposed SNR estimator enhances the sensing performance considerably.

The work can be further extended in the following aspects in our future research. First, in the proposed CHMM based spectrum sensing method, the hidden spectrum chain is modeled by first-order Markov chain, and it has not made full use of the existing historical information. In the future research, the high-order Markov chain can be used to model the spectrum sensing process to further improve the accuracy of prediction. Second, in this paper, we considered cognitive UAV networks with a single primary user. Recently, cognitive UAV networks with multi primary users [33] have been proposed to enhance the spectrum utilization. Thus, we can explore our work in multi-PUs CUAVNs in the future.

**Author Contributions:** Conceptualization and methodology, Y.F. and W.X.; software and validation, Y.F.; formal analysis, Y.F.; investigation, Y.F., W.X. and Z.Z.; resources, Y.F.; data curation, Y.F.; writing original draft preparation, Y.F.; writing—review and editing, Y.F. and F.W.; visualization,Y.F. and F.W.; supervision, W.X.; project administration, Z.Z.; funding acquisition, W.X. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported in part by the National Key Research and Development Program of China under Grant 2019YFC1511302; and in part by the National Natural Science Foundation of China under Grant 61871057.

**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.

#### **Appendix A. Formulas of the Updated CHMM Parameters**

The parameters of the CHMM [39] *λ* = (*π*, *A*, *µ*, Σ, *C*) can be updated with parameters *γt*(*i*), *ξt*(*i*, *j*) and *γt*(*j*, *m*) as follows

$$
\pi\_{\mathbf{i}} = \gamma\_{\mathbf{i}}(\mathbf{i}),\tag{A1}
$$

$$a\_{ij} = \frac{\sum\_{t=1}^{T-1} \delta\_t(i, j)}{\sum\_{t=1}^{T-1} \gamma\_t(i)},\tag{A2}$$

$$c\_{j,m} = \frac{\sum\_{t=1}^{T} \gamma\_t(j,m)}{\sum\_{t=1}^{T} \sum\_{n=1}^{2} \gamma\_t(j,n)},\tag{A3}$$

$$\mu\_{j,m} = \frac{\sum\_{t=1}^{T} \gamma\_t(j\_\prime m) o\_t}{\sum\_{t=1}^{T} \gamma\_t(j\_\prime m)},\tag{A4}$$

$$\Sigma\_{j,m} = \frac{\sum\_{t=1}^{T} \gamma\_t(j\_\prime m)(o\_t - \mu\_{j,m})(o\_t - \mu\_{j,m})^T}{\sum\_{t=1}^{T} \gamma\_t(j\_\prime m)}. \tag{A5}$$

#### **Appendix B. Proof of Distribution of the Fusion Value**

In this section, we derive the distribution of the fusion result *T*, where *T* = *T<sup>E</sup>* + *TC*. To obtain the fusion result *T*, we first add the energy detection result *TEij* and the cyclostationary detection result *TC<sup>i</sup>* , respectively, and then calculate the summation of the

two different detection results. When *M* is large enough, the sum of the energy detection *T<sup>E</sup>* follows normal distribution written as

$$f(T\_E) = \frac{1}{\sqrt{2\pi}\sigma\_E^2} e^{-\frac{(T\_E-\gamma)^2}{2\sigma\_E^2}}\text{.}\tag{A6}$$

where mean *γ* = ∑ *n <sup>i</sup>*=<sup>1</sup> ∑ *K<sup>E</sup> j*=1 *δhi* /{2(1 + *γEij*)}, variance *σ<sup>E</sup>* <sup>2</sup> = ∑ *n <sup>i</sup>*=<sup>1</sup> ∑ *K<sup>E</sup> j*=1 *δhi* 2 {2*M*(1 + *γEij*) 2 }. To obtain the distribution which the sum of cyclostationary detection follows, we do the following derivation. For simplicity, x denotes one cyclostationary detection result *TC<sup>i</sup>* , y denotes another cyclostationary detection result *TC<sup>j</sup>* , and the distribution of *x* and *y* can be represent as

$$f(\mathbf{x}) = \frac{2}{\sqrt{2\pi}\sigma} e^{-\frac{\mathbf{x}^2}{2\sigma^2}},\tag{A7}$$

$$f(y) = \frac{2}{\sqrt{2\pi}\sigma} e^{-\frac{y^2}{2\sigma^2}},\tag{A8}$$

with the same mean 0 and the same variance *σ* <sup>2</sup> = <sup>2</sup> *<sup>M</sup>* , *z* = *ax* + *by* is defined as the sum of x and y, we can derive the distribution of *z* as

$$\begin{split} f(z) &= \int\_0^{+\infty} \frac{1}{a} f\_x \left( \frac{z - by}{a} \right) f\_y(y) dy \\ &= \int\_0^{+\infty} \frac{1}{a} \frac{2}{\sqrt{2\pi}\sigma} \exp\left\{ - \frac{(\frac{z - by}{a})^2}{2\sigma^2} \right\} \frac{2}{\sqrt{2\pi}\sigma} \exp\left\{ - \frac{y^2}{2\sigma^2} \right\} dy \\ &= \left(\frac{2}{\sqrt{2a\pi\sigma}}\right)^2 \exp\left\{ - \frac{\frac{z^2}{a^2 + b^2}}{2\sigma^2} \right\} \int\_0^{+\infty} \exp\left\{ - \frac{(\sqrt{a^2 + b^2}y - \frac{bz}{\sqrt{a^2 + b^2}})^2}{2a^2\sigma^2} \right\} dy \\ &= \left(\frac{2}{\sqrt{2a\pi\sigma}}\right)^2 \int\_{-\frac{b}{\sqrt{2a\pi}}}^{+\infty} \exp\left\{ -t^2 \right\} dt \frac{\sqrt{2a}\sigma}{\sqrt{a^2 + b^2}} \exp\left\{ - \frac{z^2}{2(a^2 + b^2)\sigma^2} \right\} \\ &= \left(\frac{2}{\sqrt{2a\pi\sigma}}\right)^2 \frac{\sqrt{2a}\sigma^2}{\sqrt{a^2 + b^2}} \exp\left\{ - \frac{z^2}{2(a^2 + b^2)\sigma^2} \right\} \int\_{-\infty}^{+\infty} \exp\left\{ - t^2 \right\} dt \\ &= \frac{4}{\sqrt{2\pi(a^2 + b^2)\sigma^2}} \exp\left\{ - \frac{z^2}{2(a^2 + b^2)\sigma^2} \right\}, \end{split} \tag{A9}$$

where *t* = √ *a* <sup>2</sup>+*b* <sup>2</sup>*y*<sup>−</sup> <sup>√</sup> *bz a* 2+*b* <sup>2</sup> <sup>√</sup> 2*aσ* , *y* = √ <sup>2</sup>*aσt*+ √ *bz a* 2+*b* <sup>2</sup> <sup>√</sup> *a* <sup>2</sup>+*b* 2 . From the above proof, we can conclude that the distribution of *z* is similar with *x* and *y*, but there exist some differences: the coefficient has doubled, the variance is the sum of *σ<sup>x</sup>* and *σy*. Thus *T<sup>C</sup>* follows the following distribution

$$f(T\_{\mathbb{C}}) = \frac{2n}{\sqrt{2\pi\sigma\_{\mathbb{C}}^2}} \exp\left\{-\frac{T\_{\mathbb{C}}}{2\sigma\_{\mathbb{C}}}^2\right\}.\tag{A10}$$

After obtaining the distribution of cyclostationary detection, we can further derive the distribution of the fusion result *T* as follows,

$$\begin{split} f(T) &= \int\_{0}^{+\infty} f\_{\rm T\_{\rm C}}(T - T\_{\rm C}) f\_{\rm T\_{\rm C}}(T\_{\rm C}) dT\_{\rm C} \\ &= \frac{2n}{\sqrt{2\pi\sigma\_{\rm E}^{2}}\sqrt{2\pi\sigma\_{\rm C}^{2}}} \int\_{0}^{+\infty} \exp\left\{-\frac{T\_{\rm C}^{-2}}{2\sigma\_{\rm C}^{2}}\right\} \exp\{-\frac{\left(T - T\_{\rm C} - \gamma\right)^{2}}{2\sigma\_{\rm E}^{2}}\} dT\_{\rm C} \\ &= \frac{2n}{2\pi\sqrt{\sigma\_{\rm E}^{2}\sigma\_{\rm C}^{2}}} \exp\left\{-\frac{\left(T - \gamma\right)^{2}}{2\sigma^{2}}\right\} \int\_{0}^{+\infty} \exp\left\{\frac{-\left[\sigma T\_{\rm C} - \frac{\sigma\sigma^{2}\left(T - \gamma\right)}{\sigma}\right]^{2}}{2\sigma\_{\rm E}^{2}\sigma\_{\rm C}^{2}}\right} \right\} dT\_{\rm C} \\ &= \frac{2n}{2\pi\sqrt{\sigma\_{\rm E}^{2}\sigma\_{\rm C}^{2}}} \exp\left\{-\frac{\left(T - \gamma\right)^{2}}{2\sigma^{2}}\right\} \frac{2\sqrt{\sigma\_{\rm E}^{2}\sigma\_{\rm C}^{2}}}{\sigma} \int\_{0}^{+\infty} \exp\left\{-t^{2}\right\} dt \\ &= \frac{n}{\sqrt{2\pi}\sigma^{2}} \exp\left\{-\frac{\left(T - \gamma\right)^{2}}{2\sigma^{2}}\right\} \text{erfc}\left(\frac{\sigma\_{\rm C}(T-\gamma)}{\sqrt{2\sigma\_{\rm E}^{2}\sigma^{2}}}\right). \end{split} \tag{A11}$$

where *n* is number of the clusters, *t* = p *σE* <sup>2</sup> + *σ<sup>C</sup>* <sup>2</sup>*T<sup>C</sup>* − *σC* 2 √ (*T*−*γ*) *σE* <sup>2</sup>+*σ<sup>C</sup>* 2 / p 2*σ<sup>E</sup>* <sup>2</sup>*σ<sup>C</sup>* 2 , *σ<sup>C</sup>* <sup>2</sup> = ∑ *n* 2*Mδh<sup>i</sup>* <sup>2</sup>*γC<sup>i</sup>* , *σ* <sup>2</sup> = *σ<sup>E</sup>* <sup>2</sup> + *σ<sup>C</sup>* 2 . Finally, the distribution of the fusion result *T* under

*i*=1 hypotheses *H*<sup>0</sup> can be obtained as

$$f(T) = \frac{n}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(T-\gamma)^2}{2\sigma^2}\right\} \text{erfc}\left(\frac{\sigma\_C(T-\gamma)}{\sqrt{2\sigma\_E^2\sigma^2}}\right). \tag{A12}$$

#### **References**

