随着现代航空航天和通信技术的不断发展，卫星通信因其具有覆盖范围广、通信容量大和不受地理限制等优点，已广泛应用于广播、救援、军事等多种场景^[1-4]。然而，由于卫星通信的广域覆盖特性和广播特性使得其传输的信号极易被非法用户窃听，因此如何保证合法用户的通信安全成为了卫星通信中的关键问题。不同于传统的基于计算复杂度的上层加密算法^[5-6]，物理层安全技术利用无线信道的随机性、互易性、信道间的差异性等特点从物理层角度增强信息传输的安全性^[7-9]。Wyner^[10]和Csiszar^[11]等的研究表明了当窃听信道的质量低于合法用户信道的质量时可以实现安全传输，该成果为物理层安全技术的研究奠定了基础。

基于多天线的波束成形（beamforming，BF）技术作为实现无线通信物理层安全的一种有效手段，既可以增强合法用户对信号的接收，同时也能有效地抑制窃听过程，受到了学术界的广泛关注^[12-15]。张杰等^[16]研究了在已知理想的窃听信道状态信息的情况下，以安全速率最大化为优化目标，提出了迫零和增强信漏噪比安全BF算法，并通过仿真得出了增强信漏噪比BF算法能获得更好的安全速率。由于在实际的通信场景中，发送端获取的信道状态信息存在估计误差和反馈时延等情况，理想的信道状态信息几乎无法获取。假设已知非理想的窃听信道状态信息，Lin等^[17]研究了单窃听场景下多波束卫星通信系统的安全传输问题，提出了卫星总发射功率约束下基于安全速率最大化准则的BF算法；王舒等^[18]针对卫星总发射功率受限情况下，单窃听场景中的鲁棒性安全速率最大化问题，利用发射波束成形向量与噪声协方差的联合优化进行安全BF设计。进一步地，考虑到存在多个窃听者的情况，Ma等^[19]研究了卫星发射功率受限和安全中断概率约束下的安全速率最大化问题，在已知统计信道状态信息的情况下，使用半定松弛法求得了最优BF权矢量；Lin等^[20]研究了针对星地融合系统中存在多个未知窃听者的情况，引入人工噪声和协作干扰信号对窃听者进行干扰，提出了一种联合BF算法实现安全传输；Lin等^[21]研究了考虑窃听者间合作与非合作的情况，在基站总功率受限和主用户干扰功率约束下，以安全速率最大化为优化目标，提出了相应的BF算法；顾晨伟等^[22]研究了卫星总发射功率约束下的系统安全速率最大化问题，在仅已知窃听者大概位置的条件下设计出了相应的安全BF算法。

上述研究大都考虑在总发射功率受限情况下，如何有效提升通信系统的安全性。考虑在多波束卫星中，每个馈源配置专用的射频链路和功率放大器，因此研究单波束功率约束下的安全问题更为实际。本文针对多波束卫星通信系统，在已知统计信道状态信息条件下，考虑窃听者间合作的情况，分别以卫星发射功率最小化和安全速率最大化为优化目标，提出了相应的BF方案。第1类方案，以安全速率为约束条件，卫星总发射功率最小化为目标函数建立优化问题，运用拉格朗日乘数法推导出最优BF权矢量的表达式。第2类方案，以安全速率最大化为优化目标，分别以卫星总发射功率受限和单波束发射功率受限为约束条件建立优化问题。针对以总发射功率受限为约束条件的优化问题，通过将问题转换成广义瑞利熵的形式，推导出最优BF权矢量的表达式；针对以单波束发射功率受限为约束条件的优化问题，先引用辅助变量将非凸优化问题转换为凸优化问题，再提出了一种基于黄金分割法的BF算法用于求解最优BF权矢量。最后，通过计算机仿真验证了所提方案的可行性与有效性。

1 系统模型

图1为工作在Ka频段的多波束卫星通信系统，由一个地球同步轨道（GEO）卫星、一个信关站和多个用户组成，其中，1个合法用户，另外K个卫星覆盖区域内的网内用户有可能窃听卫星发送给合法用户的私密信号而成为窃听者。不失一般性，卫星采用多馈源单反射面形式的天线，配置有L个馈源产生N个波束（ $L \ge N$ ），合法用户和窃听者均使用高增益的抛物面天线。

假设合法用户的数据信号为 $s\left( t \right)$ ，满足 ${\rm{E}}\left[ {|s\left( t \right){|^2}} \right] = 1$ ，信关站获取合法用户的数据信号 $s\left( t \right)$ 后，先采用BF技术进行处理，再通过馈电链路将处理后的信号 ${{w}}s\left( t \right)$ 发送给卫星，假设在馈电链路中进行理想传输，无噪声^[2,14]；该信号经过卫星信道后，合法用户接收的信号和第i个窃听者接收的信号分别表示为：

${y_{\rm{s}}}\left( t \right) = {{h}}_{\rm{s}}^{\rm{H}}{{w}}s\left( t \right) + {n_{\rm{s}}}\left( t \right)$

(1)

${\;\;\;\;\;\;\;y_i}\left( t \right) = {{h}}_i^{\rm{H}}{{w}}s\left( t \right) + {n_i}\left( t \right),\;\;i \in \{ 1,{\rm{2}}, \cdots ,K\} $

(2)

式中： ${{w}} \in {\mathbb{C}^{N \times 1}}$ 为BF权矢量； ${n_m} $ （ $ m \in \{ {\rm{s}},i\} $ ）表示均值为0、方差为 $\sigma _m^2$ 的加性高斯白噪声，噪声功率为 $\sigma _m^2 = $ $ \kappa {T_m}{B_m}$ ，其中， $\kappa \approx 1.38 \times {10^{ - 23}}\;{\rm J/K}$ 为玻尔兹曼常数， ${T_m}$ 为噪声温度， ${B_m}$ 为噪声带宽，不失一般性地，令 $\sigma _{\rm{s}}^2 = $ $ \sigma _i^2{\rm{ = }}{\sigma ^2}$ ； ${{{h}}_m} \in {\mathbb{C}^{N \times 1}}$ ( $m \in \{ {\rm{s}},i\} $ )为卫星下行链路的信道矢量，根据文献[22]可将其建模为：

${\;\;\;\;\;\;\;\;\;\;\;{{h}}_m} = \sqrt {{G_m}} {{r}}_m^{ - \tfrac{1}{2}} \odot {{b}}_m^{\tfrac{1}{2}} \odot {{\tilde{ h}}_m},m \in \left\{ {{\rm{s}},i} \right\}$

(3)

式中： ${G_m}$ 为地面用户（包括合法用户和窃听者）的抛物面天线增益，具体数值参见文献[23]设置； ${{{r}}_m}$ 为降雨衰落系数矢量； ${{{b}}_m}$ 为卫星波束增益矢量，具体的元素数值参见文献[24]； ${{\tilde{ h}}_m}$ 为信道响应矢量，具体公式参见文献[24]。

由于在实际通信过程中，几乎无法获取理想的信道状态信息，本文假设仅已知统计信道信息，即 ${{{\varPhi }}_{\rm{s}}} = {\rm{E}}\left\{ {{{{h}}_{\rm{s}}}{{h}}_{\rm{s}}^{\rm{H}}} \right\}$ ， ${{{\varPhi }}_i} = {\rm{E}}\left\{ {{{{h}}_i}{{h}}_i^{\rm{H}}} \right\}$ 。根据式（1）和（2），合法用户和第i个窃听者的接收信噪比（signal to noise ratio，SNR）分别表示为：

${\gamma _{\rm{s}}} = \frac{{{{{w}}^{\rm{H}}}{{{\varPhi }}_{\rm{s}}}{{w}}}}{{{\sigma ^2}}}$

(4)

${\gamma _i} = \frac{{{{{w}}^{\rm{H}}}{{{\varPhi }}_i}{{w}}}}{{{\sigma ^2}}}$

(5)

在该系统中，K个窃听者试图合作窃取私密信息。因此，系统的安全速率表示为^[20]：

$\begin{aligned}[b] {R_{\rm{s}}} =& {\left[ {{\rm{lb}}\left( {1 + {\gamma _{\rm{s}}}} \right) - {\rm{lb}}\left( {1 + \sum\limits_{i = 1}^K {{\gamma _i}} } \right)} \right]^ + } = \\ &{\rm{ lb}}\left(1 + \frac{{{{{w}}^{\rm{H}}}{{{\varPhi }}_{\rm{s}}}{{w}}}}{{{\sigma ^2}}}\right) - {\rm{lb}}\left(1 + \frac{{\displaystyle\sum\limits_{i = 1}^K {{{{w}}^{\rm{H}}}{{{\varPhi }}_i}{{w}}} }}{{{\sigma ^2}}}\right){\rm{ = }} \\ &{\rm{ lb}}\left(\frac{{{{{w}}^{\rm{H}}}{{{\varPhi }}_{\rm{s}}}{{w}} + {\sigma ^2}}}{{{{{w}}^{\rm{H}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}{{w}} + {\sigma ^2}} }}\right) \end{aligned} $

(6)

式中， $ {\left[ x \right]^ + } = \max \left( {x,0} \right) $ 。

基于上述系统模型，下文分别以卫星总发射功率最小化和系统安全速率最大化为目标对BF权矢量 ${{w}}$ 进行优化设计。

2 基于卫星发射功率最小化准则的BF方案

考虑到星上资源的有限性，在保证卫星通信系统安全性的同时，面临能量稀缺的问题。因此，需要研究安全速率约束下的卫星发射功率最小化问题，该优化问题在数学上表示为：

$\begin{array}{l} \mathop {\min }\limits_{{w}} \;\;{{{w}}^{\rm{H}}}{{w}} \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;{R_{\rm{s}}} \ge {R_{{\rm{th}}}} \\ \end{array} $

(7)

式中， $ R_{\rm th} $ 为安全速率门限值。

将式（6）代入式（7），优化问题转换为：

${\;\;\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\min }\limits_{{w}} \;\;\;\;\;\;\;\;\;\;{{{w}}^{\rm{H}}}{{w}} \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\dfrac{{{{{w}}^{\rm{H}}}{{{\varPhi }}_{\rm{s}}}{{w}} + {\sigma ^2}}}{{{{{w}}^{\rm{H}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}{{w}} + {\sigma ^2}} }} \ge {2^{{R_{{\rm{th}}}}}} \\ \end{array}} $

(8)

进一步将该问题写为：

${\;\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\min }\limits_{{w}} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{{w}}^{\rm{H}}}{{w}} \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;{{{w}}^{\rm{H}}}\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}}\right) {{w}} \ge ({2^{{R_{{\rm{th}}}}}}{\rm{ - }}1){\sigma ^2} \end{array}} $

(9)

需要注意的是， ${R_{{\rm{th}}}}$ 的取值要保证 $\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}}\right) $ 为半正定矩阵，否则，该优化问题无可行解。由于式（9）中的 ${{{w}}^{\rm{H}}}\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}} \right){{w}}$ 关于 ${{w}}$ 是单调递增的，因此，其不等式约束条件在取等号时取得最优解^[25]，即：

${\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\min }\limits_{{w}} \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{{{w}}^{\rm{H}}}{{w}} \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;{{{w}}^{\rm{H}}}\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}}\right) {{w}}{\rm{ = }}({2^{{R_{{\rm{th}}}}}}{\rm{ - }}1){\sigma ^2} \\ \end{array}} $

(10)

针对问题（10），由拉格朗日乘数法得：

$\begin{array}{l} L\left( {{{w}},\lambda } \right) \buildrel \Delta \over = {{{w}}^{\rm{H}}}{{w}} - \lambda \left( {{{{w}}^{\rm{H}}}\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}}\right) {{w}} - } {\;({2^{{R_{{\rm{th}}}}}}{\rm{ - }}1){\sigma ^2}} \right) \end{array} $

(11)

式中， $\lambda $ 为拉格朗日乘数。通过式（11）对 ${{{w}}^{\rm{H}}}$ 求偏导得：

${\;\;\;\;\;\;\;\;\;\;\;\frac{{\partial L\left( {{{w}},\lambda } \right)}}{{\partial {{{w}}^{\rm{H}}}}}}{\rm{ = }}{{w}}{\rm{ - }}\lambda \left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\sum\limits_{i = 1}^K {{{{\varPhi }}_i}}\right) {{w}}$

(12)

令 $\dfrac{{\partial L\left( {{{w}},\lambda } \right)}}{{\partial {{{w}}^{\rm{H}}}}}{\rm{ = }}0$ ，可得：

$\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\sum\limits_{i = 1}^K {{{{\varPhi }}_i}}\right) {{w}}{\rm{ = }}\frac{1}{\lambda }{{w}}$

(13)

由式（13）可以看出， ${{w}}$ 为 ${{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}}$ 的一个特征向量， $1/\lambda $ 为其对应的特征值。将式（13）两边同时乘以 $\lambda {{{w}}^{\rm{H}}}$ 得：

${\;\;\;\;\;\;\;\;\;\;\;\;\;{{w}}^{\rm{H}}}{{w}}{\rm{ = }}\lambda {{{w}}^{\rm{H}}}\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\sum\limits_{i = 1}^K {{{{\varPhi }}_i}}\right) {{w}}$

(14)

结合问题（10）中的约束条件，式（14）可进一步表示为：

${{{w}}^{\rm{H}}}{{w}}{\rm{ = }}\lambda ({2^{{R_{{\rm{th}}}}}}{\rm{ - }}1){\sigma ^2}$

(15)

将式（15）代入式（10），优化问题（10）的目标函数等同于：

$\mathop {\min }\limits \;\lambda $

(16)

因此，优化问题（10）的最优BF权矢量 ${{{w}}_{{\rm{opt}}}}$ 表示为：

${{{w}}_{{\rm{opt}}}} = \mu {{q}}$

(17)

式中： ${{q}}$ 为矩阵 $\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}}\right) $ 最大特征值对应的归一化特征向量； $\;\mu $ 为满足问题（10）约束条件的系数，表达式为：

${\;\;\;\;\;\;\;\;\;\mu} = \sqrt {\frac{{({2^{{R_{{\rm{th}}}}}}{\rm{ - }}1){\sigma ^2}}}{{\;\;{{{q}}^{\rm{H}}}\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}} \right)\;{{q}}\;}}} $

(18)

最优BF权矢量 ${{{w}}_{{\rm{opt}}}}$ 具体表示为：

$\begin{aligned}[b] {\;\;\;\;\;\;\;\;{{w}}_{{\rm{opt}}}}&{\rm{ = }}\sqrt {\dfrac{{({2^{{R_{{\rm{th}}}}}}{\rm{ - }}1){\sigma ^2}}}{{{{{q}}^{\rm{H}}}\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}} \right)\;{{q}}}}} \cdot \\ &\left\{ {\max .\;{\rm{eigenvector}}\left({{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}}\right) } \right\} \end{aligned} $

(19)

式中， $\max .\;{\rm{eigenvector( }} \cdot {\rm{ }})$ 表示求取矩阵的最大特征值所对应的归一化特征向量。

进一步，得到卫星总发射功率的最小值为：

${\;\;\;\;\;\;\;\;\;\;P_{\rm{T}}^{\min }}\left( {{R_{{\rm{th}}}}} \right) = \frac{{({2^{{R_{{\rm{th}}}}}}{\rm{ - }}1){\sigma ^2}}}{{{\lambda _{\max }}\left( {{{{\varPhi }}_{\rm{s}}}{\rm{ - }}{2^{{R_{{\rm{th}}}}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}} } \right)}}$

(20)

式中， ${\lambda _{\max }}({\rm{ }} \cdot {\rm{ }})$ 表示求取矩阵的最大特征值。

3 基于安全速率最大化准则的BF方案

为了提升卫星与合法用户通信的安全性，研究不同功率约束下的安全速率最大化问题，并分别提出相应的BF算法。

3.1 以总发射功率受限为约束条件

以卫星总发射功率受限为约束条件、安全速率最大化为目标函数建立优化问题，其数学表示为：

$\begin{array}{l} \mathop {\max }\limits_{{w}} \;\;{R_{\rm{s}}} \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;{{{w}}^{\rm{H}}}{{w}} \le {P_{\max }} \end{array} $

(21)

式中， ${P_{\max }}$ 为卫星总发射功率最大值。结合式（6），该优化问题可表示为：

${\;\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\max }\limits_{{w}} \;\;{\rm{lb}}\left(\dfrac{{{{{w}}^{\rm{H}}}{{{\varPhi }}_{\rm{s}}}{{w}} + {\sigma ^2}}}{{{{{w}}^{\rm{H}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}{{w}} + {\sigma ^2}} }}\right) \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\;{{{w}}^{\rm{H}}}{{w}} \le {P_{\max }} \end{array}} $

(22)

令 ${{w}} = \sqrt P {\tilde{ w}}$ 且 ${{\tilde{ w}}^{\rm{H}}}{\tilde{ w}}{\rm{ = }}1$ ，由于 ${\rm{lb}}\left( {{\rm{ }} \cdot {\rm{ }}} \right)$ 为单调递增函数，式（22）可等价转换为：

${\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\max }\limits_{{w}} \;\;\dfrac{{P{{{\tilde{ w}}}^{\rm{H}}}{{{\varPhi }}_{\rm{s}}}{\tilde{ w}} + {\sigma ^2}}}{{P{{{\tilde{ w}}}^{\rm{H}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}{\tilde{ w}} + {\sigma ^2}} }} \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\;{{{\tilde{ w}}}^{\rm{H}}}{\tilde{ w}}{\rm{ = }}1,P \le {P_{\max }} \\ \end{array}} $

(23)

可以看出，式（23）的目标函数是关于P的单调递增函数，在 $P{\rm{ = }}{P_{\max }}$ 时取得最大值。则优化问题进一步转换为：

$ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\max }\limits_{{w}} \;\;\dfrac{{{{{\tilde{ w}}}^{\rm{H}}}({P_{\max }}{{{\varPhi }}_{\rm{s}}}{\rm{ + }}{\sigma ^2}{{I}}){\tilde{ w}}}}{{{\tilde{ w}}}^{\rm{H}}\left(\displaystyle\sum\limits_{i = 1}^K {P_{\max }}{{{\varPhi }}_i}{\rm{ + }}{\sigma ^2}{{I}}\right){\tilde{ w}}} \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\;{{{\tilde{ w}}}^{\rm{H}}}{\tilde{ w}}{\rm{ = }}1 \end{array}}$

(24)

根据广义瑞利熵的性质可知：式（24）中的目标函数的最大值是由矩阵 $\left(\displaystyle\sum\limits_{i = 1}^K {P_{\max }}{{{\varPhi }}_i}{\rm{ + }}{\sigma ^2}{{I}}\right)^{{\rm{ - }}1} ({P_{\max }}{{{\varPhi }}_{\rm{s}}}{\rm{ + }}{\sigma ^2}{{I}})$ 的最大特征值所限定。因此，可得最优BF权矢量 ${{{w}}_{{\rm{opt}}}} = \sqrt {{P_{\max }}} {{\tilde{ w}}_{{\rm{opt}}}}$ 。其中， ${{\tilde{ w}}_{{\rm{opt}}}}$ 为矩阵 $\left(\displaystyle\sum\limits_{i = 1}^K {P_{\max }}{{{\varPhi }}_i}{\rm{ + }}{\sigma ^2}{{I}}\right)^{{\rm{ - }}1} \cdot $ $ ({P_{\max }}{{{\varPhi }}_{\rm{s}}}{\rm{ + }}{\sigma ^2}{{I}})$ 的最大特征值对应的归一化特征向量。

则最优BF矢量 ${{{w}}_{{\rm{opt}}}}$ 的具体表示为：

$\begin{aligned}[b] {{{w}}_{{\rm{opt}}}} = &\sqrt {{P_{\max }}} \;\cdot \\ &\Bigg\{ { {\max .\;{\rm{eigenvector}}\left( {\left(\displaystyle\sum\limits_{i = 1}^K {P_{\max }}{{{\varPhi }}_i}{\rm{ + }}{\sigma ^2}{{I}}\right)^{{\rm{ - }}1} \cdot } \right.} \Bigg.} \Bigg.\\ &\Bigg. {\Bigg. {({P_{\max }}{{{\varPhi }}_{\rm{s}}}{\rm{ + }}{\sigma ^2}{{I}})} \Bigg)} \Bigg\}\\[-8pt] \end{aligned}$

(25)

进一步，得到系统安全速率的最大值为：

$\begin{aligned}[b] {\;\;\;\;\;\;\;\;\;\;R_{\max }}\left( {{P_{\max }}} \right) =& {\rm{lb}}\Bigg( {\lambda _{\max }}\Bigg\{ {\left(\displaystyle\sum\limits_{i = 1}^K {P_{\max }}{{{\varPhi }}_i}{\rm{ + }}{\sigma ^2}{{I}}\right)^{{\rm{ - }}1}\cdot } \Bigg.\Bigg.\\ &\Bigg. {\Bigg. {({P_{\max }}{{{\varPhi }}_{\rm{s}}}{\rm{ + }}{\sigma ^2}{{I}})} \Bigg\}} \Bigg) \end{aligned}$

(26)

3.2 以单波束发射功率受限为约束条件

在多波束卫星通信系统中，每个馈源配置专用的射频链路和功率放大器，每个馈源都存在相应的功率限制。因此，有必要研究单波束发射功率受限情况下的安全速率最大化问题，该优化问题的数学表示为：

${\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\max }\limits_{{w}} \;\;{R_{\rm{s}}} \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\left[ {{W}} \right]{_{n,n}} \le {P_n},\; n{\rm{ = }}1,2, \cdots ,N \end{array}} $

(27)

式中， ${{W}} = {{w}}{{{w}}^{\rm{H}}}$ ， $\left[ {{W}} \right]{_{n,n}}$ 为第n个波束的发射功率， ${P_n}$ 为第n个波束发射功率最大值。结合式（6），优化问题（27）转换为：

${\;\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\max }\limits_{{w}} \;\;{\rm{lb}}\left(\dfrac{{{{{w}}^{\rm{H}}}{{{\varPhi }}_{\rm{s}}}{{w}} + {\sigma ^2}}}{{{{{w}}^{\rm{H}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}{{w}} + {\sigma ^2}} }}\right) \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\left[ {{W}} \right]{_{n,n}} \le {P_n},\; n{\rm{ = }}1,2, \cdots ,N; \\ \;\;\;\;\;\;\;\;\;\;{\rm{rank}}\;{{W}} = 1,\;{{W}} \succeq 0 \end{array}} $

(28)

式中， ${{W}} \succeq 0 $ 表示 ${{W}}$ 为半正定矩阵。由于 ${\rm{lb}}\left( {{\rm{ }} \cdot {\rm{ }}} \right)$ 具有单调递增特性，优化问题（28）可进一步写为：

${\;\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\max }\limits_{{w}} \;\dfrac{{{{{w}}^{\rm{H}}}{{{\varPhi }}_{\rm{s}}}{{w}} + {\sigma ^2}}}{{{{{w}}^{\rm{H}}}\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}{{w}} + {\sigma ^2}} }} \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\left[ {{W}} \right]{_{n,n}} \le {P_n},\; n{\rm{ = }}1,2, \cdots ,N; \\ \;\;\;\;\;\;\;\;\;\;\;{\rm{rank}}\;{{W}} = 1,\;{{W}}\succeq 0 \end{array}} $

(29)

可等价转换为：

${\;\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\max }\limits_{{W}} \;\dfrac{{{\rm{tr}}({{{\varPhi }}_{\rm{s}}}{{W}}) + {\sigma ^2}}}{{\rm tr}\left(\displaystyle\sum\limits_{i = 1}^K {{{\varPhi }}_i}{{W}}\right)+ {\sigma ^2} } \\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;\left[ {{W}} \right]{_{n,n}} \le {P_n},\; n{\rm{ = }}1,2, \cdots ,N; \\ \;\;\;\;\;\;\;\;\;\;{\rm{rank}}\;{{W}} = 1,\;{{W}}\succeq 0 \end{array}} $

(30)

式中， ${\rm{tr}}\left( {{\rm{ }} \cdot {\rm{ }}} \right)$ 表示矩阵的迹。由于问题（30）的最后一个约束条件是非凸形式的，无法进行直接求解。因此，忽略 ${{W}}$ 的秩1约束，通过引入辅助变量 $t$ 将非凸优化问题转为凸优化问题，表示为：

${\;\;\;\;\;\;\;\;\;\;\;\begin{array}{l} \mathop {\max }\limits_{{{W}},t} \;\;\;\;\;t\\ \;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\;{\rm{tr}}\left(\left({{{\varPhi }}_{\rm{s}}} - t\displaystyle\sum\limits_{i = 1}^K {{{{\varPhi }}_i}} \right){{W}}\right) \ge (t - 1){\sigma ^2};\\ \;\;\;\;\;\;\;\;\;\; {\rm{tr}}({{W}}{{{X}}_n})\; \le 1,\;n{\rm{ = }}1,2, \cdots ,N;\\ \;\;\;\;\;\;\;\;\;\;{{W}}\succeq 0,\;t \ge 0 \end{array}}$

(31)

式中， ${{{X}}_n}$ 表示主对角线元素为 $1/{P_n}$ 而其余元素全为0的矩阵。问题（31）是对t的1维搜索问题，提出基于黄金分割法的BF算法求得最优BF权矢量 ${{{w}}_{{\rm{opt}}}}$ ，具体步骤如下：

1）输入精度 $\varepsilon $ ，初始搜索区间 $a\;{\rm{ = }}\;0 $ ， $ b\;{\rm{ = }}\;f\left( {{P_{\max }}} \right) $ ，其中， ${P_{\max }} = \displaystyle\sum\limits_{n = 1}^N {{P_n}} $ ， $f(t) = \dfrac{{{\rm{tr}}({{{\varPhi }}_{\rm{s}}}{{W}}) + {\sigma ^2}}}{{\rm tr}\left(\displaystyle\sum\limits_{i = 1}^K {{{\varPhi }}_i}{{W}}\right)+ {\sigma ^2}} $ 。

2）计算 ${t_1} = a + 0.382(b - a)$ ， ${t_2} = a + 0.618(b - a)$ 。

3）将 ${t_1}$ 和 ${t_2}$ 代入到问题（31）中利用MATLAB软件的CVX凸优化工具包求解 $f\left( {{t_1}} \right)$ 、 $f\left( {{t_2}} \right)$ 。

4）比较 $f\left( {{t_1}} \right)$ 和 $f\left( {{t_2}} \right)$ 的大小：如果 $f\left( {{t_1}} \right) < f\left( {{t_2}} \right)$ ，令 $a = {t_1}$ ， ${t_1}= {t_2}$ ， $f\left( {{t_1}} \right) = f\left( {{t_2}} \right)$ ， ${t_2} = a + 0.618(b - a)$ ，将 ${t_2}$ 代入到问题（31）中利用CVX凸优化工具包求解 $f\left( {{t_2}} \right)$ ；否则，令 $b = {t_2}$ ， ${t_2} = {t_1}$ ， $f\left( {{t_2}} \right) = f\left( {{t_1}} \right)$ ， ${t_1} = a + 0.382(b - a)$ ，将 ${t_1}$ 代入到问题（31）中利用CVX凸优化工具包求解 $f\left( {{t_1}} \right)$ 。

5）如果 $\;\left( {b - a} \right)\; < \varepsilon $ ，停止迭代，执行下一步；否则，返回步骤4）继续迭代。

6）令 ${t^*} = (a + b)/2$ ，将 ${t^*}$ 代入到问题（31）中利用CVX凸优化工具包得到 ${{{W}}^*}$ 。

7）对 ${{{W}}^*}$ 特征值分解，得到最优BF权矢量 ${{{w}}_{{\rm{opt}}}}$ 。

该算法的计算复杂度与优化变量的个数、线性矩阵不等式（linear matrix inequality，LMI）约束的个数和大小有关。从问题（31）可以看出，该问题中包括 $N^2$ 个设计变量和1个辅助变量、1个大小为 $N$ 的LMI约束、 $N+2$ 个大小为1的LMI约束。因此，该算法的计算复杂度为 $O\left( {\sqrt {2N + 2} \cdot m \cdot } \right.\left[ {{N^2}\left( {N + m} \right)} \right. + \left( {N{\rm{ + }}2} \right) \cdot \left( {1{\rm{ + }}m} \right)+ $ $ \left. {\left. { {m^2}} \right]} \right),\;m = O\left( {{N^2} + 1} \right)$ 。

4 仿真与分析

下面通过计算机仿真验证本文所提BF方案的可行性与有效性，分析不同参数对系统安全性能的影响。假设已知合法用户和窃听者的地理位置（ $x,y$ ）（单位为km）分别为：合法用户（0 km，0 km）、窃听者1（100 km，–300 km）、窃听者2（100 km，300 km）。仿真主要参数如表1所示。

4.1 基于卫星发射功率最小化准则的BF方案

图2为卫星发射功率最小化准则下波束数为7时的归一化信噪比分布。从图2中可知，该方案计算出的BF权矢量的主瓣对准了合法用户，同时，在所有窃听者处产生了至少–30 dB的零陷，证明了所提BF方案可以有效抑制窃听过程，保障了卫星与合法用户间的安全通信。

图3为波束数分别为3、4、7的情况下，采用不同BF方案所得卫星总发射功率随安全速率变化的曲线。从图3中可以看出，所提BF方案与迫零（zero force，ZF）BF方案的卫星总发射功率均随着安全速率门限值的增大而增加，同时，随着波束数的增加而减小。这是因为当安全速率门限值增大时，卫星需要消耗更多的发射功率以保证卫星与合法用户间通信的安全性；当波束数增加时，合法用户获得更高的天线增益，因此在安全速率门限值一定的情况下，波束数的增加可有效地减少卫星发射功率的消耗。同时，由图3可以看出所提BF方案的性能要优于ZF BF方案，这是因为与所提BF方案相比，ZF BF方案需要消耗更多的发射功率来抑制窃听者的窃听。

4.2 基于安全速率最大化准则的BF方案 4.2.1 以总发射功率受限为约束条件

图4为总发射功率约束下波束数为7时的归一化信噪比分布。

从图4中可以看出，合法用户的接收信噪比在0 dB左右，两个窃听者的接收信噪比在–30 dB以下，验证了该方案可以有效抑制窃听过程。

图5为总发射功率约束下安全速率与波束数量的关系直方图。从图5中可看出：安全速率随波束数的增加而增大；当波束数固定时，安全速率随着总发射功率门限值的增大而增长。因此，在总发射功率门限值一定的情况下，可以通过增加波束数以提高系统的安全速率。

图6为波束数为7时不同BF方案下的安全速率随总发射功率变化的曲线。从图6中可以看出：所提BF方案的性能要优于ZF BF方案和最大比发射（maximal ratio transmission，MRT）BF方案；并且，随着总发射功率的增大，ZF BF方案的性能与所提BF方案性能不断接近，而MRT BF方案的性能逐渐趋于稳定。

4.2.2 以单波束发射功率受限为约束条件

图7为单波束发射功率约束下波束数为7时的归一化信噪比分布。由图7可知，所提方案能够在所有窃听者处产生至少–30 dB的零陷，表明了该方案可以有效抑制窃听过程。

图8为单波束发射功率约束下的安全速率与波束数量的关系直方图。由图8可以看出，增加波束数可以提高系统安全速率，这是因为波束数的增加能够提高用户的天线增益，故提高了系统性能。

图9为波束数为7时总发射功率受限情况与单波束功率受限情况的对比。

从图9中可以看出，单波束功率约束下的安全速率小于总功率约束下的安全速率，这是因为总功率约束条件是单波束功率约束条件的放松形式，因此总功率约束条件下的方案可以获得更好的安全性能。但考虑到总发射功率受限情况与单波束功率受限情况的性能存在一定的差距，为了更贴近多波束卫星通信系统的实际情况，应采用单波束功率约束下的BF方案进行安全波束成形设计。

5 结　论

本文研究多窃听场景下多波束卫星系统的安全传输问题，在已知统计信道状态信息条件下，根据不同的性能准则提出了多种BF方案。首先，针对安全速率约束下的卫星发射功率最小化问题，提出了基于拉格朗日乘数法的BF方案并推导出最优BF权矢量的解析表达式。其次，针对卫星总发射功率和单波束发射功率受限条件下的安全速率最大化问题，分别提出了基于广义瑞利熵和黄金分割法的BF算法，并推导出相应的最优BF权矢量。最后，计算机仿真验证了所提BF方案的安全性能均优于现有的ZF BF方案和MRT BF方案，所提的BF方案均可以通过增加波束数以提高系统的安全性能。在后续的工作中，将进一步研究在非合作的窃听模式下，多波束卫星通信系统的安全传输问题。