CN104020439A - 基于空间平滑协方差矩阵稀疏表示的波达方向角估计方法 - Google Patents

Abstract

本发明公开了一种基于空间平滑协方差矩阵稀疏表示的波达方向角估计方法，主要解决现有技术运算量大，处理相干信号源性能差，造成无源定位估计误差大的问题，其实现步骤是：1)采用天线接收机形成均匀线阵；2)采用空间平滑技术计算阵列输出的空间平滑协方差矩阵；3)将空间平滑协方差矩阵矢量化，得到稀疏模型向量；4)将空域网格划分，构造超完备基；5)根据稀疏模型向量和超完备基的稀疏表示关系，建立约束优化方程；6)采用凸优化方法求解约束优化方程得到最优估计；7)根据最优估计值绘制幅度谱图，获得波达方向角度值。本发明提高了无源测向的运算速度及低信噪比下对相干信号源的估计性能，可用于目标侦察和无源定位。

Description

基于空间平滑协方差矩阵稀疏表示的波达方向角估计方法

技术领域

本发明属于信号处理技术领域，特别涉及一种基于均匀线阵的波达方向角估计方法，可用于目标侦察与无源定位。

背景技术

波达方向角DOA估计是利用处于空间不同位置的天线阵列接收多个不同方向的信号源发出的信号，运用现代信号处理方法快速准确的估计出信号源的方向，在雷达、声纳、无线通信等领域具有重要应用价值。针对该问题的研究中，出现较早、应用较为广泛的是多重信号分类MUSIC等基于子空间的模型，之后的大部分DOA估计算法都是利用该模型生成的。近年来，由Donoho等提出的压缩感知理论为DOA估计问题提供了一种新思路，从而产生出一类基于稀疏表示模型的DOA估计算法。

目前，基于稀疏信号表示的DOA估计方法中最经典的是L1_SVD方法。L1_SVD方法利用阵列接收数据奇异值分解得到的信号子空间构造稀疏表示模型，然后通过二阶锥规划对L1范数约束模型进行求解。近年来又出现了基于阵列协方差向量稀疏表示的L1_SRACV算法，联合稀疏逼近JLZA算法等，但这些算法均运算量较大，且在低信噪比情况下，角度分辨率不理想，对抗相干信号性能不强。而在实际应用中，目标侦察与无源定位均需要在角度估计的基础上进行，以上算法中的缺陷将造成目标侦察和无源定位反应速度慢和估计误差较大的不足。

发明内容

本发明的目的在于针对上述已有技术的不足，提出一种基于空间平滑协方差矩阵稀疏表示的波达方向角度估计方法，以在降低运算量的情况下，提高目标侦察和无源定位在低信噪比下的检测成功率和对相干信号的估计能力，避免因角度估计误差引起的目标侦察失误。

为实现上述目的，本发明的实现步骤包括如下：

1)采用M个天线接收机形成均匀线性阵列，并假设有K个信号入射到该均匀线性阵列，各天线接收机间距均为d，每个天线接收机称为一个阵元，其中，M≥2，K≥1，0＜d≤λ/2，λ为入射窄带信号波长；

2)利用均匀线性阵列的平移不变性，采用前后向空间平滑方法计算均匀线阵输出的空间平滑协方差矩阵R^fb；

3)根据空间平滑协方差矩阵R^fb构造稀疏模型向量y：

y＝vec(R^fb)，

其中，vec(·)表示向量化运算；

4)构造超完备基Φ并定义一个空域稀疏向量u：

4a)根据入射信号的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个区域，定义为波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，q＝1,2,…,Q，Q＞＞M，构造一个信号稀疏化后对应的(2m-1)×Q维导向矩阵B(θ)：

B(θ)＝[b(θ₁),…,b(θ_q),…,b(θ_Q)]，

其中，b(θ_q)表示角度θ_q对应的导向向量：

$b\left({\mathrm{\θ}}_q\right)={[e^{(m-1)\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q},...,e^{\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q},1,e^{-\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q},...,e^{-(m-1)\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q}]}^T,$

其中，表示两个相邻阵元间的相位差，(·)^T表示矩阵转置运算；

4b)计算选择矩阵G：

$G=[\mathrm{vec}\left(J_{m-1}\right),...,\mathrm{vec}\left(J_1\right),\mathrm{vec}\left(J_0\right),\mathrm{vec}\left(J_1^T\right),...,\mathrm{vec}\left(J_{m-1}^T\right)],$

其中，J₀,J₁,…,J_m-1按下式计算：

$J_1=\left[\begin{array}{cc}{0}_{m-l,l}& I_{m-l}\\ 0_{l,l}& 0_{l,m-l}\end{array}\right],l=\mathrm{0,1},...,m-1;$

4c)根据选择矩阵G和导向矩阵B(θ)，得到超完备基Φ：

Φ＝GB(θ)；

4d)定义一个Q×1维的空域稀疏向量u：

u＝[u₁,u₂…,u_q,…,u_Q]^T，

其中，u_q为u中第q个元素，u中元素均为未知变量；

5)通过稀疏重构获得空域稀疏向量u的最优估计

5a)根据空间平滑协方差矩阵R^fb计算加权矩阵W：

$W=1/N{\left(R^{\mathrm{fb}}\right)}^T\⊗R^{\mathrm{fb}},$

其中，N表示采样点数，表示矩阵的Kronecker积运算；

5b)利用稀疏表示的思想，根据步骤(3)和(4)得到的结果将波达方向角估计问题转化为求解约束优化方程式：

$\underset{\hat{u}}{\min }{\left|\right|\hat{u}\left|\right|}_{1}s.t.{\left|\right|W^{-1/2}(y-\mathrm{\Φ}\hat{u}-{\mathrm{\σ}}_n^2\mathrm{vec}\left(I_m\right))\left|\right|}_2\≤\mathrm{\ϵ}$

其中||·||₁和||·||₂分别表示求向量l₁范数和l₂范数，s.t.表示约束关系，Φ为超完备基，为噪声方差，I_m表示m×m维单位矩阵，ε表示误差的允许值， 表示自由度为m²、置信度为ρ的卡方分布，ρ的取值为0.001；

5c)利用凸优化方法求解上述约束优化方程，得到空域稀疏向量的最优估计

6)以波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]的值为x轴坐标，以稀疏向量的最优解的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

本发明与现有技术相比具有以下优点：

1)本发明采用稀疏表示的思想将波达方向角估计问题转化为稀疏重构问题，是新理论技术与传统问题的结合，利用入射信号源的空域稀疏特性进行建模，避免了传统算法的角度搜索或角度匹配过程，同时突破了阵列分辨率的瑞丽限，提高了测角精度。

2)本发明构建了矢量化协方差矩阵的稀疏表示模型，该构建将多测量矢量MMV问题转化为单测量矢量SMV模型，在稀疏重构过程中大大降低了运算量。

3)本发明采用前后向空间平滑技术得到阵列输出的协方差矩阵，兼顾相干和非相干信号源入射角的分辨，尤其对相干信号源有良好的角度估计性能，在现实环境中具有更实际的应用价值。

附图说明

图1是本发明的实现流程图；

图2是本发明与现有两种波达方向角估计方法运算时间对比图；

图3是本发明与现有两种波达方向角估计方法在不同信噪比条件下处理非相干信号的发现概率对比图；

图4是本发明与现有两种波达方向角估计方法在不同信噪比条件下处理相干信号的发现概率对比图；

具体实施方式

以下参照附图，对本发明的技术方案和效果作进一步的详细说明。

参照图1，本发明的实现步骤如下：

步骤1：利用天线接收机形成均匀线阵。

每隔间距d放置1个天线接收机，共放置M个，形成一个均匀线性天线阵列，每个天线接收机称为一个阵元，假设有K个远场窄带信号入射到该均匀线阵上，且信号在传播过程中加入了均值为零的复高斯白噪声，其中，M≥2，K≥1，0＜d≤λ/2，λ为入射窄带信号波长。

步骤2：采用前后向空间平滑方法计算天线阵列输出的空间平滑协方差矩阵R^fb。

2a)利用均匀线阵的平移不变性，将整个均匀线阵划分成相互交错的P个子阵，每个子阵的阵元数为m，即以均匀线阵中第1个到第m个阵元构成第一个子阵，第2个到第m+1个阵元构成第二个子阵，以此类推，第P个阵元到第M个阵元构成最后一个子阵，将这P个子阵依次排列，称为前向子阵，其中，K≤m≤M，P＝M-m+1；

2b)计算前向平滑协方差矩阵R^f：

计算第i个前向子阵的协方差矩阵

$R_i^f=E\left\{X_i^f\left(t\right){\left(X_i^f\left(t\right)\right)}^H\right\},$

其中，i＝1,2,…,P，上标f表示前向，表示前向子阵中第i个子阵的输出信号数据，E(·)表示求数学期望，(·)^H表示共轭转置运算，对所有前向子阵的协方差矩阵求平均，得到前向平滑协方差矩阵R^f：

$R^f=\frac{1}{P}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{R}_i^f;$

2c)以均匀线阵中第M个到第M-m+1个阵元构成第一个子阵，第M-1个到第M-m个构成第二个子阵，以此类推，第m个到第1个阵元构成最后一个子阵，将这P个子阵依次排列，称为后向子阵；

2d)计算后向平滑协方差矩阵R^b：

计算第i个后向子阵的协方差矩阵

$R_i^b=E\left\{X_i\left(t\right){\left(X_i^b\left(t\right)\right)}^H\right\},$

其中，i＝1,2,…,P，上标b表示后向，表示后向子阵中第i个子阵的输出信号数据，对所有后向子阵的协方差矩阵求平均，得到后向平滑协方差矩阵R^b：

$R^b=\frac{1}{P}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{R}_i^b;$

2e)对前向平滑协方差矩阵R^f和后向平滑协方差矩阵R^b求平均，得到空间平滑协方差矩阵R^fb：

R^fb＝1/2(R^f+R^b)。

步骤3：对空间平滑协方差矩阵R^fb进行矢量化运算，得到稀疏模型向量y：

y＝vec(R^fb)，

其中，vec(·)表示向量化运算。

步骤4：构造超完备基Φ。

根据稀疏信号重构理论，任意信号都可以由一个基矩阵线性表示，在这里，构造超完备基Φ矩阵的目的就是将稀疏重构向量y通过矩阵的形式表示出来，便于构建单测量矢量模型，其构造步骤如下：

4a)根据入射信号源所具有的空域稀疏特性，对观测空域进行空间网格划分处理，即将观测空域[-90°,90°]等间隔划分成Q个区间，θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，θ表示波达方向角范围，θ_q为第q个角度区间，q＝1,2,…,Q，Q＞＞M，网格划分间隔的取值根据期望达到的角度估计精度进行设定，网格划分间隔越小，则最终得到的角度估计值精度越高。

4b)构造一个空域稀疏化后对应的(2m-1)×Q维导向矩阵B(θ)：

B(θ)＝[b(θ₁),…,b(θ_q),…,b(θ_Q)]，

其中，b(θ_q)表示B(θ)的第q列，b(θ_q)为角度θ_q对应的2m-1维导向向量：

$b\left({\mathrm{\θ}}_q\right)={[e^{(m-1)\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q},...,e^{\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q},1,e^{-\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q},...,e^{-(m-1)\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q}]}^T,$

其中，表示两个相邻阵元间的相位差，(·)^T表示矩阵转置运算；

4c)计算选择矩阵G：

$G=[\mathrm{vec}\left(J_{m-1}\right),...,\mathrm{vec}\left(J_1\right),\mathrm{vec}\left(J_0\right),\mathrm{vec}\left(J_1^T\right),...,\mathrm{vec}\left(J_{m-1}^T\right)],$

其中，J₀,J₁,…,J_m-1按下式计算：

$J_1=\left[\begin{array}{cc}{0}_{m-l,l}& I_{m-l}\\ 0_{l,l}& 0_{l,m-l}\end{array}\right],l=\mathrm{0,1},...,m-1$

4d)根据选择矩阵G和导向矩阵B(θ)，得到超完备基Φ：

Φ＝GB(θ)。

步骤5：定义一个空域稀疏向量u。

假设每个角度区间θ_q都对应一个入射信号，这样就构造一个Q×1维的空域稀疏向量u＝[u₁,u₂,…,u_q,…,u_Q]^T，u_q为u中第q个元素，q＝1,2,…,Q，u中元素均为未知变量。

步骤6：通过稀疏重构获得空域稀疏向量u的最优估计

6a)根据空间平滑协方差矩阵R^fb计算加权矩阵W：

$W=1/N{\left(R^{\mathrm{fb}}\right)}^T\⊗R^{\mathrm{fb}},$

其中，N表示采样点数，表示矩阵的Kronecker积运算；

6b)根据稀疏表示的思想，根据步骤(3)和(4)得到的结果将波达方向角估计问题转化为求解以下的约束优化方程式：

$\underset{\hat{u}}{\min }{\left|\right|\hat{u}\left|\right|}_{1}s.t.{\left|\right|W^{-1/2}(y-\mathrm{\Φ}\hat{u}-{\mathrm{\σ}}_n^2\mathrm{vec}\left(I_m\right))\left|\right|}_2\≤\mathrm{\ϵ}$

其中||·||₁和||·||₂分别表示求向量l₁范数和l₂范数，s.t.表示约束关系，Φ为超完备基，为噪声方差，I_m表示m×m维单位阵，ε表示误差的允许值，表示自由度为m²、置信度为ρ的卡方分布，ρ的取值为0.001；

6c)采用凸优化方法求解上述约束优化方程式，凸优化是一种比较特殊的优化，是指目标函数和约束函数均为凸函数的优化问题，凸优化问题有一套非常完备的解决算法，在此采用现有针对凸优化问题的软件包CVX(Grant M,Boyd S.CVX:Matlabsoftware for disciplined convex programming[J].2008[Online]Available:http://stanfordedu/～boyd/cvx)来求解，通过该方法能够快速地得到空域稀疏向量的最优估计

步骤7：根据最优估计绘制幅度谱图，得到波达方向角的估计值。

步骤6得到的最优估计为一个K稀疏向量，即其中只有K个值为非零值，其余值均为零，这K个非零值对应的空间角度区间就是入射信号源的方向，因此，以波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]的值为x轴坐标，以稀疏向量的最优估计的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

本发明的效果可通过以下仿真说明：

1.仿真条件与方法：

采用多个天线接收机形成均匀线阵，天线接收机间距d为半波长，划分子阵个数为3个，信噪比SNR为0dB，采样点数为500，观测空域角度范围为[-90°,90°]，空间网格划分间隔为1°。

仿真方法为本发明和现有的L1_SVD、L1_SRACV算法。

2.仿真内容与结果：

仿真1：假设有两个非相干窄带信号分别以角度16°和26°入射到均匀线阵，均匀线阵的阵元数由6个增加到12个，利用本发明和现有的L1_SVD、L1_SRACV算法分别进行100次独立的波达方向角估计实验，并统计运算时间，结果如图2所示。图2中横坐标表示阵元数，纵坐标表示运算时间。

从图2可以看出，本发明与现有L1_SVD、L1_SRACV算法相比，大幅度降低了角度估计的运算量，即使阵元数增加，运算量仍基本保持稳定不变。

仿真2：设定均匀线阵的阵元数为10，假设有两个非相干窄带信号分别以角度16°和46°入射到均匀线阵，信噪比SNR由-10dB增加到10dB。利用本发明和现有的L1_SVD、L1_SRACV算法分别进行100次独立的波达方向角度估计实验，分别计算不同信噪比下三种方法的发现概率，结果如图3所示，图3中横坐标表示信噪比值，纵坐标表示发现概率。

从图3可以看出，在低信噪比的情况下，本发明的发现概率明显高于其它两个算法。

仿真3：设定均匀线阵的阵元数为10，假设两个相干窄带信号分别以角度16°和46°入射到均匀线阵，信噪比SNR由-10dB增加到10dB。利用本发明和现有的L1_SVD、L1_SRACV算法分别进行100次独立的波达方向角度估计实验，分别计算不同信噪比下的发现概率，结果如图4所示，图4中横坐标表示信噪比值，纵坐标表示发现概率。

从图4可以看出，在低信噪比的情况下，本发明的发现概率明显高于其它两种方法，说明本发明在处理相干信号时有明显的优势。

综上，本发明在降低了角度估计运算量的同时，提高了低信噪比情况下角度估计的发现概率，并且兼具非相干和相干信号源的分辨能力，保证了目标侦察和无源定位的快速反应和准确有效，避免了相干信号源背景下角度估计失误的问题。

Claims (2)

1.一种基于空间平滑协方差矩阵稀疏表示的波达方向角估计方法，包括以下步骤：

1)采用M个天线接收机形成均匀线性阵列，并假设有K个信号入射到该均匀线性阵列，各天线接收机间距均为d，每个天线接收机称为一个阵元，其中，M≥2，K≥1，0＜d≤λ/2，λ为入射窄带信号波长；

2)利用均匀线性阵列的平移不变性，采用前后向空间平滑方法计算均匀线阵输出的空间平滑协方差矩阵R^fb；

3)根据空间平滑协方差矩阵R^fb构造稀疏模型向量y：

y＝vec(R^fb)，

其中，vec(·)表示向量化运算；

4)构造超完备基Φ并定义一个空域稀疏向量u：

4a)根据入射信号的空域稀疏特性，采用空间网格划分方法，将观测空域[-90°,90°]等间隔划分成Q个区域，定义为波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]，q＝1,2,…,Q，Q＞＞M，构造一个信号稀疏化后对应的(2m-1)×Q维导向矩阵B(θ)：

B(θ)＝[b(θ₁),…,b(θ_q),…,b(θ_Q)]，

其中，b(θ_q)表示角度θ_q对应的导向向量：

$b\left({\mathrm{\θ}}_q\right)={[e^{(m-1)\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q},...,e^{\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q},1,e^{-\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q},...,e^{-(m-1)\frac{j2\mathrm{\πd}}{\mathrm{\λ}}\sin {\mathrm{\θ}}_q}]}^T,$

其中，表示两个相邻阵元间的相位差，(·)^T表示矩阵转置运算；

4b)计算选择矩阵G：

$G=[\mathrm{vec}\left(J_{m-1}\right),...,\mathrm{vec}\left(J_1\right),\mathrm{vec}\left(J_0\right),\mathrm{vec}\left(J_1^T\right),...,\mathrm{vec}\left(J_{m-1}^T\right)],$

其中，J₀,J₁,…,J_m-1按下式计算：

$J_1=\left[\begin{array}{cc}{0}_{m-l,l}& I_{m-l}\\ 0_{l,l}& 0_{l,m-l}\end{array}\right],l=\mathrm{0,1},...,m-1;$

4c)根据选择矩阵G和导向矩阵B(θ)，得到超完备基Φ：

Φ＝GB(θ)；

4d)定义一个Q×1维的空域稀疏向量u：

u＝[u₁,u₂…,u_q,…,u_Q]^T，

其中，u_q为u中第q个元素，u中元素均为未知变量；

5)通过稀疏重构获得空域稀疏向量u的最优估计

5a)根据空间平滑协方差矩阵R^fb计算加权矩阵W：

$W=1/N{\left(R^{\mathrm{fb}}\right)}^T\⊗R^{\mathrm{fb}},$

其中，N表示采样点数，表示矩阵的Kronecker积运算；

5b)利用稀疏表示的思想，根据步骤(3)和(4)得到的结果将波达方向角估计问题转化为求解约束优化方程式：

$\underset{\hat{u}}{\min }{\left|\right|\hat{u}\left|\right|}_{1}s.t.{\left|\right|W^{-1/2}(y-\mathrm{\Φ}\hat{u}-{\mathrm{\σ}}_n^2\mathrm{vec}\left(I_m\right))\left|\right|}_2\≤\mathrm{\ϵ}$

其中||·||₁和||·||₂分别表示求向量l₁范数和l₂范数，s.t.表示约束关系，Φ为超完备基，为噪声方差，I_m表示m×m维单位矩阵，ε表示误差的允许值，表示自由度为m²、置信度为ρ的卡方分布，ρ的取值为0.001；

5c)利用凸优化方法求解上述约束优化方程，得到空域稀疏向量的最优估计

6)以波达方向角范围θ＝[θ₁,θ₂,…,θ_q,…,θ_Q]的值为x轴坐标，以稀疏向量的最优解的幅度值为y轴坐标，绘制幅度谱图，从该幅度谱图中按照从高到低的顺序寻找幅值较大的前K个谱峰，这些谱峰的峰值点所对应的x轴坐标即为所求的波达方向角度值。

2.根据权利要求1所述的基于空间平滑协方差矩阵稀疏表示的波达方向角估计方法，其中步骤2)所述的采用前后向空间平滑方法计算均匀线阵输出的空间平滑协方差矩阵R^fb，按如下步骤进行：

2a)将整个均匀线阵划分成相互交错的P个子阵，每个子阵的阵元数为m，即以均匀线阵中第1个到第m个阵元构成第一个子阵，第2个到第m+1个阵元构成第二个子阵，以此类推，第P个阵元到第M个阵元构成最后一个子阵，将这P个子阵依次排列，称为前向子阵，其中，K≤m≤M，P＝M-m+1；

2b)计算前向平滑协方差矩阵R^f：

计算第i个前向子阵的协方差矩阵

$R_i^f=E\left\{X_i^f\left(t\right){\left(X_i^f\left(t\right)\right)}^H\right\},$

其中，i＝1,2,…,P，上标f表示前向，表示前向子阵中第i个子阵的输出信号数据，E(·)表示求数学期望，(·)^H表示共轭转置运算，对所有前向子阵的协方差矩阵求平均，得到前向平滑协方差矩阵R^f：

$R^f=\frac{1}{P}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{R}_i^f$

2c)以均匀线阵中第M个到第M-m+1个阵元构成第一个子阵，第M-1个到第M-m个构成第二个子阵，以此类推，第m个到第1个阵元构成最后一个子阵，将这P个子阵依次排列，称为后向子阵；

2d)计算后向平滑协方差矩阵R^b：

计算第i个后向子阵的协方差矩阵

$R_i^b=E\left\{X_i\left(t\right){\left(X_i^b\left(t\right)\right)}^H\right\},$

其中，i＝1,2,…,P，上标b表示后向，表示后向子阵中第i个子阵的输出信号数据，对所有后向子阵的协方差矩阵求平均，得到后向平滑协方差矩阵R^b：

$R^b=\frac{1}{P}\underset{i=1}{\overset{P}{\mathrm{\Σ}}}{R}_i^b;$

2e)对前向平滑协方差矩阵R^f和后向平滑协方差矩阵R^b求平均，得到空间平滑协方差矩阵R^fb：

R^fb＝1/2(R^f+R^b)。