Spectrum Situation Awareness for Space–Air–Ground Integrated Networks Based on Tensor Computing

The spectrum situation awareness problem in space–air–ground integrated networks (SAGINs) is studied from a tensor-computing perspective. Tensor and tensor computing, including tensor decomposition, tensor completion and tensor eigenvalues, can satisfy the application requirements of SAGINs. Tensors can effectively handle multidimensional heterogeneous big data generated by SAGINs. Tensor computing is used to process the big data, with tensor decomposition being used for dimensionality reduction to reduce storage space, and tensor completion utilized for numeric supplementation to overcome the missing data problem. Notably, tensor eigenvalues are used to indicate the intrinsic correlations within the big data. A tensor data model is designed for space–air–ground integrated networks from multiple dimensions. Based on the multidimensional tensor data model, a novel tensor-computing-based spectrum situation awareness scheme is proposed. Two tensor eigenvalue calculation algorithms are studied to generate tensor eigenvalues. The distribution characteristics of tensor eigenvalues are used to design spectrum sensing schemes with hypothesis tests. The main advantage of this algorithm based on tensor eigenvalue distributions is that the statistics of spectrum situation awareness can be completely characterized by tensor eigenvalues. The feasibility of spectrum situation awareness based on tensor eigenvalues is evaluated by simulation results. The new application paradigm of tensor eigenvalue provides a novel direction for practical applications of tensor theory.


Introduction 1.Space-Air-Ground Integrated Networks and Situation Awareness
The space-air-ground integrated networks (SAGINs), working as network infrastructures, provide ubiquitous, collaborative, and efficient information services for various network applications in large-scale space.The SAGINs are composed of three kinds of networks, including space networks, air networks, and ground networks.The ground networks are the main body of SAGINs.The space networks and air networks serve primarily as supplements and extensions.Through the deep integration of multidimensional networks, the SAGINss can effectively utilize various resources comprehensively, carry out intelligent network control and information processing so as to flexibly cope with the network services with different demands, and realize the functions of an integrated communication system, which is called communication-computing-cache.
Under the framework of SAGINs, the space networks consist of various satellite systems that form the space-backbone network and space-access network, providing global coverage, ubiquitous connectivity, broadband access, and communication services, mainly for rural areas and the regions where ground networks are difficult to be deployed.The space networks are composed of high-altitude communication platforms and UAV networks, which can enhance coverage, enable edge network access, and provide additional network services in emergency situations.The ground networks are mainly composed of ground internet and mobile communication networks, which are responsible for the network service in the business-intensive areas like cities or areas requiring high communication quality.

Space Networks
The space networks consist of different types of satellites, constellations, and the corresponding ground infrastructures, including the ground stations and the control centers.Satellites can be divided into GEO, MEO, LEO [1], and VLEO [2], according to the altitudes where they are working.GEO satellites, also called high elliptical orbit satellites, work at an altitude of more than 36,000 km and are often used for international long-distance communication.MEO satellites orbit at an altitude of 2000-36,000 km, and LEO satellites orbit at an altitude of 400-2000 km.Because their orbits are relatively low, the former are often used as positioning systems, such as GPS and Beidou navigation satellite systems, while the latter are often used for earth observation, earth surveys, and space stations.VLEO satellites have the lowest orbits, so they can be used to provide high-speed data communication, precise positioning, and other services.
Since satellites have a high position and wide coverage, they are regarded as the main method to enable global mobile communications.Iridium, for example, expects to provide global voice and data communication by adopting VLEO satellites.Starlink creates a low-cost, high-coverage space network that supports global communications and can provide internet access services at the speeds of optical fiber by increasing the density of satellites.Meanwhile, due to the high position of the satellites, transmission links are long and the propagation delays are large, making it easy to be destroyed by random or deliberate nodes attacking of malicious nodes during signal transmission.So, it is difficult to guarantee the QoS of real-time information interactive applications.

Air Networks
The air networks consist of stratospheric airships, high-altitude balloons, drones, helicopters, and other high-altitude platforms.As an intermediate layer for SAGINs, it can provide data routing between the platforms and the ground networks and exchange data information with space networks.High-altitude platforms (HAPs) are the main body of the air networks, which are located 2-20 km above the stratosphere and upper troposphere [3].HAPs, compared with satellite communication platforms, have the advantages of simple deployment, short communication response time, and low cost and are often used in temporary high-bandwidth communication scenarios, such as emergency communication and disaster relief activities.
However, due to the rapid changes in HAPs network topology and communication links, in order to ensure the security and stability of the network and provide users with high throughput and low delay network access services, it is necessary to design effective coordination mechanisms such as channel allocation, which may result in a complex network structure or allocation mechanism and high allocation cost.

Ground Networks
The ground networks are composed of many sub-networks, such as Wireless Local Area Networks (WLAN), cellular networks, and Ad Hoc Networks.Most of the common communication modes in our daily lives can be classified into the ground networks.Meanwhile, the ground networks can provide users with communication services with high data transmission rates and throughputs.
However, as the construction of base stations and other infrastructure is more limited by complex terrain, the high construction cost leads to poor network services provided by the ground networks in remote areas with sparse populations.In addition, ground infrastructure is vulnerable to extremely severe weather, man-made damage, and other factors, resulting in communication disruption.To sum up, it is difficult to depend on the ground networks singly to satisfy the rapidly increasing and complicating communication needs.

Spectrum Situation Awareness
Spectrum situational awareness originates from spectrum sensing technology and the concept of situational awareness.Spectrum sensing was first proposed in [4], where Mitola coined a new concept of cognitive radio (CR) to realize dynamic spectrum access.In CR networks, the secondary users (SUs) can access the idle channels that are not occupied by the primary users (PUs) to totally improve spectrum efficiency.Spectrum sensing technology can obtain spectrum usage information in wireless systems through various signal detection methods, detect spectrum holes, and prevent interference to the PUs.The concept of situational awareness originated in the 1980s.It is mainly used to analyze information about the environment to acquire current and future situations and make corresponding judgments and decisions.Now, it is extended to perceive the environmental elements in a certain time and a certain space, understand the meaning of these elements and predict their future state.
There is still no exact definition of spectrum situational awareness at the moment.However, according to the core purpose of spectrum sensing and situational awareness, it can be understood as the acquisition of various state information from the current spectrum space, such as spectrum busy state and spectrum radiated power.On this basis, various parameters and the development trend of spectrum space can be analyzed.Its core technologies are summarized as wide-area spectrum situational awareness, dynamic spectrum situation generation, and spectrum situation utilization [5].

Tensor Computing and Tensor Eigenvalues
Tensor computing is a new concept mainly including tensor decomposition, tensor completion, and tensor eigenvalue.Tensor decomposition and completion have been widely used in signal processing [6], machine learning [7], big data analysis [8] and other fields [9,10], while engineering applications of tensor eigenvalue in the field of communications are relatively lacking.

Tensor Decomposition and Tensor Completion
Tensor decomposition originated with Hitchcock in 1927 [11], and the conception of the multiway model was proposed by Cattell in 1944 [12].It can be viewed as the higher-order generalization of matrix factorization, that is, converting higher-order data into a combination of lower-order and lower-dimension data.The two most common types of tensor decomposition are CP decomposition and Tucker decomposition.
Based on the definition of the rank-one tensor, Hitchcock proposed that the tensor was divided into a finite number of rank-one tensors for the first time, named as polyadic form, which was the rudiments of CP decomposition.The name of CP decomposition has not been clearly given; parallel factors [13] and CANDECOMP [14] both were the former names of CP decomposition until Kiers called this form of tensor representation method CANDECOMP/PARAFAC: that is, CP decomposition [15].For example, given a third-order tensor X ∈ R I×J×K , it can be re-written by CP decomposition as where the positive integer R is called CP-rank, and a r ∈ R I , b r ∈ R J , and c r ∈ R K for r = 1, . . ., R; the notation "•" represents the outer product.Tucker decomposition, first proposed by Tucker in 1963 [16], can be regarded as a higher-order extension of principal component analysis (PCA).Like CP decomposition, Tucker decomposition also has many other names, such as three-mode PCA (3MPCA) [17] and higher-order SVD (HOSVD) [18].For example, for a third-order tensor X ∈ R I×J×K , we can rewrite it by Tucker decomposition as where A ∈ R I×L , B ∈ R J×M , and C ∈ R K×N are factor matrices.The tensor G ∈ R L×M×N is called the core tensor of Tucker decomposition, and a l ∈ R I , b m ∈ R J and c n ∈ R K for l = 1, . . ., L, m = 1, . . ., M, n = 1, . . ., N. The notation "× n " represents the n-mode product.
Tensor completion is used to complete the tensor by estimating the values according to the relationship between the existing data, structural properties of data, and the missing elements, which is often used in pattern recognition [19], compressed sensing [20] and other fields.Tensor completion is mainly divided into two kinds of methods [21].One is based on the low-rank property, called the nuclear norm minimization method, like [22,23].The other is based on low-rank tensor decomposition, like [24,25].

Tensor Eigenvalue
Tensor eigenvalue is the essential part of tensor computing, which is developed from the concept of matrix eigenvalue.Unlike a singular value or eigenvalue after the matricization of a tensor, which is a matrix singular value or eigenvalue, respectively, here, we discuss the tensor eigenvalue, which is calculated by regarding the tensor as a whole unit.For different research purposes and application backgrounds, scholars have put forward many different concepts of tensor eigenvalue, but the origin of tensor eigenvalue was proposed by professors Qi and Lim based on their respective research directions independently.The former treats the tensor eigenvalues as roots of multidimensional polynomials [26], while the latter proposed the concept of tensor eigenvalue by analogy with the Rayleigh quotient of symmetric matrix eigenvalues and the constraint variational method [27].Although their starting points are different, the two concepts are essentially the same.Since the introduction of tensor eigenvalues, with a particular focus on Z-eigenvalues, they have been extensively utilized in a wide range of fields.These applications include global optimization problems [28,29], hypergraph theory [30,31], homogeneous polynomial system stability problems [32], and many others.
For supersymmetric tensors, where tensor elements have the same value if their indexes belong to the full permutation set of indexes, Qi called a number λ an N-eigenvalue of A ∈ R [m,n] if λ is the solution of the following homogeneous polynomial equation where , and he called the solution x an N-eigenvector of A associated with the N-eigenvalue λ.If λ ∈ R, x ∈ R n , λ is called an H-eigenvalue and x is called an H-eigenvector.Lim called them l m -eigenvalue and l m -eigenvector, which were defined as follows where I n ∈ R n is an all-one vector.E-eigenvalues and Z-eigenvalues satisfy In order to solve the problem of eigenvalue computing of non-symmetric tensors and different structure tensors as well as summarize different types of tensor eigenvalues, B-eigenpairs and B R -eigenpairs [33,34], where they are called B R -eigenpairs if B-eigenpairs are all real, are proposed and defined as Therein, B ∈ R [m ′ ,n] is a symmetric tensor, which has different forms depending on the different types of tensor eigenvalues we want to express.The notation "A {k} " means to find the eigenvalues of A in the k-th direction.Because A is a non-symmetric tensor, different order directions lead to different eigenpairs, and k is needed to indicate the order direction and the B-eigenpair can also be called a mode-k B-eigenpair.
Four types of eigenpairs described above can be represented by B-eigenpair and B R -eigenpair.If B is the unit tensor I ∈ R [m,n] , and m = m ′ , the mode-1 B-eigenpairs are the N-eigenpairs and the mode-1 B R -eigenpairs are the H-eigenpairs.If B is the unit matrix I ∈ R n×n , and m ′ = 2, the mode-1 B-eigenpairs are the E-eigenpairs and the mode-1 B R -eigenpairs are the Z-eigenpairs.
The main contributions of the paper can be summarized as follows: 1.
A novel tensor-based spectrum situational awareness model is proposed to store and process multidimensional, heterogeneous, and massive spectral big data from space-air-ground integrated networks.

2.
Two eigenvalue computing schemes, including the semidefinite relaxation algorithm and the homotopy continuation algorithm, are included to calculate the eigenvalues of spectrum situational awareness tensors.The computing performances of two algorithms are evaluated, and the superiority of the homotopy algorithm for tensor eigenvalue estimation is demonstrated.

3.
A novel spectrum situational awareness scheme based on tensor eigenvalues is proposed, where the tensor eigenvalue distribution is used to evaluate the state of the spectrum.The new application paradigm of a tensor eigenvalue provides a novel direction for practical applications of tensor eigenvalues, especially using tensor eigenvalue distributions to construct a spectrum situational awareness scheme.
The remainder of the paper is provided as follows.Section 2 introduces the spaceair-ground integrated network model and the tensor-based spectrum situation awareness model.Section 3 introduces the problem of tensor eigenvalue computing and two classical tensor eigenvalue computing algorithms, a semidefinite relaxation algorithm and a homotopy continuation algorithm.The performance of the two algorithms are compared in terms of convergence ratio, accuracy and CPU time.Section 4 introduces a tensor-eigenvaluebased spectrum sensing algorithm and simulation results with theoretical analysis.In Section 5, the paper is concluded.

Space-Air-Ground Integrated Network Model
The model of SAGINs is shown in Figure 1, where space networks, air networks and ground networks are three layers.As shown in the figure, the satellites including GEO, MEO, LEO and VLEO exchange information and realize key communication functions through intra-segment wireless links.In air networks, high-altitude platforms, such as aircrafts, airships, balloons, and UAVs, constitute corresponding sublayer networks.The air networks can provide special functions for emergency communications and other tasks.In ground networks, 5G networks can provide basic network infrastructures and basic internet access services.Many services of ground networks are implemented through intersegment communication, like using air networks for Internet access enhancement services in emergency communications and using space networks for global communications through ground gateways.
In this model, sublayers within each network are independent of each other so as to give full play to their characteristics and advantages and facilitate the modification of the topology structure.Inter-segment networks are interwoven and fused with each other to realize the integration of space networks, air networks, and ground networks to complete the core purpose, accurate acquisition, rapid processing, and efficient transmission of information.A whole model and unified heterogeneous big data model are required!

Tensor-Based Spectrum Situation Awareness Model
For the proposed unified SAGINs model, the corresponding big data model is required to store and process heterogeneous big data.We try to propose a spectrum situational awareness model, adopt the unified data tensor model to solve the above problem, and use the data of this model to deal with the spectrum utilization issues from the perspective of spectrum usage.
For multidimensional, heterogeneous, and massive spectral big data, each order of tensors is used to represent a specific dimension, such as time, space, and frequency.The smallest unit of the tensor is constant; that is, the element of the tensor whose index is greater than or equal to the order is constant.In practical application, it is often necessary to collect multiple data in one location.Since these data do not belong to the same variable, they are often placed separately in multiple tensors of the same size during tensor modeling; that is, tensors are added to these data, respectively.Different tensors are represented by different notations, resulting in heterogeneous data.The huge amount of symbols and the heterogeneous data tensors pose a challenge to the final form of the same representation of big data.
To overcome this challenge, the tensor representation model of vector elements is proposed.Using vectors instead of scalars as the smallest unit of the tensor, all information can be expressed without error as long as the vector elements are ordered in the same order.The tensor representation model of vector elements is represented by the symbol A i 1 i 2 ...i m (j) , where m is the order of tensor A. For example, the temperature and humidity in a certain space are measured, and three sampling points are taken for each length, width, and height.The tensor model is a three-order three-dimensional tensor A with two variables, and the symbol is denoted as A xyz(a) , where A xyz = (t, h) and x, y, z, t, and h represent the length, width, height, temperature, and humidity, respectively.a is the sequence number of the variable in the vector.A xyz(1) represents the temperature tensor of the space and A 111(1) represents the specific temperature value at sampling point 111.In other words, the uniform variable is formed as the order of the tensor, and the remaining variables are converted into vectors and stored in the tensor.
Using the tensor representation model of vector elements above, we can obtain the tensor-based spectrum situation awareness model.For the space-air-ground integrated spectrum situational awareness, user-centered situational awareness is carried out, and the awareness model is a five-order tensor composed of time, frequency, longitude, latitude, and altitude, where the elements of the tensor are data vectors.A five-order tensor A ∈ R T×F×X×Y×Z , composed of time, frequency, longitude, latitude, and altitude, is taken as an example of the awareness model here, and the model is shown in Figure 2

Related Work of Tensor Eigenvalue Calculation
Unlike matrix eigenvalue generation, the eigenvalue calculation problem of third or higher-order tensors can be considered as an NP-hard problem [35] due to the socalled "curse of dimensionality".Nonetheless, several algorithms computing one or some eigenvalues of a tensor have been developed recently.Most algorithms are flawed, and these algorithms are designed for tensors of specific types, such as non-negative or real symmetric tensors.
For non-negative tensors, Ng, Qi, and Zhou proposed a power-type method to compute the largest H-eigenvalue of a non-negative tensor, which was called the Ni-Qi-Zhou method based on the Perron-Frobenius theorem [36].For real symmetric tensors, Hu, Huang, and Qi proposed a sequential semidefinite programming method for computing extreme Z-eigenvalues [37].Hao proposed a sequential subspace projection method for a similar purpose [38].Kolda and Mayo proposed a shifted power method (SSHOPM) for computing Z-eigenvalues [39].Han proposed an unconstrained optimization method for computing a real generalized eigenvalue for every order real symmetric tensor [40].Lv and Ma proposed a Levenberg-Marquardt method to obtain all the H-eigenvalues of real semi-symmetric tensors [41].In addition, for all symmetric tensor eigenvalues, there are many tensor eigenvalue algorithms such as NCM (Newton correction method), O-NCM [42] and FNS (fast Newton-Shultz-type iterative method) [43].
In recent years, the research on the eigenvalue computing of general tensors has made a breakthrough.For general tensors, Cui, Dai, and Nie proposed a novel method for computing all real eigenvalues of symmetric tensors by semidefinite relaxation [33] and then extended it to general tensors [44].Chen, Han, and Zhou also proposed a homotopy method for computing all eigenvalues [34].For solving tensor equations with applications, Liang, Zheng, and Zhao proposed alternating iterative methods based on ADMM, such as G-ADMM (Gauss-Seidel scheme) and TT-ADMM (tensor-train) [45].Chen et al. provided a new idea to compute tensor eigenvalues by using digital signal processing technology and proved its feasibility from the perspective of a continuous-mode system [46].
To sum up, since the concept of the tensor eigenvalue was proposed, scholars have paid much attention to solving the tensor eigenvalue problem and put forward a mass of tensor eigenvalue numerical approximate algorithms, but most of them are based on the Newton iteration method to estimate values, just changing the initial conditions or the iteration equation to accelerate the rate of equation convergence, and a few use convex optimization algorithms to solve eigenvalue problems with special structures.

Semidefinite Relaxation Algorithm
The semidefinite relaxation method, which is one of the first few algorithms to try to compute all eigenvalues of a tensor, is of great significance for the proposal and improvement of the subsequent algorithm, although it has some deficiencies like low convergence rate, long convergence time, inability to compute complex eigenvalues, and others.For the semidefinite relaxation method, to calculate Z-eigenvalues or E-eigenvalues, the corresponding algorithm is different, but the main idea is identical.Below, we take its Z-eigenvalue algorithm as an example to introduce the semidefinite relaxation method.
Let a eigenpair (λ, x) be a Z-eigenpair of A ∈ R [m,n] .We can derive that Z-eigenvalue λ is Ax m and the corresponding eigenvector x needs to satisfy Ax m−1 = (Ax m )x, x T x = 1.So, we can obtain a polynomial function Then, x is a Z-eigenvector of A when p z (x) = 0.The Z-eigenvalues can be calculated from the smallest to the largest.
Firstly, compute the smallest Z-eigenvalue λ 1 of A. Define the polynomial optimization problem min f (x where the optimal solution of ( 8) is λ 1 .According to Lasserre's hierarchy [47], using the semidefinite relaxation method, the polynomial optimization problem (8) can be converted to where ⟨•⟩, y, L (k) p (y), and M k (y) are defined in [47].Then, with k = k 0 as the initial point, the solutions of the semidefinite relaxation (9) can be obtained, where k 0 = ⌈(m + 1)/2⌉.If there is no solution for k = k 0 , A has no Z-eigenvalue; otherwise, solve (9) again with an optimizer y * .If y * satisfies rank M t−k 0 (y * ) = rank M t (y * ), we obtain λ 1 = f k 1 ; otherwise, let k = k + 1 and repeat the above procedures.In the end, we can obtain the smallest Z-eigenvalue λ 1 .
Secondly, we need to know whether the next larger Z-eigenvalue λ i+1 exists or not and then compute λ i+1 if it exists.Let δ be a positive number that is close to zero.To find larger eigenvalues, (9) is modified to the following formula Like (9), we can obtain Lasserre's hierarchy of semidefinite relaxations If the semidefinite relaxation (11) converges, which is checked by condition (10), the Z-eigenvectors x 1 , . . ., x r can be computed by the method in [48].Then, consider the optimization problems v + (λ i , δ) := max f (x) and v − (λ i , δ) := min f (x) The solutions of ( 12) and ( 13) are used to determine whether the eigenvalue is an isolated eigenvalue in order to generate δ, which is used in (11).
The implementation of the algorithm is shown in Algorithm 1 [44].
Step 1. Solve the semidefinite relaxation (9) by using Sedumi [49].If there is no solution, A has no Z-eigenvalue and stop; if not, compute a minimizer y * .
Step 6. Solve the relaxation (11) by using Sedumi.If there is no solution, λ i is the largest Z-eigenvalue and stop; if not, compute a minimizer y * for it.
A concrete example below is provided to illustrate Algorithm 1.Consider the symmetric tensor A ∈ R [4,3]  Using Algorithm 1, we can obtain all the Z-eigenvalues of A. The Z-eigenvalues are shown in Table 1, where λ (n) k means λ k has n repeated roots and n = rank M t (y * ).The computation task takes about 4.56 s.

Homotopy Continuation-Type Algorithm
Compared with other calculation algorithms, the homotopy method has its own advantages, like less restrictive conditions, fewer operation costs, more universal methods to compute all common types of tensor eigenvalue, and so on.In order to use the homotopy method, we define the equivalence class where T are unknown, while a 1 , . . . ,a n , b are random numbers.
The main idea of the homotopy method is to convert the general polynomial system T(x) = 0 into another polynomial system Q(x) = 0 that is easy to solve.In particular situations, the homotopy H(x, t) = 0 has smooth solution paths parameterized by t for t ∈ [0, 1), and all the isolated solutions of P(x) = 0 can be reached by tracing these solution paths.A useful homotopy is the linear homotopy [50]: where γ is a random nonzero complex number.Another common homotopy is the polyhedral homotopy [51,52] based on Bernstein's theorem: where h i (x, t) = (1 − t)γQ i (x) + tT i (x) and T i (x) is the i-th polynomial in a polynomial system.We mainly focus on the homotopy method to solve the tensor E-eigenvalue computing problem.The algorithm for computing the E-eigenvalue by the homotopy method is provided as follows.Firstly, with the polyhedral homotopy method, we can obtain the equivalence class T(λ, x).The solution of T(λ, x) is relevant to the tensor eigenpairs, and (λ, x) is called equivalence eigenpairs.The polynomial expression T(λ, x) for the E-eigenvalue is where A ∈ R [m,n] , B is the identity matrix I n ∈ R n×n .Then, we obtain all equivalent eigenvalues and eigenvectors from the equivalence class by using NAClab [53] based on the PSolve [54] method, which is widely applied in sparse matrix factorization such as [55,56].NAClab, a Matlab toolbox realizing the PSolve method, provides us a package in the numerical solution of polynomial systems by the homotopy continuation method.Finally, we find all E-eigenpairs by using the correspondence between equivalence eigenpairs and E-eigenpairs.The implementation of the algorithm is shown in Algorithm 2 [34].
Some concrete examples are given below.Consider the symmetric tensor A ∈ R [4,3] , which is same as the example A given in the semidefinite relaxation algorithm.Using Algorithm 2, we likewise obtain all the E-eigenvalues and E-eigenvectors of A. Then, we can filter the real part of E-eigenpairs to obtain Z-eigenpairs.The Z-eigenvalues are shown in Table 2.The computation task takes about 0.28 s, which is a significant reduction in calculation time compared with the previous algorithm.

Algorithm 2 Compute all E-eigenpairs of A
Step 1.Using modified PSolve to obtain all equivalent eigenvalues and eigenvectors (λ, x) of T(λ, x).
The computation task takes about 0.37 s.

Calculation Performance Analysis
In this section, we will evaluate the calculation performances of the above two algorithms from several different aspects, such as the ratio of convergence, accuracy, and CPU time.All the numerical experiments were conducted on a PC with an Intel(R) core (TM) CPU at 3.00 GHz, 8 GB of RAM, and Windows 10.The packages Tensor-Toolbox [9] and TenEig-2.0[34] were run using Matlab 2020a.
Firstly, the function "tenrand" in Tensor-Toolbox was used to generate the thirdorder two-dimension tensors to the fifth-order three-dimension tensors with 5000 sets across a total of 50,000 sets, in which the values of tensors follow the standard normal distribution.The generated dataset used Algorithm 2, "eeig" in TenEig-2.0,to calculate all the E-eigenpairs.As a matter of convenience, the non-convergence ratio, which is used to characterize the ratio of convergence, is defined as the number of error tensors, reporting the error during the computing over the total number of tensors.Among the 50,000 sets of data in this experiment, there were four tensors that did not converge.The non-convergence ratio is 8 × 10 −5 , and the specific occurrence of non-convergence is shown in the following Table 4.
The accuracy of the algorithm is divided into an estimation bias, which is known as residual, and upper bound bias.The residual is defined as the difference between the actual truth value and the estimated fit value, that is, the value of X x m i − λ i , where X ∈ R [m,n] is an arbitrary tensor and (λ i , x i ) is the eigenpair of X .For E-eigenpairs of X ∈ R [m,n] , the upper bound is ((m − 1) n − 1)/(m − 2) [57], which means that X has ((m − 1) n − 1)/(m − 2) equivalent eigenpairs if the number of E-eigenpairs are finite, and we use the notation E[m, n] to represent the upper bound.We use the PSolve method in Algorithm 2 to obtain equivalent eigenpairs directly, and in the third step of the algorithm, the equivalence eigenpairs are transformed into eigenpairs, so the actual upper bound in theory is 2E[m, n], which is called E ac [m, n].Based on the previous ratio of convergence test, the residual and upper limit data are obtained at the same time.
From Table 5, we can see that the residual is roughly 15 decimal places, which is negligible.The upper bound in theory is the same as the calculated upper bound.Finally, the computing costs represented by CPU time are shown in Table 6.
For small-scale tensors, where the order and dimension are less than four, the computing cost is linearly dependent on the upper bound with each eigenpair costing approximately 0.01 s.For large-scale tensors, affected by the "curse of dimension", the upper bound increases rapidly, and the cost of computation increases extremely fast.The upper bound of the fifth-order five-dimension tensor is 682, and the CPU time is 30.28 s, while the upper bound of the sixth-order six-dimension tensor is 7812 and the CPU time is increased to about 2878 s, where the unit cost increases to 0.37 s.

Spectrum Situation Awareness Based on Tensor Eigenvalues
Spectrum situation awareness can be divided into three stages: perception (sensing), understanding, and prediction.Spectrum perception is the primary task to know the spectrum usage situation, which can be fulfilled by various sensors.Due to the unideal environment with noise, it is necessary to use heterogeneous information to improve sensing performance.
For SAGINs with the tensor big data model, it is indispensable to solve two significant problems before the spectrum situation awareness; one is the storage overhead problem, and the other is the data missing problem.These two problems can be well solved by tensor decomposition and tensor completion, which were mentioned above.For a big data tensor A ∈ R n 1 ×n 2 ×...×n m , the required storage space is n m , and the overhead is unacceptable when m and n are large.Based on this, CP decomposition can greatly reduce the storage overhead.The specific algorithm [9] is as Algorithm 3.

Algorithm 3 CP decomposition algorithm for big data tensors
Input: the big data tensor A, the CP-rank R, the tolerance γ 0 , and the maximium number of iterations N 0 .
In the actual scenario, it is inevitable to avoid partial data loss because of sensor failure, transmission loss, and other reasons, which is particularly common for big data.For spectral big data tensor X , the missing value problem can be solved by tensor completion based on CP decomposition, Tucker decomposition, and the minimum trace norm.Representation by the optimization problem can be written min where Θ is the set of nonzero-valued indexes in X , and || • || * is the tensor trace norm, which is defined as ||X || * = ∑ m i=1 α i ||X (i) || * with ∑ m i=1 α i = 1, which can be regarded as the weighted sum of the n-mode unfolding matrix traces.All of these problems can be solved using the block coordinate descent algorithm, using the solution of (19) as an example to illustrate the tensor completion algorithm, as shown below Algorithm 4.

Algorithm 4 CP-based completion algorithm for big data tensors
Input: the big data tensor Y ∈ R n 1 ×•••×n m , the tolerance γ 0 , and the maximium number of iterations N 0 .
Step 0. Initialize A 1 • • • A m by using random numbers, Θ and its complementary set Θ, let X Θ = Y Θ , N = 0 and k = 1, and compute Output: the completed tensor X .After solving the above problems, we try to use tensor eigenvalues to construct a spectrum situation awareness scheme.To the best of our knowledge, it is the first to use tensor eigenvalues to evaluate spectrum situation awareness.Similar to classical signal detection methods, the proposed situation awareness scheme is based on a binary hypothesis test where x(t) denotes the received signal, s(t) is the target signal, and n(t) is the noise.Hypothesis test results are H 1 and H 0 .
Based on the spectrum situation awareness model, the signal tensor is generated with the target signal and noise.The eigenvalues of the signal tensor are used to construct the detectors, and the sensing results are generated by comparing the detector with the given threshold.For a given specific false alarm, the thresholds can be determined by the signal tensor with only noise H 0 .In polynomial time, the eigenvalue of the signal tensor is calculated with the homotopy method.Then, the detection performances can be evaluated by comparing such an eigenvalue with the given threshold.If the detector is greater than the threshold, the state is H 1 , indicating that there is a target signal.Otherwise, the state is H 0 , indicating no signal.
In order to demonstrate the performance of the algorithm, we simulated 52 sets of spectrum tensors with different tensor sizes and SNRs by Matlab; each set consisted of 10,000 received signal tensors and 1000 noise tensors.The detection performances of Algorithm 5 are shown in Figure 3 with varying tensor sizes, SNRs, and P α .In Figure 3a, the probability of detection P d over 25 samples is plotted against the SNR under the various tolerated probabilities of false-alarm P α .It is found that with the increase of SNR, P d gradually increases and finally reaches 100%.For the same threshold, as the SNR increases, the maximum eigenvalue increases, and the detection probability also increases.This phenomenon shows that this algorithm is effective for spectrum detection.When the SNR is sufficient, the detection success probability is close to 100%; that is, the maximum eigenvalue can be used to represent the existence of average energy in the tensor.Moreover, we noticed at the same time that P α mainly affects the speed of P d to reach 100%.The larger P α is, the smaller SNR is when P d reaches 100%.When P α = 0.01, SNR = 5 dB, that P d reaches 100% for the first time, but when P α = 0.10, SNR = 2 dB.The reason for this phenomenon is that P α directly affects the value of the threshold.When P α increases, the threshold decreases, resulting in the tensor eigenvalue of a lower SNR being higher than the threshold.

Algorithm 5 Algorithm of eigenvalue-based spectrum situation awareness
Input: the noise tensors N , the signal tensors X , and the tolerated probability of false alarm P α .
Step 0. Compute the E-eigenvalues of N by Algorithm 2, and obtain the noise eigenvalue distribution.
Step 1. Compute the threshold T for a given P α based on noise eigenvalue distribution.
Step 2. Compute the E-eigenvalues of X and obtain the maximum eigenvalue (module value) S.
Step 3. Compare T and S. If T ≥ S, consider X to be the noise tensor.Otherwise, X is the signal tensor.Output: the result of spectrum situation awareness.Based on the above findings, the detection probability will increase if the tolerated probability of false alarm increases and vice versa.In order to illustrate this effect, the algorithm is compared in terms of the probability of detection P d as a function of the tolerated probability of false alarm P α in Figure 3b.When P α < 0.15, P d increases rapidly, while when P α > 0.15, P d increases gradually, and the curve of the larger SNR is above that of the smaller.The former is because, in hypothesis testing, the threshold distribution obeys the Gaussian distribution.In the first part, the threshold decreases rapidly with the increase of P α , and the detection success probability increases rapidly with the same SNR.In the second part, due to the concentrated distribution of the threshold, the threshold changes insignificantly with the increase of P α , resulting in a slow change of detection probability.The latter is because the larger SNR makes the desired threshold for P d = 100% higher, which is self-consistent with the conclusion drawn by Figure 3a.
Figure 3c,d are mainly to illustrate the effect of tensor size on the algorithm.Instead of using tensors X ∈ R [3,3] in Figure 3a,b, we use tensors Y ∈ R [3,5] in Figure 3c,d, where each tensor grows from 27 elements to 243 elements.From the comparison, it is easy to find that the effect of the algorithm for Y is significantly better than that for X .When P α = 0.01, the SNR just has to be equal to 2 dB in order for P d to reach 100%.The slope of the front part of the curve is significantly higher than that of Figure 3a.
Through subsequent simulation experiments on tensors of different order and different dimensions, it seems to be concluded that for tensors of the same order, the larger the dimension is, the better the algorithm effect will be.For the same dimensional tensor, the larger the order is, the worse the algorithm is.Explaining the cause of this phenomenon is the key problem to be solved in the following work.However, in a word, the eigenvalue detection method for spectrum sensing can effectively solve the problem of spectrum signals.

Conclusions
A novel tensor-based spectrum situational awareness scheme has been proposed, and the tensor was used to model and process multidimensional, heterogeneous, and massive spectral big data from space-air-ground integrated networks.In particular, the tensor eigenvalues have been utilized to construct the spectrum situation awareness, in which the distributions of E-eigenvalues were included to formulate spectrum sensing algorithms.Two tensor eigenvalue calculation schemes have also been provided to generate tensor eigenvalues.Simulation results have evaluated the correctness of the proposed situational awareness scheme.The situation awareness scheme can detect the spectrum state quickly, accurately and coarse grained, providing valuable support for subsequent understanding and prediction.However, as the complexity of tensor eigenvalue computation increases catastrophically with the increase of tensor size, this scheme can only deal with local small-scale data of spectral big data.In the future, the research will focus on exploring the correlation between tensor decomposition and eigenvalue, aiming to enable large-scale tensor eigenvalue computation.Additionally, further investigation in tensor computing will be conducted to enhance the subsequent operations of situational awareness.Overall, the novel tensor-based spectrum situational awareness scheme has provided a new application paradigm for tensor theory.

Figure 2 .
Figure 2. Tensor representation model and tensor-based spectrum situation awareness model.

Figure 3 .
Figure 3. Probability of detection of spectrum sensing algorithms based on tensor eigenvalue, as a function of the SNR and the tolerated probability of false-alarm.(a) Performance against SNR; (b) performance against P α ; (c) performance against SNR; (d) performance against P α .

Table 4 .
The non-convergence number in 5000 sets of tensors by homotopy continuation algorithm.

Table 5 .
The accuracy in 5000 sets of tensors by the homotopy continuation algorithm.

Table 6 .
The time of computing in 5000 sets of tensors by homotopy continuation algorithm.