Abstract

In this paper, the optimal beamforming problem of multi-input single-output (MISO) cognitive radio (CR) downlink networks with simultaneous wireless information and power transfer is studied. Due to the nonconvexity of the objective function, the considered nonconvex optimization problem is firstly transformed to an equivalent subtraction problem and then an approximated convex optimization problem is obtained by using the successive convex approximation (SCA). When the instantaneous channel state information (CSI) of the eavesdropping link is unknown to the legitimate transmitter, another interruption-constrained energy efficiency optimization problem is proposed and the Bernstein-type inequality (BTI) is used to conservatively approximate the probability constraint. The paper proposes a two-level iterative algorithm based on Dinkelbach to find the optimal solution of the EE maximization problem. Numerical results validate the effectiveness and convergence of the proposed algorithm.

1. Introduction

In order to alleviate the spectrum scarcity, CR has been proposed as one of the most promising solutions. As long as the interference caused by the secondary user (SU) is tolerable to the primary user (PU), the secondary user can make use of the licensed spectrum to guarantee the quality of service [1]. As cognitive radio can improve spectrum efficiency, it has been widely used in relay networks and cellular networks [2].

Different from the traditional energy resources, in the energy-constrained wireless communication system, the radio frequency (RF)-based simultaneous wireless information and power transfer (SWIPT) [3] technology is more suitable for scavenging energy from the pervasive radio frequency signals. Radio frequency signals can simultaneously transmit energy and information by the electromagnetic waves. Energy efficiency (EE) is considered to be an effective metric for balancing the power use and spectrum efficiency [4], so a lot of work has contributed to the maximization of EE in CRNs. In [5], the actual nonlinear energy collection model is considered in a nonorthogonal multiple-access cognitive radio networks. Y. Wang et al. developed a multiobjective resource optimization problem to maximize the sum of the collected power from all the receivers.

However, due to the open broadcasting feature of SWIPT cognitive radio, the confidential information transmitted is easy to be eavesdropped by energy-harvesting users, so that the network security is threatened by these potential eavesdroppers. Therefore, physical layer security (PLS) based on information theory is hence proposed to guarantee the security of wireless communication systems [69]. In order to further protect the legitimate communication using PLS, artificial noises and beamforming based on zero forcing are the popular approaches.

Under the constraints of information leakage and energy efficiency, the high security of full-duplex SWIPT system is ensured by formulating and maximizing information transmission rate in [10]. Ouyang et al. in [11] investigated a tradeoff between the secrecy rate and energy efficiency by designing zero-forcing beamforming. In [12, 13], Zhang et al. designed a beamforming technology in MISO cognitive radio networks and studied different secure energy efficiency resource optimization problems under the constraints of interference power leakage and the constraints of energy collection in the secondary network.

The security resource allocation problem is designed in the case of ideal CSI in [1013]. However, in practice, due to the channel estimation errors, it is difficult to achieve the confidentiality requirements without the security scheme design in the physical layer. In particular, in secure cognitive radio networks with SWIPT, imperfect CSI will significantly degrade the performance of beamforming. Therefore, the academic contributions based on the CSI error model have received a lot of attention [14]. At present, there are totally two types of beamforming designs to achieve robust security. The first type bounds the quantization error of CSI into a confined range. Although it is of low complexity, the actual performance may be underestimated and is not suitable for delay-sensitive cognitive radio networks.

The other type is the statistical CSI error model, which can well approximate to the constraint based on the outage probability and is also an effective measure in realizing the robust security beamforming. In [15], the problem of secret energy efficiency maximization under perfect CSI and statistical CSI is studied, respectively, and a robust beamforming design based on the ellipsoidal channel uncertainty is considered, but the computational complexity is too high to be realized in practice.

This paper focuses on the radio frequency-based wireless power transmission and wireless energy harvesting in MISO CR networks. The communication spectrums are shared between the primary users and secondary users. The multiantenna secondary user simultaneously transmits confidential information and energy to the secondary receiver and several energy receivers. Based on the omnidirectional nature of wireless transmission, the energy receiver (ER) can be regarded as a potential eavesdropper. The paper proposes an energy efficiency maximization problem through optimizing artificial noise-assisted beamforming, i.e., to maximize energy efficiency while satisfying the requirements of quality of service. The main contributions of this paper are summarized as follows:(1)Given the assumption that CSI is perfect, the beamforming design problem can be incorporated into the formulated EE maximization problem, which satisfies the constraints, such as confidentiality, interference leakage, and transmission power. Due to the objective function being a nonconvex fraction, the global optimal solution cannot be explored with ease directly. In order to solve the nonconvex optimization problem, it needs to be firstly transformed to an equivalent subtraction form and then to be applied with the successive convex approximation. A two-layer iterative algorithm is developed for solving the transformed optimization problem.(2)When the perfect CSI is unavailable, the robust design based on the statistical CSI error model is considered and the EE maximization problem is constructed under the constraints of required secrecy rate and the confined interference probability. Because the probability constraint has no closed expression, the formulated optimization problem is still nonconvex. Therefore, the original problem is transformed to an equivalent problem with closed expression through the Bernstein inequality. By introducing new instrumental variables, a two-stage optimization algorithm is designed to calculate the corresponding optimal beamforming vector.

The rest of the paper is organized as follows: Section 2 describes the network model. Section 3 formulates the EE maximization problem considering the perfect CSI model. In Section 4, the problem of robust EE maximization adopting the statistical CSI error model is studied. In Section 5, the simulation results are given to evaluate the effectiveness of the proposed scheme. The conclusion is summarized in Section 6.

2. Network Model and Problem Formulation

2.1. Downlink Channel Model

As shown in Figure 1, the paper considers a downlink underlying CR network with multiple primary users and eavesdroppers. The primary network composes of a primary transmitter (PT) and several primary receivers (PRs). When using the spectrum sharing paradigm, the primary network is coexisted with the secondary networks [16], which consists of a secondary transmitter (ST), a secondary receiver (SR), and secondary energy receivers (SERs). The PT is equipped with antennas, while the ST is equipped with antennas. PRs, SR, and SERs are all equipped with single antenna. This article assumes that ST transmits confidential information to the SR and SERs harvest RF energy from ST’s transmission using wireless power scavenging technology. However, due to the omnidirectional propagation of electromagnetic wave, SERs are also the potential eavesdroppers to intercept and monitor the transmitted information from ST to SR. To improve the security of the secondary network, ST attaches artificial noises to the transmitted data to degrade the SERs’ channel quality.

The secondary transmitter will simultaneously transmit artificial noises and information signals in order to ensure communication security and to promote effective power transmission. The transmitted signals from PT and ST can be separately modeled as follows:where and are the data signal of PT and the corresponding beamforming vector, respectively. and are the confidential information of ST and the corresponding beamforming vector, respectively. Without loss of generality, we can obtain that and . is the noise vector generated by ST for improving the secrecy rate of SR. The AN vector is usually modeled as a circularly symmetric complex Gaussian vector with mean 0 and covariance matrix .

Thus, the received signals at PR, SR, and SERs can be given as follows, respectively:where and . represents the channel-fading coefficient vector for the links from the ST to the PRs, SR, and SERs. In addition, represents the channel-fading coefficient vectors for the links from PT to PRs, SR, and SERs. Moreover, denotes the additive white Gaussian noise (AWGN) with zero mean and variance at the PRs, SR and SERs, respectively, which include the thermal noise and signal processing noise at all types of receivers i.e., PRs, SR, and SERs.

2.2. Performance Metric

The instantaneous rate of the link between the ST and SR can be expressed as follows:

According to [17], the channel data rates of the link between the STs and SERs are given as follows:

The secrecy rate is usually introduced according to the service quality requirement of system design and is used to ensure secure communication [18]. In particular, the secrecy rate is defined as the maximum of authorized data rate that the eavesdroppers cannot decode even though its power supply is enough. Thus, the secrecy rates transmitted from ST to SR for ensuring the confidential information [19] can be expressed as the following equation:

In addition, the total transmission power at the ST for downlink transmission can be expressed as follows:where represents the circuit power consumed of each transmit antenna at the ST and represents ST’s fixed power consumed.

Finally, to reach a good tradeoff between the sum achievable rate and the total power consumption of the secondary CR networks, EE is used to evaluate the degree of such balance. The energy efficiency in the paper is defined as the ratio of the transmission rate of the system to the total power consumption [20] and is given by the following equation:

3. EE Maximization with Perfect CSI

In this section, the perfect channel state information (CSI) of all terminals is assumed to be available at SR [21, 22]. According to the secure transmission requirements of SR, the interference power control of PRs, and the power restraint of ST, the optimal transmission covariance matrix can be obtained by solving the EE maximization problem marked as P1.where is the secrecy rate threshold required by the ST. The interference power at PR should be lower than to achieve the demanded quality of service for PRs, and represents the maximal transmission power of the ST. Equation (11b) guarantees that the secrecy rate of SR should be no less than . Equation (11c) implies that the interference power control demand for protecting the quality of service of PRs. Equation (11d) indicates the maximal allowed transmission power of the ST.

This problem is nonconvex due to the fractional form of the objective function shown as equation (11a). In order to solve the problem P1, it is firstly converted to parametric programming problem with respect to [23], shown as P2.

Theorem 1. Let be the unique zero solution to , where and denote the optimal secrecy rate and the sum power consumption, respectively. has the same optimal solution to the original optimization problem P1. The maximum of energy efficiency can be given in equation (13), if and only if is satisfied.

Proof 1. Please refer to Appendix A for the proof of Theorem 1.

Theorem 1 shows that if satisfies , then the fractional optimization problem P1 can be solved using the optimal solution of problem P2 shown as equations (12a)-(12b). To obtain the optimal solution to problem P2, inequality constraint (11b) can be further separated to two inequalities, i.e., (14b) and (14c), so the optimization problem P2 can be equivalently transformed to the following problem P3:where is the maximum required signal-to-interference-plus-noise ratio at the k-th SER in the secondary network and and .

Through the variable substitution by defining , , and , problem P3 can be further reorganized as the following problem P4:

Due to the logarithmic functions being concave, the objective function remarked as equation (15a) is now converted to the difference of two concave functions (DC) and is still a nonconvex function.

Let

Then, in order to solve the nonconvex objective function, the DC approximation method [24] can be used to approximate at the feasible solution. By applying the first-order Taylor series expansion, as the gradient of at can be given by the following equation:

Then, according to the Taylor series expansion of , equation (18) can be obtained:

Finally, by substituting into equation (15a), the problem P4 can be transformed to an approximated optimization problem P5 shown as follows:

However, due to the existence of rank-one constraint marked as equation (15g), problem P5 is still a nonconvex optimization problem. Theorem 2 is proposed herein to guarantee the existence of the optimal solution to problem P5.

Theorem 2. If problem P5 is feasible, the rank-one optimal solution always exists.

Proof 2. Please refer to Appendix B for the proof of Theorem 2.

Based on the conclusion disclosed by Theorem 2, the rank-one constraint numbered as equation (15g) can be reasonably ignored, so the approximated problem P5 is a convex optimization problem. Even the rank-one constraint equation (15g) is removed from P5, and it has no effects on the optimal solution. Hence, the optimal solution can be obtained by solving problem P6 shown as follows:

Therefore, it is straightforward to obtain the optimal solution of P6 through the existing CVX software tools [25]. Algorithm 1 presented in the pseudo code diagram summarizes the algorithm procedures based on the DC programming approach.

(1)Initialize and choose an initial value ;
(2)Repeat—Outer loop (Dinkelbach’s algorithm)
(3)  Initial , and , ;
(4)  Repeat—Inner loop (DC programming)
(5)    Find the optimal solution of problem P6 for obtained
(6)   , and by using CVX;
(7)    Intermediate variable update ;
(8)    Set ;
(9)  Until;
(10)  Update ;
(11)  Set ;
(12)Until, and is the convergence threshold;
(13)Return ;
(14)Output .

4. EE Maximization with Statistical CSI

For high-rate applications, data services such as high-fidelity video and real-time voice data transmission are both sensitive to time delay, so the secure beamforming design that relies on the perfect CSI model is not suitable in such scenarios anymore. In this section, ST has difficulty to obtain perfect CSI for PRs, SRs, and ERs, of which they also do not acquire the channel state information before their deployments.

As to the statistical CSI model, in contrast with perfect CSI, the channel state error is random and follows a certain distribution. In the paper, the complex circular Gaussian CSI error model is applied [26] and then the channel vectors can be given as follows:where denotes the estimation of actual CSI, which is known at the ST, and is the corresponding statistical errors. and are independent to each other for any . and are also independent to each other for any . Moreover, denotes the covariance matrices of the corresponding channel error vectors, where [27]. is the variance of the corresponding channel uncertainty. Such assumption for channel uncertainty has been widely adopted in literature studies about the design in physical layer security.

4.1. Outage-Constrained Robust Problem Design

This section aims at maximizing the EE of CRN while satisfying the given outage constraints. The optimization problem can be formulated as the following problem P7:where indicates the maximum of outage probabilities associated with the secrecy rate at SR. The outage secrecy rate constraint (24b) guarantees that the probability of the secrecy rate above should be no less than . And, is the maximum tolerable probability that the interference power caused to the -th PR. The outage interference power constraint presented by inequality (24c) means that the probability of the interference power of the PRs lower than should be no less than .

In the secondary network, the channel quality is guaranteed by using the orthogonal frequency division multiplexing. is designed to be orthogonal to and , respectively, and hence, and can be obtained. By introducing an intermediate variable , problem P7 can be reshaped as the following problem P8:

Since the objective function is a fraction, of which the numerator is in the form of a logarithmic function, problem P8 is also a nonconvex optimization problem. However, P8 is more challenging to be handled with, not only since equations (25b)–(25d) have no closed forms but also a rank-one constraint shown as equation (25g) also exists.

4.2. Beamforming Design with the BTI-Based Approach

The outage probability indicated in equation (25b) can be seen as the coupling of links between SERs. All the SERs’ channels are assumed to be independent, the secrecy rate in the constraint of outage probability shown as equation (25b) can be decoupled as inequality inwhere . The terminal inequality in equation (26) can be regarded as the constraint of secrecy outage probability for each SER [28]. Equation (26) can be further transformed to the following equation:where and .

By equivalent transformations, equation (28) can be obtained and shown as follows:where , .

Since the occurrence of link outage is probabilistically constrained, it is hard to derive a closed-form expression for equations (25a)–(25g). In order to transform the probability constraint into a closed form, Bernstein’s inequality is used to approximate the uncertain constraint. It is worthy to note that the optimal solution obtained by such approximation always satisfies the interruption constraint [28].

Lemma 1. (the Bernstein-type inequality). Suppose the following constraint:where is a set of optimization variables, , and is a fixed constant. By introducing two slack variables , inequality (29) can be equivalently represented by a group of inequalities [29]:

Based on the result of Lemma 1, the approximated and convex form for inequality constraint (28) is now transformed to a set of inequalities, shown as the set of inequalities in where both and are the slack variables.

Let , in which . Applying Lemma 1 to constraint equation (25c), the outage interference power constraint can be approximated to a set of inequalities:where and are the slack variables, , and .

Similarly, let , in which . According to the conclusion of Lemma 1, the outage interference power constraint, i.e., inequality constraint (25c), can be approximated as the following set of inequalities:where and are the slack variables.

Although the interruption constraint of the original problem has been transformed through the above approximation of the probability constraints, the problem is still nonconvex due to the existence of the fractional objective function and its nonconvexity caused by the coupling relations of optimization variables. Now, the original problem P8 can be rewritten as problem P9:

Because the objective function in P9 is a fraction with nonconvex terms, the problem P9 is still a nonconvex optimization problem. Following a similar approach revealed by Theorem 1, which has been used to transform equation (11a), the problem P9 can be transformed to a parametric problem with respect to , problem P10, as follows:

By introducing a new relaxation variable , the optimization problem P10 can be equivalently transformed to problem P11.

The problem P11 can be further reorganized as a two-stage optimization problem by introducing a new relaxation variable , which makes P11 decomposable. The outer problem is redefined as problem P12.where is the minimal transmission power that makes problem P12 feasible. In addition, the inner problem determines for a fixed , which can be marked as problem P13.

As can be seen from equations (37a)-(37b), the outer problem is a single-variable optimization problem, of which the optimal objective function value can be obtained by performing a one-dimensional line search, such as the golden section search method.

Then, by applying semidefinite relaxation (SDR) approach to equation (34d), the nonconvex rank-one constraint can be removed [30] and then the relaxed inner problem P13 can be rewritten as an optimization problem P14.

Remark 1. Provided that problem P14 is feasible, the SDR approach is tight and the global optimal solution satisfies , which can be proved to be rank one in a way similar to [30].

In the following section, the inner problem after relaxation technique can be solved by using Charnes–Cooper transformation [31] and S-procedure [32].

4.3. Reformulation of the Inner Problem

By introducing the Charnes–Cooper transformation shown asand a slack variable , the inner problem P14 can be transformed to problem P15 shown as follows:

Lemma 2 is necessary to transform P14 to a linear matrix inequality (LMI) constraint.

Lemma 2. (S-procedure). Define the functionwhere , and . Provided that there exists a vector such that , the implication holds if and only if there exists a satisfying the following equation:

By applying Lemma 2, inequality (41b) can be written as follows:where and denotes the channel uncertainty region, which can be modeled by the inverse cumulative distribution function (CDF) of the Chi-square distribution with degrees of freedom, e.g., . Therefore, by applying the Cauchy–Schwarz inequality to the channel uncertainty region, inequality can be obtained.

According to the S-procedure given by Lemma 2, inequality constraint (41c) can also be converted to the following inequality:where .

Synthesizing the above analysis, the inner problem P15 can be reshaped as problem P16.

At present, the specific numerical simulation toolbox (such as CVX) can be used to effectively obtain the optimal solution of the problem P16. The detailed procedures for solving the optimization problem P7 are summarized as the pseudo code of Algorithm 2.

(1)Initialize and choose an initial value ;
(2)Repeat
(3)  Perform the one-dimensional line search over to obtain the optimal variable values of the outer problem P12, wherein each is obtained by solving the inner problem P16 via the interior-point method;
(4)  Update the corresponding through equation (40);
(5)  Compute ;
(6)  Update ;
(7)  Set ;
(8)Until satisfied;
(9)Return the optimal value of the inner problem in P16 to the outer problem in P12;
(10)Return as the optimal solution.

5. Simulation Results and Performance Evaluations

In this section, the simulation results are given here to evaluate the performances of the two proposed schemes. The secure beamforming scheme with perfect CSI and the robust secure beamforming scheme with the statistical CSI error model are thoroughly compared in different network scenarios.

In the considered underlying MISO CR networks, the primary transmitter and the secondary transmitter are equipped with antennas, respectively, while SR, PR , and ER are uniformly equipped with a single antenna. In the secure beamforming scheme with perfect CSI, the channel coefficient on each wireless link is generated by CSCG and its mean and unit variance are zero. Other simulation parameters are configured as follows: , , and . The convergence tolerances and are both set to be .

In the robust secure beamforming scheme with the statistical CSI error model, all channel coefficients are modeled as being obeyed with Rayleigh fading, which can be modeled as follows:where denotes the line-of-sight (LOS) deterministic component with and indicates the Rayleigh fading component. stands for the simplified large-scale fading mode, where is the reference distance, represents the actual distance between the receiver and the transmitter in the secondary network, and the path loss exponent . is the Rician factor [33]. In addition, the secrecy rate outage probability and interference power outage probability remain constant by setting .

Figure 2 shows the convergence tendency of the two proposed schemes, i.e., optimal secure beamforming with the perfect CSI model and with the statistical CSI model. Experimental results show that the proposed two schemes can converge to their own optimal EE after a limited number of iterations. Due to the complexities of the two schemes being distinct, the perfect CSI model scheme is first to converge to its optimal solution. Although the beamforming scheme with the statistical CSI model needs more iterations, the energy efficiency can also converge to its optimum. In addition, as can be seen from Figure 2, higher energy efficiency can be achieved when the basic power consumption of the two algorithms is higher. This means that the increase in basic power consumption of ST can improve the metric EE in these two schemes.

Figure 3 illustrates the optimal EE versus the number of ST’s transmitting antennas by varying the secrecy rate. As shown in Figure 3, the optimal EE in the two proposed schemes, i.e., perfect CSI and statistical CSI, decreases when the secrecy rate rises. The higher security rate requires more transmission power to improve the physical layer security, resulting in a lower optimal EE achieved in the CRNs. In particular, when the secrecy rate is set as 3.3–3.5 bps/Hz, the optimal EE in the statistical CSI schemes degrades to zero, which means that the optimization problem is not feasible anymore as the required secrecy rate is too high to be achieved. On the other hand, it can also be seen from Figure 3 that the value of optimal EE in the perfect CSI improves when increasing the number of transmitting antennas at ST. This is because the diversity gain is directly proportional to the number of transmitting antennas.

Figure 4 shows that the achieved optimal EE of the considered two different CSI models can converge by changing the total transmission power. It can be seen in Figure 4 that when the total transmission power ranges from 20 dBm to 30 dBm, the optimal EE in two schemes increases when the total transmission power rises and remains unchanged after the maximum of EE is reached. This is because in order to maintain the maximally achievable EE, the proposed two schemes reduce the energy consumption by canceling the continued allocation of transmission power. In addition, the average EE of the proposed scheme increases when the number of ST’s transmitting antennas grows. This can be explained by the fact that by increasing the diversity gain of the channel through the increase of , more transmission power is used to meet the secrecy rate requirement and the achievable EE trivially increases finally.

Figure 5 shows the relationship between the available maximal EE and the number of SERs in the two CSI model schemes. It can be seen from Figure 5 that the optimal energy efficiency of CRNs decreases as the number of SERs grows. This is because the proposed scheme needs to be designed to meet the secrecy rate constraint and the secrecy rate is affected by the number of SERs. The results show that when the amount of SER is no less than 6 and , the problem of maximizing EE with the perfect CSI model is infeasible anymore, because the network cannot meet the security rate constraint of the secondary receiver and the transmission power already exhausts.

6. Conclusion

In this paper, the artificial noise technology is applied to improve the proposed physical layer security for CRNs. The EE maximization problem in the underlying MISO CR networks with the models of perfect CSI and statistical CSI is formulated in a network scenario, where multiple PRs and SERs are involved. Due to the existence of the fractional objective function, the original optimization problem is nonconvex. By using the nonlinear fractional programming and the DC programming, the considered nonconvex problem is transformed to an approximated convex optimization problem. When the statistical CSI error model is applied, a robust secure beamforming strategy is designed. The outage probability constraint is converted to a closed expression, and a two-stage iterative algorithm is proposed by using SDR and BTI to explore the approximated solution. Simulation results verify the convergence of the proposed algorithms; hence, the optimal security beamforming vector and the optimal AN covariance matrix can always be obtained.

Appendix

Proof of Theorem 1

It is obvious that the problems P1 and P2 have the same feasible region when given the constraints defined in equation (12b). Firstly, define , , and as the optimal solution to the problem P1 and as the set of feasible solutions. Then, (A.1) can be derived as

Based on the fact that , there exist the following conditions:

From the inequality in (A.2), equation (A.3) is obtainedand from (A.2), it is clear that the maximum value of the objective function of the problem P2 can be achieved at the optimal solutions , , and . The first part of the proof has been completed, which satisfies the sufficient condition.

As the proof by contradiction, the correctness of necessary condition can be achieved by the following procedures. Let be the optimal solution of (A.3) and , (A.4) can be concluded as follows:

Through equivalent operations, the following formula can be obtained:

From (A.5), it is obvious to find that is also the optimal solution and is the maximum objective value of the problem P1. The proof of Theorem 1 is validated.

Proof of Theorem 2

It can be proved that the problem P5 is a convex optimization problem with respect to the optimization variable and satisfies Slater’s condition, as well as its dual gap is zero [34]. The relaxed Lagrangian function of the problem P5 is given aswhere are the Lagrangian multipliers associated with constraints (15b)–(15g), respectively. According to the Karush–Kuhn–Tucker (KKT) conditions [35], equation (B.2) is hence obtained.

Based on these KKT conditions, the following equality holds:

Then, substituting in of equation (B.2), (B.4) can be obtained after equivalent transformations.

Hence, the following rank relation holds:where the first inequality is due to the basic rank inequality . By eliminating the trivial solution , holds. Similarly, can also be derived.

Data Availability

The original experimental data used to support all the charts of this study are contained in Section 5 of the article. The algorithm data used to support the findings of this study are contained in Sections 3 and 4 of the article. For details, see Algorithms 1 and 2. Readers can implement verification according to the pseudo code provided. The source program data used to support the findings of this study cannot be made freely available due to privacy or ethical restrictions. Reasonable request for access to these data should be made to the corresponding author.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was partially supported by the Natural Science Foundation of Hebei Province under Grant No. F2021402009 and the Scientific Research Plan Projects of Hebei Education Department under Grant No. BJ2014019. The work of M. Zhang was supported in part by the National Nature Science Foundation of China under Grant 62101080, the Science and Technology Research Program of Chongqing Municipal Education Commission of China under Grant KJQN202100738, and the Research Start Up Funding of Chongqing Jiaotong University under Grants 2020020070 and XJ2021000501.