CN106383349A - 一种基于x波段多普勒雷达降水估测系统及方法 - Google Patents

Abstract

本发明提供一种基于X波段多普勒雷达降水估测系统及方法，系统包括X波段多普勒雷达、雨滴谱仪、监控主机和显示系统，所述X波段多普勒雷达包括外机和内机两个部分。本发明通过X波段多普勒雷达垂直指向探测降水区，测量得到多普勒速度谱廓线，导出降水粒子的滴谱廓线，从而能够计算降水粒子的衰减系数廓线和反射率因子廓线，就可以建立衰减系数k和反射率因子Z之间的关系k＝c*Zd，并将关系应用于60公里范围内的降水，解决X波段多普勒雷达的衰减订正问题，扩大了X波段多普勒雷达在降水方面的应用范围。

Description

一种基于X波段多普勒雷达降水估测系统及方法

技术领域

本发明属于大气科学技术领域，尤其涉及一种基于X波段多普勒雷达降水估测系统及方法。

背景技术

X波段多普勒雷达体积小、重量轻，探测范围一般在60公里半径，广泛应用于中小尺度灾害性天气预警和局地山洪、城市内涝的应急预警。X波段的雷达相对于S波段和C波段的雷达，波长较短，大气中的降水(雨、雪、雹等)粒子对雷达波的衰减不容忽视，尤其是在出现局地强天气时，大粒子的衰减可能影响雷达的正常探测。

通常，常规雷达的衰减订正采用经验关系式k＝c*Zd，其中，Z为雷达测量的降水粒子的反射率因子，k为衰减系数，c和d为常数，在假定雨滴谱为Gamma分布的基础上，理论上，可以计算c和d的近似值。然而，在实际应用中，由于缺乏雨滴谱数据，一般雷达都避免使用衰减波段的雷达测量降水，而是使用如S波段或C波段雷达，或者使用双极化或双波长雷达来测量降水，衰减订正成为制约X波段多普勒雷达应用的重要因素。

雷达之所以能用于降水估测是因为雷达测量的反射率因子Z和雨滴谱有关，而降水也取决于雨滴谱，在一定的假定下，两者之间可以建立简单的关系式Z＝A*Rb，其中，R为雨强，A和b为系数。由于雷达的波束在大气中传播，而降水是落到地面的雨滴，雨滴从空中到地面的过程中，受到风、温、湿、压的影响，可能增大、缩小、合并、破碎，因此，地面降水和空中降水存在时间、空间和强度上的差异，这在降水估测中需要特别关注。

发明内容

为了解决上述技术问题，本发明提供一种基于X波段多普勒雷达降水估测系统及方法，旨在解决现有的基于X波段多普勒雷达降水估测系统及方法存在的衰减订正、应用范围受限、效率低下等的问题。

本发明的目的旨在提供一种基于X波段多普勒雷达的降水估测方法，所述基于X波段多普勒雷达的降水估测方法包括以下步骤：

1)、用无地物的最低仰角作平面位置扫描，得到雷达测量的反射率因子场；最低仰角≤18°；

2)、X波段多普勒雷达作垂直指向探测降水区，判断探测降水区即雷达处是否有降水；

3)当探测降水区有降水时，测量降水粒子的反射率因子和多普勒速度谱廓线，并由降水粒子的多普勒速度谱廓线计算降水粒子的滴谱廓线，转步骤4)；

4)当探测降水区没有降水而某个雨滴谱仪处有降水时，雨滴谱仪得到雨滴谱，转步骤3)；

5)当探测降水区和雨滴谱仪所在位置都没有降水时，用缺省的系数c和d，进行衰减订正，并用缺省的系数A和b进行降水估测，转步骤6)；通过最小二乘法拟合确定系数c和d；所述系统性偏差修正为在雨滴谱仪所在格点周边，计算雨滴谱得到的反射率因子与雷达测量得到的反射率因子的差值，视最小的差值为最佳匹配，这样，每个雨滴谱仪都有一个最小差值，计算这些最小差值的平均值，视此平均值为雷达测量的反射率因子的系统性偏差，并对雷达测量的反射率因子进行订正；

6)计算降水区的衰减系数k、降水粒子的反射率因子Z和雨强R；

7)根据衰减系数k和反射率因子Z，拟合k＝c*Zd确定系数c和d，将此关系应用于降水区域，实现雷达测量的反射率因子的衰减订正，并根据反射率因子Z和雨强R，拟合Z＝A*Rb确定系数A和b；

8)根据衰减订正以及将A和b应用于经过系统性偏差修正的雷达测量的反射率因子，估测雷达探测范围内的降水；所述的衰减订正方法建立的关系式k＝c*Zd应用于60公里范围内的降水；所述的X波段多普勒雷达的衰减订正方法包括以下步骤：

a)、X波段多普勒雷达垂直指向探测降水区，测量降水粒子的反射率因子和多普勒速度谱廓线；多普勒速度谱廓线还包括强度谱廓线和谱宽，以及根据速度谱廓线、强度谱廓线导出的二次产品，如回波顶高、回波底高、反射率因子、径向速度、谱宽、反射率因子剖面、径向速度剖面、谱宽剖面、组合反射率、垂直液态水含量、1小时累计含水量、3小时累计含水量、最大反射率、雨强因子，综合分析以上因子得出可能出现的降水量级；

b)、由降水粒子的多普勒速度谱廓线计算降水粒子的滴谱廓线；

c)、根据滴谱廓线计算米散射下的衰减系数廓线和降水粒子的反射率因子廓线；

d)、根据计算得到的衰减系数廓线和降水粒子的反射率因子廓线建立衰减系数k和反射率因子Z之间的关系，拟合确定系数c和d，从而确定衰减订正的关系式k＝c*Zd；

9)降水估测结束；

所述的基于X波段多普勒雷达降水估测方法的系统，包括X波段多普勒雷达、雨滴谱仪、监控主机和显示系统，所述X波段多普勒雷达包括外机和内机两个部分，所述的外机部分包括天线伺服分系统、发射分系统、接收分系统、监控标校分系统和配电系统；所述的内机部分包括信号处理器和用户分系统。

进一步，所述监控主机设置有同步正交跳频信号盲源分离模块，所述步正交跳频信号盲源分离的信号处理方法包括：

步骤一，利用含有M个阵元的阵列天线接收来自多个同步正交跳频电台的跳频信号，对每一路接收信号进行采样，得到采样后的M路离散时域混合信号m＝1,2,…,M；采集阵列天线节点间不同时间片的交互次数，根据得到的数据建立时间序列，通过三次指数平滑法来预测节点间下一个时间片的交互次数，将交互次数预测值与实际值的相对误差作为节点的直接信任值；直接信任值的具体计算步骤为：采集网络观测节点i与节点j之间的n个时间片的交互次数：选取一定时间间隔t作为一个观测时间片，以观测节点i和被测节点j在1个时间片内的交互次数作为观测指标，真实交互次数，记作y_t，依次记录n个时间片的y_n，并将其保存在节点i的通信记录表中；预测第n+1个时间片的交互次数：根据采集到的n个时间片的交互次数建立时间序列，采用三次指数平滑法预测下一个时间片n+1内节点i和j之间的交互次数，预测交互次数，记作计算公式如下：

${\hat{y}}_{n+1}=a_n+b_n+c_n$

预测系数a_n、b_n、c_n的取值可由如下公式计算得到：

$a_n=3{\hat{y}}_{n+1}^{\left(1\right)}-3{\hat{y}}_{n+1}^{\left(2\right)}+{\hat{y}}_{n+1}^{\left(3\right)}$

$b_n=\frac{\mathrm{\α}}{2{(1-\mathrm{\α})}^2}\[(6-5\mathrm{\α}){\hat{y}}_{n+1}^{\left(1\right)}-2(5-4\mathrm{\α}){\hat{y}}_{n+1}^{\left(2\right)}+(4-3\mathrm{\α}){\hat{y}}_{n+1}^{\left(3\right)}\]$

$c_n=\frac{{\mathrm{\α}}^2}{2{(1-\mathrm{\α})}^2}\[{\hat{y}}_{n+1}^{\left(1\right)}-2{\hat{y}}_{n+1}^{\left(2\right)}+{\hat{y}}_{n+1}^{\left(3\right)}\]$

其中：分别是一次、二次、三次指数平滑数，由如下公式计算得到：

${\hat{y}}_{n+1}^{\left(1\right)}=\mathrm{\α}\×y_n+(1-\mathrm{\α})\×{\hat{y}}_n^{\left(1\right)}$

${\hat{y}}_{n+1}^{\left(2\right)}=\mathrm{\α}\×{\hat{y}}_{n+1}^{\left(1\right)}+(1-\mathrm{\α})\×{\hat{y}}_n^{\left(2\right)}$

${\hat{y}}_{n+1}^{\left(3\right)}=\mathrm{\α}\×{\hat{y}}_{n+1}^{\left(2\right)}+(1-\mathrm{\α})\×{\hat{y}}_n^{\left(3\right)}$

是三次指数平滑法的初始值，其取值为

${\hat{y}}_0^{\left(1\right)}={\hat{y}}_0^{\left(2\right)}={\hat{y}}_0^{\left(3\right)}=\frac{y_1+y_2+y_3}{3}$

α是平滑系数(0＜α＜1)，体现信任的时间衰减特性，即离预测值越近的时间片的y_t权重越大，离预测值越远的时间片的y_t权重越小；一般地，如果数据波动较大，且长期趋势变化幅度较大，呈现明显迅速的上升或下降趋势时α应取较大值(0.6～0.8)，可以增加近期数据对预测结果的影响；当数据有波动，但长期趋势变化不大时，α可在0.1～0.4之间取值；如果数据波动平稳，α应取较小值(0.05～0.20)；

计算直接信任值：

节点j的直接信任值TD_ij为预测交互次数和真实交互次数y_n+1的相对误差，

采用多路径信任推荐方式而得到的计算式计算间接信任值；收集可信节点对节点j的直接信任值：节点i向所有满足TD_ik≤φ的可信关联节点询问其对节点j的直接信任值，其中φ为推荐节点的可信度阈值，根据可信度的要求精度，φ的取值范围为0～0.4；计算间接信任值：综合计算所收集到的信任值，得到节点j的间接信任值TR_ij，其中，Set(i)为观测节点i的关联节点中与j节点有过交互且其直接信任值满足TD_ik≤φ的节点集合；

由直接信任值和间接信任值整合计算得出综合信任值；综合信任值(T_ij)的计算公式如下：T_ij＝βTD_ij+(1-β)TR_ij，其中β(0≤β≤1)表示直接信任值的权重，当β＝0时，节点i和节点j没有直接交互关系，综合信任值的计算直接来自于间接信任值，判断较客观；当β＝1时，节点i对节点j的综合信任值全部来自于直接信任值，在这种情况下，判断较为主观，实际计算根据需要确定β的取值；

步骤二，对M路离散时域混合信号进行重叠加窗短时傅里叶变换，得到M个混合信号的时频域矩阵p＝0,1,…,P-1,q＝0,1,…,N_fft-1，其中P表示总的窗数，N_fft表示FFT变换长度；p，q)表示时频索引，具体的时频值为这里N_fft表示FFT变换的长度，p表示加窗次数，T_s表示采样间隔，f_s表示采样频率，C为整数，表示短时傅里叶变换加窗间隔的采样点数，C＜N_fft，且K_c＝N_fft/C为整数，也就是说采用的是重叠加窗的短时傅里叶变换；

步骤三，对步骤二中得到的跳频混合信号时频域矩阵进行预处理；

步骤四，利用聚类算法估计每一跳的跳变时刻以及各跳对应的归一化的混合矩阵列向量、跳频频率；在p(p＝0，1，2，…P-1)时刻，对表示的频率值进行聚类，得到的聚类中心个数表示p时刻存在的载频个数，个聚类中心则表示载频的大小，分别用表示；对每一采样时刻p(p＝0，1，2，…P-1)，利用聚类算法对进行聚类，同样可得到个聚类中心，用表示；对所有求均值并取整，得到源信号个数的估计即：

$\hat{N}=round\left(\frac{1}{p}\underset{p=0}{\overset{P-1}{\Σ}}{\hat{N}}_p\right);$

找出的时刻，用p_h表示，对每一段连续取值的p_h求中值，用表示第l段相连p_h的中值，则表示第l个频率跳变时刻的估计；根据估计得到的以及第四步中估计得到的频率跳变时刻估计出每一跳对应的个混合矩阵列向量具体公式为：

${\hat{a}}_n\left(l\right)=\left\{\begin{array}{cc}\frac{1}{{\stackrel{\‾}{p}}_h\left(1\right)}\·\underset{p=1,p\≠p_h}{\overset{{\stackrel{\‾}{p}}_h\left(1\right)}{\Σ}}{b}_{n,p}^0& l=1,\\ \frac{1}{{\stackrel{\‾}{p}}_h\left(l\right)-{\stackrel{\‾}{p}}_h(l-1)}\·\underset{p={\stackrel{\‾}{p}}_h(l-1)+1,p\≠p_h}{\overset{{\stackrel{\‾}{p}}_h\left(l\right)}{\Σ}}{b}_{n,p}^0& l>1,\end{array}\right.,n=1,2,...,\hat{N};$

这里表示第l跳对应的个混合矩阵列向量估计值；估计每一跳对应的载频频率，用表示第l跳对应的个频率估计值，计算公式如下：

${\hat{f}}_{c,n}\left(l\right)=\left\{\begin{array}{cc}\frac{1}{{\stackrel{\‾}{p}}_h\left(1\right)}\·\underset{p=1,p\≠p_h}{\overset{{\stackrel{\‾}{p}}_h\left(1\right)}{\Σ}}{f}_o^n\left(p\right)& l=1,\\ \frac{1}{{\stackrel{\‾}{p}}_h\left(l\right)-{\stackrel{\‾}{p}}_h(l-1)}\·\underset{p={\stackrel{\‾}{p}}_h(l-1)+1,p\≠p_h}{\overset{{\stackrel{\‾}{p}}_h\left(l\right)}{\Σ}}{f}_o^n\left(p\right)& l>1,\end{array}\right.,n=1,2,...,\hat{N};$

步骤五，根据步骤四估计得到的归一化混合矩阵列向量估计时频域跳频源信号；

步骤六，对不同跳频点之间的时频域跳频源信号进行拼接；

步骤七，根据源信号时频域估计值，恢复时域跳频源信号；对每一采样时刻p(p＝0，1，2，…)的频域数据Y_n(p,q),q＝0,1,2,…,N_fft-1做N_fft点的IFFT变换，得到p采样时刻对应的时域跳频源信号，用y_n(p,q_t)(q_t＝0,1,2,…,N_fft-1)表示；对上述所有时刻得到的时域跳频源信号y_n(p,q_t)进行合并处理，得到最终的时域跳频源信号估计，具体公式如下：

$s_n\&lsqb;kC:(k+1)C-1\&rsqb;=\left\{\begin{array}{cc}\underset{m=0}{\overset{k}{\mathrm{\&Sigma;}}}{y}_n\&lsqb;m,(k-m)C:(k-m+1)C-1\&rsqb;& k<K_c\\ \underset{m=k-K_c+1}{\overset{k}{\mathrm{\&Sigma;}}}{y}_n\&lsqb;m,(k-m)C:(k-m+1)C-1\&rsqb;& k\&GreaterEqual;K_c\end{array}\right.,k=0,1,2,...$

这里K_c＝N_fft/C，C为短时傅里叶变换加窗间隔的采样点数，N_fft为FFT变换的长度。

进一步，所述信号处理器设置有求解单元，所述求解单元的信号处理方法包括：

第一步，保证次级用户网络能长时间工作，需对次级用户的发射功率进行限制，保证次级用户网络的平均发射功率低于限定值：

E{α₀P₀+α₁P₁+β₀P₀+β₁P₁}≤P_av；

式中P_av是次级用户发射机SU-Tx的最大平均发射功率，这的平均是指信道衰减系数h_i,g_ss,g_sp随机变量的期望；

第二步，认知无线电网络的首要任务是保护主用户网络的服务质量，故对网络的干扰功率进行了限制；根据基于合作感知的频谱共享网络模型，知道干扰只在主用户PU处于忙状态时发生，所以平均干扰功率约束写成如下形式：

E{g_sp(β₀P₀+β₁P₁)}≤Q_av；

第三步，确保各个节点处的检测概率和网络的整体检测概率分别不低于各自的目标检测概率，关于检测概率的限制条件如下：

P_d≥P_th,P_di≥p_th,i＝1,2…k；

第四步，根据上述限制条件下，建立以最大化次级网络的平均吞吐量为目标函数的最优化问题：

$\begin{array}{cc}\underset{\{{\mathrm{\&tau;}}_s,\mathrm{\&epsiv;},\{{\mathrm{\&epsiv;}}_i\},P_0,P_1\}}{\max }& R\\ subjectto& \left(2\right),\left(3\right),\left(4\right),P_0\&GreaterEqual;0,P_1\&GreaterEqual;0\\ & 0\&le;{\mathrm{\&tau;}}_s\&le;T\end{array}---\left(\mathrm{Pr}oblem1\right)$

第五步，求解所建立的最优化问题，选择使得次级网络的吞吐量最大的合作感知的感知周期和次级用户的信号发射功率作为该频谱感知模型的感知参数；

具体包括以下的步骤：

1)，对不等式约束条件组(4)取等号，简化Problem 1为Problem 2；

$\begin{array}{cc}\underset{\{{\mathrm{\&tau;}}_s,P_0,P_1\}}{\max }& R\\ subjectto& \left(2\right),\left(3\right),P_0\&GreaterEqual;0,P_1\&GreaterEqual;0\\ & 0\&le;{\mathrm{\&tau;}}_s\&le;T\end{array}---\left(\mathrm{Pr}oblem2\right)$

2)，弱化对感知周期τ_s的求解，重点求解使平均吞吐量最大化的信号发射功率P₀，P₁；关于发射功率P₀和P₁的拉格朗日函数如下：

$\begin{array}{c}L(P_0,P_1,\mathrm{\&lambda;},\mathrm{\&mu;})=E\{\frac{T-{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_s}{T}\&lsqb;{\mathrm{\&alpha;}}_0{r}_{00}+{\mathrm{\&alpha;}}_1{r}_{01}+{\mathrm{\&beta;}}_0{r}_{10}+{\mathrm{\&beta;}}_1{r}_{11}\&rsqb;\}\\ -\mathrm{\&lambda;}\&lsqb;E\{{\mathrm{\&alpha;}}_0{P}_0+{\mathrm{\&alpha;}}_1{P}_1+{\mathrm{\&beta;}}_0{P}_0+{\mathrm{\&beta;}}_1{P}_1\}-P_{av}\&rsqb;\\ -\mathrm{\&mu;}\&lsqb;E\left\{g_{sp}({\mathrm{\&beta;}}_0{P}_0+{\mathrm{\&beta;}}_1{P}_1)\right\}-Q_{av}\&rsqb;\end{array}.$

所以Problem 2的拉格朗日对偶优化问题为：

$\begin{array}{cc}\underset{\mathrm{\&lambda;}\&GreaterEqual;0,\mathrm{\&mu;}\&GreaterEqual;0}{\min imize}& g(\mathrm{\&lambda;},\mathrm{\&mu;})\end{array}---\left(\mathrm{Pr}oblem3\right)$

其中表示拉格朗日对偶函数；证明优化问题Problem 2与Problem 3的最优值差值为零，说明优化问题Problem 2与其拉格朗日对偶优化问题Problem 3之间是等价的，故只需求Problem 3的最优解即可；该问题是一个关于双变量P₀P₁的联合规划问题，为此将分解成两个子优化问题：

SP1：

SP2：

看出SP1和SP2分别是关于P₀P₁的无约束凸优化问题，此时运用拉格朗日函数以及KKT条件，便得到当检测到主用户PU处于闲状态时次级用户发射机SU-Tx的最优发射功率：

$P_0={\&lsqb;\frac{A_0+\sqrt{{\mathrm{\&Lambda;}}_0}}{2}\&rsqb;}^{+};$

其中：

${\mathrm{\&Lambda;}}_0=A_0^2-\frac{4}{g_{ss}}\{\frac{{\mathrm{\&sigma;}}_u^4+{\mathrm{\&sigma;}}_u^2{h}_k{P}_p}{g_{ss}}-\frac{{\log }_2\left(e\right)\&lsqb;{\mathrm{\&alpha;}}_0({\mathrm{\&sigma;}}_u^2+h_k{P}_p)+{\mathrm{\&beta;}}_0{\mathrm{\&sigma;}}_u^2\&rsqb;}{\mathrm{\&lambda;}({\mathrm{\&alpha;}}_0+{\mathrm{\&beta;}}_0)+{\mathrm{\&mu;\&beta;}}_0{g}_{sp}}\}$

当检测到主用户PU处于忙状态时，次级用户发射机SU-Tx的最优发射功率为：

$P_1={\&lsqb;\frac{A_1+\sqrt{{\mathrm{\&Lambda;}}_1}}{2}\&rsqb;}^{+};$

其中：

${\mathrm{\&Lambda;}}_1=A_1^2-\frac{4}{g_{ss}}\{\frac{{\mathrm{\&sigma;}}_u^4+{\mathrm{\&sigma;}}_u^2{h}_k{P}_p}{g_{ss}}-\frac{{\log }_2\left(e\right)\&lsqb;{\mathrm{\&alpha;}}_1({\mathrm{\&sigma;}}_u^2+h_k{P}_p)+{\mathrm{\&beta;}}_1{\mathrm{\&sigma;}}_u^2\&rsqb;}{\mathrm{\&lambda;}({\mathrm{\&alpha;}}_1+{\mathrm{\&beta;}}_1)+{\mathrm{\&mu;\&beta;}}_1{g}_{sp}}\}$

式中[x]⁺＝max{0,x}；λ≥0，μ≥0是式E{α₀P₀+α₁P₁+β₀P₀+β₁P₁}≤P_av和E{g_sp(β₀P₀+β₁P₁)}≤Q_av的拉格朗日乘子。

本发明通过X波段多普勒雷达垂直指向探测降水区，测量得到多普勒速度谱廓线，导出降水粒子的滴谱廓线，从而能够计算降水粒子的衰减系数，结合计算的反射率因子廓线，就可以进一步建立衰减系数k和反射率因子之间的关系k＝c*Zd，并将这一关系应用于60公里范围内的降水，解决X波段多普勒雷达的衰减订正问题，扩大了X波段多普勒雷达在降水方面的应用范围。进一步，从滴谱廓线或雨滴谱仪测量的雨滴谱计算得到雨强，拟合Z＝A*Rb确定系数A和b，从而实现雷达估测降水。整个过程从多普勒雷达实际测量的速度数据出发，不需要太多的理论假定，只是在正常观测过程中，增加一次垂直指向探测就可以实现。本发明在不知道任何信息的条件下，仅根据接收到的多个跳频信号的混合信号，估计出跳频源信号，能在接收天线个数小于源信号个数的条件下，对多个跳频信号进行盲估计，仅仅利用了短时傅里叶变换，计算量小，容易实现，该方法在对跳频信号进行盲分离的同时，还能对部分参数进行估计，实用性强，具有较强的推广与应用价值。

附图说明

图1是本发明实施例提供的基于X波段多普勒雷达降水估测系统的结构示意图；

图2是本发明实施例提供的基于X波段多普勒雷达降水估测方法的流程示意图；

图3是本发明实施例提供的基于X波段多普勒雷达衰减订正的方法示意图。

图中：1、X波段多普勒雷达；1-1、外机；1-1-1、天线伺服分系统；1-1-2、发射分系统；1-1-3、接收分系统；1-1-4、监控标校分系统；1-1-5、配电系统；1-2、内机；1-2-1、信号处理器；1-2-2、用户分系统；2、雨滴谱仪。

具体实施方式

为能进一步了解本发明的发明内容、特点及功效，兹例举以下实施例，并配合附图详细说明如下。

请参阅附图1至图3：

本发明提供一种基于X波段多普勒雷达降水估测系统及方法，该基于X波段多普勒雷达降水估测系统包括X波段多普勒雷达1、雨滴谱仪2、监控主机和显示系统，所述X波段多普勒雷达1包括外机1-1和内机1-2两个部分，所述的外机1-1部分包括天线伺服分系统1-1-1、发射分系统1-1-2、接收分系统1-1-3、监控标校分系统1-1-4和配电系统1-1-5；所述的内机1-2部分包括信号处理器1-2-1和用户分系统1-2-2；

一种基于X波段多普勒雷达的降水估测方法，利用一部X波段多普勒雷达和周边的雨滴谱仪联合估测雷达有效探测范围内的降水场，包括以下步骤：

S1、用无地物的最低仰角作平面位置扫描，得到雷达测量的反射率因子场；

S2、X波段多普勒雷达1作垂直指向探测降水区，判断本站即雷达处是否有降水，

(1)当本站有降水时，测量降水粒子的反射率因子和多普勒速度谱廓线，并由降水粒子的多普勒速度谱廓线计算降水粒子的滴谱廓线，转步骤3)；

(2)当本站没有降水而某个雨滴谱仪2处有降水时，雨滴谱仪2得到雨滴谱，转步骤3)；

(3)当本站和雨滴谱仪2所在位置都没有降水时，用缺省的系数c和d，进行衰减订正，并用缺省的系数A和b进行降水估测，转步骤6)；

S3、计算降水区的衰减系数k、降水粒子的反射率因子Z和雨强R；

S4、根据衰减系数k和反射率因子Z，拟合k＝c*Zd确定系数c和d，将此关系应用于降水区域，实现雷达测量的反射率因子的衰减订正，并根据反射率因子Z和雨强R，拟合Z＝A*Rb确定系数A和b；

S5、根据衰减订正以及将A和b应用于经过系统性偏差修正的雷达测量的反射率因子，估测雷达探测范围内的降水；

S6、降水估测结束。

进一步，在步骤S1中，所述最低仰角≤18°。

进一步，在步骤S4中，通过最小二乘法拟合确定系数c和d。

进一步，在步骤S5中，所述系统性偏差修正为在雨滴谱仪2所在格点周边，计算雨滴谱得到的反射率因子与雷达测量得到的反射率因子的差值，视最小的差值为最佳匹配，这样，每个雨滴谱仪2都有一个最小差值，计算这些最小差值的平均值，视此平均值为雷达测量的反射率因子的系统性偏差，并对雷达测量的反射率因子进行订正。

进一步，所述的天线伺服分系统1-1-1采用旋转抛物面天线，并设有旋转抛物面保护罩；

该分系统的参数为：

1)、频率：X波段；

2)、极化方式：水平、垂直的双线性极化；

3)、直径：1.8m；

4)、增益：≥40dB；

5)、副瓣电平：≤-27dB；

6)、电轴指向波束主轴方向差：≤0.1。；

7)、天线扫描方式：PPI、RHI、体扫、扇扫、定点；

8)、天线扫描范围：方位0～360。；俯仰0～90。；

9)、天线扫描速度：方位1～4rpm；俯仰1～2rpm。

进一步，所述的X波段多普勒雷达1的衰减订正方法包括以下步骤：

S21、X波段多普勒雷达垂直指向探测降水区，测量降水粒子的反射率因子和多普勒速度谱廓线；多普勒速度谱廓线还包括强度谱廓线和谱宽，以及根据速度谱廓线、强度谱廓线导出的二次产品，如回波顶高、回波底高、反射率因子、径向速度、谱宽、反射率因子剖面、径向速度剖面、谱宽剖面、组合反射率、垂直液态水含量、1小时累计含水量、3小时累计含水量、最大反射率、雨强因子，综合分析以上因子得出可能出现的降水量级；

S22、由降水粒子的多普勒速度谱廓线计算降水粒子的滴谱廓线；

S23、根据滴谱廓线计算米散射下的衰减系数廓线和降水粒子的反射率因子廓线；

S24、根据计算得到的衰减系数廓线和降水粒子的反射率因子廓线建立衰减系数k和反射率因子Z之间的关系，拟合确定系数c和d，从而确定衰减订正的关系式k＝c*Zd。

进一步，所述的衰减订正方法建立的关系式k＝c*Zd应用于60公里范围内的降水。

进一步，所述监控主机设置有同步正交跳频信号盲源分离模块，所述步正交跳频信号盲源分离的信号处理方法包括：

步骤一，利用含有M个阵元的阵列天线接收来自多个同步正交跳频电台的跳频信号，对每一路接收信号进行采样，得到采样后的M路离散时域混合信号m＝1,2,…,M；采集阵列天线节点间不同时间片的交互次数，根据得到的数据建立时间序列，通过三次指数平滑法来预测节点间下一个时间片的交互次数，将交互次数预测值与实际值的相对误差作为节点的直接信任值；直接信任值的具体计算步骤为：采集网络观测节点i与节点j之间的n个时间片的交互次数：选取一定时间间隔t作为一个观测时间片，以观测节点i和被测节点j在1个时间片内的交互次数作为观测指标，真实交互次数，记作y_t，依次记录n个时间片的y_n，并将其保存在节点i的通信记录表中；预测第n+1个时间片的交互次数：根据采集到的n个时间片的交互次数建立时间序列，采用三次指数平滑法预测下一个时间片n+1内节点i和j之间的交互次数，预测交互次数，记作计算公式如下：

${\hat{y}}_{n+1}=a_n+b_n+c_n$

预测系数a_n、b_n、c_n的取值可由如下公式计算得到：

$a_n=3{\hat{y}}_{n+1}^{\left(1\right)}-3{\hat{y}}_{n+1}^{\left(2\right)}+{\hat{y}}_{n+1}^{\left(3\right)}$

$b_n=\frac{\mathrm{\&alpha;}}{2{(1-\mathrm{\&alpha;})}^2}\&lsqb;(6-5\mathrm{\&alpha;}){\hat{y}}_{n+1}^{\left(1\right)}-2(5-4\mathrm{\&alpha;}){\hat{y}}_{n+1}^{\left(2\right)}+(4-3\mathrm{\&alpha;}){\hat{y}}_{n+1}^{\left(3\right)}\&rsqb;$

$c_n=\frac{{\mathrm{\&alpha;}}^2}{2{(1-\mathrm{\&alpha;})}^2}\&lsqb;{\hat{y}}_{n+1}^{\left(1\right)}-2{\hat{y}}_{n+1}^{\left(2\right)}+{\hat{y}}_{n+1}^{\left(3\right)}\&rsqb;$

其中：分别是一次、二次、三次指数平滑数，由如下公式计算得到：

${\hat{y}}_{n+1}^{\left(1\right)}=\mathrm{\&alpha;}\&times;y_n+(1-\mathrm{\&alpha;})\&times;{\hat{y}}_n^{\left(1\right)}$

${\hat{y}}_{n+1}^{\left(2\right)}=\mathrm{\&alpha;}\&times;{\hat{y}}_{n+1}^{\left(1\right)}+(1-\mathrm{\&alpha;})\&times;{\hat{y}}_n^{\left(2\right)}$

${\hat{y}}_{n+1}^{\left(3\right)}=\mathrm{\&alpha;}\&times;{\hat{y}}_{n+1}^{\left(2\right)}+(1-\mathrm{\&alpha;})\&times;{\hat{y}}_n^{\left(3\right)}$

是三次指数平滑法的初始值，其取值为

${\hat{y}}_0^{\left(1\right)}={\hat{y}}_0^{\left(2\right)}={\hat{y}}_0^{\left(3\right)}=\frac{y_1+y_2+y_3}{3}$

α是平滑系数(0＜α＜1)，体现信任的时间衰减特性，即离预测值越近的时间片的y_t权重越大，离预测值越远的时间片的y_t权重越小；一般地，如果数据波动较大，且长期趋势变化幅度较大，呈现明显迅速的上升或下降趋势时α应取较大值(0.6～0.8)，可以增加近期数据对预测结果的影响；当数据有波动，但长期趋势变化不大时，α可在0.1～0.4之间取值；如果数据波动平稳，α应取较小值(0.05～0.20)；

计算直接信任值：

节点j的直接信任值TD_ij为预测交互次数和真实交互次数y_n+1的相对误差，

采用多路径信任推荐方式而得到的计算式计算间接信任值；收集可信节点对节点j的直接信任值：节点i向所有满足TD_ik≤φ的可信关联节点询问其对节点j的直接信任值，其中φ为推荐节点的可信度阈值，根据可信度的要求精度，φ的取值范围为0～0.4；计算间接信任值：综合计算所收集到的信任值，得到节点j的间接信任值TR_ij，其中，Set(i)为观测节点i的关联节点中与j节点有过交互且其直接信任值满足TD_ik≤φ的节点集合；

由直接信任值和间接信任值整合计算得出综合信任值；综合信任值(T_ij)的计算公式如下：T_ij＝βTD_ij+(1-β)TR_ij，其中β(0≤β≤1)表示直接信任值的权重，当β＝0时，节点i和节点j没有直接交互关系，综合信任值的计算直接来自于间接信任值，判断较客观；当β＝1时，节点i对节点j的综合信任值全部来自于直接信任值，在这种情况下，判断较为主观，实际计算根据需要确定β的取值；

步骤二，对M路离散时域混合信号进行重叠加窗短时傅里叶变换，得到M个混合信号的时频域矩阵p＝0,1,…,P-1,q＝0,1,…,N_fft-1，其中P表示总的窗数，N_fft表示FFT变换长度；p，q)表示时频索引，具体的时频值为这里N_fft表示FFT变换的长度，p表示加窗次数，T_s表示采样间隔，f_s表示采样频率，C为整数，表示短时傅里叶变换加窗间隔的采样点数，C＜N_fft，且K_c＝N_fft/C为整数，也就是说采用的是重叠加窗的短时傅里叶变换；

步骤三，对步骤二中得到的跳频混合信号时频域矩阵进行预处理；

步骤四，利用聚类算法估计每一跳的跳变时刻以及各跳对应的归一化的混合矩阵列向量、跳频频率；在p(p＝0，1，2，…P-1)时刻，对表示的频率值进行聚类，得到的聚类中心个数表示p时刻存在的载频个数，个聚类中心则表示载频的大小，分别用表示；对每一采样时刻p(p＝0，1，2，…P-1)，利用聚类算法对进行聚类，同样可得到个聚类中心，用表示；对所有求均值并取整，得到源信号个数的估计即：

$\hat{N}=round\left(\frac{1}{p}\underset{p=0}{\overset{P-1}{\&Sigma;}}{\hat{N}}_p\right);$

找出的时刻，用p_h表示，对每一段连续取值的p_h求中值，用表示第l段相连p_h的中值，则表示第l个频率跳变时刻的估计；根据估计得到的以及第四步中估计得到的频率跳变时刻估计出每一跳对应的个混合矩阵列向量具体公式为：

${\hat{a}}_n\left(l\right)=\left\{\begin{array}{cc}\frac{1}{{\stackrel{\&OverBar;}{p}}_h\left(1\right)}\&CenterDot;\underset{p=1,p\&NotEqual;p_h}{\overset{{\stackrel{\&OverBar;}{p}}_h\left(1\right)}{\&Sigma;}}{b}_{n,p}^0& l=1,\\ \frac{1}{{\stackrel{\&OverBar;}{p}}_h\left(l\right)-{\stackrel{\&OverBar;}{p}}_h(l-1)}\&CenterDot;\underset{p={\stackrel{\&OverBar;}{p}}_h(l-1)+1,p\&NotEqual;p_h}{\overset{{\stackrel{\&OverBar;}{p}}_h\left(l\right)}{\&Sigma;}}{b}_{n,p}^0& l>1,\end{array}\right.,n=1,2,...,\hat{N};$

这里表示第l跳对应的个混合矩阵列向量估计值；估计每一跳对应的载频频率，用表示第l跳对应的个频率估计值，计算公式如下：

${\hat{f}}_{c,n}\left(l\right)=\left\{\begin{array}{cc}\frac{1}{{\stackrel{\&OverBar;}{p}}_h\left(1\right)}\&CenterDot;\underset{p=1,p\&NotEqual;p_h}{\overset{{\stackrel{\&OverBar;}{p}}_h\left(1\right)}{\&Sigma;}}{f}_o^n\left(p\right)& l=1,\\ \frac{1}{{\stackrel{\&OverBar;}{p}}_h\left(l\right)-{\stackrel{\&OverBar;}{p}}_h(l-1)}\&CenterDot;\underset{p={\stackrel{\&OverBar;}{p}}_h(l-1)+1,p\&NotEqual;p_h}{\overset{{\stackrel{\&OverBar;}{p}}_h\left(l\right)}{\&Sigma;}}{f}_o^n\left(p\right)& l>1,\end{array}\right.,n=1,2,...,\hat{N};$

步骤五，根据步骤四估计得到的归一化混合矩阵列向量估计时频域跳频源信号；

步骤六，对不同跳频点之间的时频域跳频源信号进行拼接；

步骤七，根据源信号时频域估计值，恢复时域跳频源信号；对每一采样时刻p(p＝0，1，2，…)的频域数据Y_n(p,q),q＝0,1,2,…,N_fft-1做N_fft点的IFFT变换，得到p采样时刻对应的时域跳频源信号，用y_n(p,q_t)(q_t＝0,1,2,…,N_fft-1)表示；对上述所有时刻得到的时域跳频源信号y_n(p,q_t)进行合并处理，得到最终的时域跳频源信号估计，具体公式如下：

$s_n\&lsqb;kC:(k+1)C-1\&rsqb;=\left\{\begin{array}{cc}\underset{m=0}{\overset{k}{\mathrm{\&Sigma;}}}{y}_n\&lsqb;m,(k-m)C:(k-m+1)C-1\&rsqb;& k<K_c\\ \underset{m=k-K_c+1}{\overset{k}{\mathrm{\&Sigma;}}}{y}_n\&lsqb;m,(k-m)C:(k-m+1)C-1\&rsqb;& k\&GreaterEqual;K_c\end{array}\right.,k=0,1,2,...$

这里K_c＝N_fft/C，C为短时傅里叶变换加窗间隔的采样点数，N_fft为FFT变换的长度。

进一步，所述信号处理器设置有求解单元，所述求解单元的信号处理方法包括：

第一步，保证次级用户网络能长时间工作，需对次级用户的发射功率进行限制，保证次级用户网络的平均发射功率低于限定值：

E{α₀P₀+α₁P₁+β₀P₀+β₁P₁}≤P_av；

式中P_av是次级用户发射机SU-Tx的最大平均发射功率，这的平均是指信道衰减系数h_i,g_ss,g_sp随机变量的期望；

第二步，认知无线电网络的首要任务是保护主用户网络的服务质量，故对网络的干扰功率进行了限制；根据基于合作感知的频谱共享网络模型，知道干扰只在主用户PU处于忙状态时发生，所以平均干扰功率约束写成如下形式：

E{g_sp(β₀P₀+β₁P₁)}≤Q_av；

第三步，确保各个节点处的检测概率和网络的整体检测概率分别不低于各自的目标检测概率，关于检测概率的限制条件如下：

P_d≥P_th,P_di≥p_th,i＝1,2…k；

第四步，根据上述限制条件下，建立以最大化次级网络的平均吞吐量为目标函数的最优化问题：

$\begin{array}{cc}\underset{\{{\mathrm{\&tau;}}_s,\mathrm{\&epsiv;},\{{\mathrm{\&epsiv;}}_i\},P_0,P_1\}}{\max }& R\\ subjectto& \left(2\right),\left(3\right),\left(4\right),P_0\&GreaterEqual;0,P_1\&GreaterEqual;0\\ & 0\&le;{\mathrm{\&tau;}}_s\&le;T\end{array}---\left(\mathrm{Pr}oblem1\right)$

第五步，求解所建立的最优化问题，选择使得次级网络的吞吐量最大的合作感知的感知周期和次级用户的信号发射功率作为该频谱感知模型的感知参数；

具体包括以下的步骤：

1)，对不等式约束条件组(4)取等号，简化Problem 1为Problem 2；

$\begin{array}{cc}\underset{\{{\mathrm{\&tau;}}_s,P_0,P_1\}}{\max }& R\\ subjectto& \left(2\right),\left(3\right),P_0\&GreaterEqual;0,P_1\&GreaterEqual;0\\ & 0\&le;{\mathrm{\&tau;}}_s\&le;T\end{array}---\left(\mathrm{Pr}oblem2\right)$

2)，弱化对感知周期τ_s的求解，重点求解使平均吞吐量最大化的信号发射功率P₀，P₁；关于发射功率P₀和P₁的拉格朗日函数如下：

$\begin{array}{c}L(P_0,P_1,\mathrm{\&lambda;},\mathrm{\&mu;})=E\{\frac{T-{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_s}{T}\&lsqb;{\mathrm{\&alpha;}}_0{r}_{00}+{\mathrm{\&alpha;}}_1{r}_{01}+{\mathrm{\&beta;}}_0{r}_{10}+{\mathrm{\&beta;}}_1{r}_{11}\&rsqb;\}\\ -\mathrm{\&lambda;}\&lsqb;E\{{\mathrm{\&alpha;}}_0{P}_0+{\mathrm{\&alpha;}}_1{P}_1+{\mathrm{\&beta;}}_0{P}_0+{\mathrm{\&beta;}}_1{P}_1\}-P_{av}\&rsqb;\\ -\mathrm{\&mu;}\&lsqb;E\left\{g_{sp}({\mathrm{\&beta;}}_0{P}_0+{\mathrm{\&beta;}}_1{P}_1)\right\}-Q_{av}\&rsqb;\end{array}.$

所以Problem 2的拉格朗日对偶优化问题为：

$\begin{array}{cc}\underset{\mathrm{\&lambda;}\&GreaterEqual;0,\mathrm{\&mu;}\&GreaterEqual;0}{\min imize}& g(\mathrm{\&lambda;},\mathrm{\&mu;})\end{array}---\left(\mathrm{Pr}oblem3\right)$

其中表示拉格朗日对偶函数；证明优化问题Problem 2与Problem 3的最优值差值为零，说明优化问题Problem 2与其拉格朗日对偶优化问题Problem 3之间是等价的，故只需求Problem 3的最优解即可；该问题是一个关于双变量P₀P₁的联合规划问题，为此将分解成两个子优化问题：

SP1：

SP2：

看出SP1和SP2分别是关于P₀P₁的无约束凸优化问题，此时运用拉格朗日函数以及KKT条件，便得到当检测到主用户PU处于闲状态时次级用户发射机SU-Tx的最优发射功率：

$P_0={\&lsqb;\frac{A_0+\sqrt{{\mathrm{\&Lambda;}}_0}}{2}\&rsqb;}^{+};$

其中：

${\mathrm{\&Lambda;}}_0=A_0^2-\frac{4}{g_{ss}}\{\frac{{\mathrm{\&sigma;}}_u^4+{\mathrm{\&sigma;}}_u^2{h}_k{P}_p}{g_{ss}}-\frac{{\log }_2\left(e\right)\&lsqb;{\mathrm{\&alpha;}}_0({\mathrm{\&sigma;}}_u^2+h_k{P}_p)+{\mathrm{\&beta;}}_0{\mathrm{\&sigma;}}_u^2\&rsqb;}{\mathrm{\&lambda;}({\mathrm{\&alpha;}}_0+{\mathrm{\&beta;}}_0)+{\mathrm{\&mu;\&beta;}}_0{g}_{sp}}\}$

当检测到主用户PU处于忙状态时，次级用户发射机SU-Tx的最优发射功率为：

$P_1={\&lsqb;\frac{A_1+\sqrt{{\mathrm{\&Lambda;}}_1}}{2}\&rsqb;}^{+};$

其中：

${\mathrm{\&Lambda;}}_1=A_1^2-\frac{4}{g_{ss}}\{\frac{{\mathrm{\&sigma;}}_u^4+{\mathrm{\&sigma;}}_u^2{h}_k{P}_p}{g_{ss}}-\frac{{\log }_2\left(e\right)\&lsqb;{\mathrm{\&alpha;}}_1({\mathrm{\&sigma;}}_u^2+h_k{P}_p)+{\mathrm{\&beta;}}_1{\mathrm{\&sigma;}}_u^2\&rsqb;}{\mathrm{\&lambda;}({\mathrm{\&alpha;}}_1+{\mathrm{\&beta;}}_1)+{\mathrm{\&mu;\&beta;}}_1{g}_{sp}}\}$

式中[x]⁺＝max{0,x}；λ≥0，μ≥0是式E{α₀P₀+α₁P₁+β₀P₀+β₁P₁}≤P_av和E{g_sp(β₀P₀+β₁P₁)}≤Q_av的拉格朗日乘子。

本发明通过X波段多普勒雷达垂直指向探测降水区，测量得到多普勒速度谱廓线，导出降水粒子的滴谱廓线，从而能够计算降水粒子的衰减系数，结合计算的反射率因子廓线，就可以进一步建立衰减系数k和反射率因子之间的关系k＝c*Zd，并将这一关系应用于60公里范围内的降水，解决X波段多普勒雷达的衰减订正问题，扩大了X波段多普勒雷达在降水方面的应用范围。进一步，从滴谱廓线或雨滴谱仪测量的雨滴谱计算得到雨强，拟合Z＝A*Rb确定系数A和b，从而实现雷达估测降水。整个过程从多普勒雷达实际测量的速度数据出发，不需要太多的理论假定，只是在正常观测过程中，增加一次垂直指向探测就可以实现。

利用本发明所述的技术方案，或本领域的技术人员在本发明技术方案的启发下，设计出类似的技术方案，而达到上述技术效果的，均是落入本发明的保护范围。

Claims (3)

1.一种基于X波段多普勒雷达的降水估测方法，其特征在于，所述基于X波段多普勒雷达的降水估测方法包括以下步骤：

1)、用无地物的最低仰角作平面位置扫描，得到雷达测量的反射率因子场；最低仰角≤18°；

2)、X波段多普勒雷达作垂直指向探测降水区，判断探测降水区即雷达处是否有降水；

3)当探测降水区有降水时，测量降水粒子的反射率因子和多普勒速度谱廓线，并由降水粒子的多普勒速度谱廓线计算降水粒子的滴谱廓线，转步骤4)；

4)当探测降水区没有降水而某个雨滴谱仪处有降水时，雨滴谱仪得到雨滴谱，转步骤3)；

5)当探测降水区和雨滴谱仪所在位置都没有降水时，用缺省的系数c和d，进行衰减订正，并用缺省的系数A和b进行降水估测，转步骤6)；通过最小二乘法拟合确定系数c和d；所述系统性偏差修正为在雨滴谱仪所在格点周边，计算雨滴谱得到的反射率因子与雷达测量得到的反射率因子的差值，视最小的差值为最佳匹配，这样，每个雨滴谱仪都有一个最小差值，计算这些最小差值的平均值，视此平均值为雷达测量的反射率因子的系统性偏差，并对雷达测量的反射率因子进行订正；

6)计算降水区的衰减系数k、降水粒子的反射率因子Z和雨强R；

7)根据衰减系数k和反射率因子Z，拟合k＝c*Zd确定系数c和d，将此关系应用于降水区域，实现雷达测量的反射率因子的衰减订正，并根据反射率因子Z和雨强R，拟合Z＝A*Rb确定系数A和b；

8)根据衰减订正以及将A和b应用于经过系统性偏差修正的雷达测量的反射率因子，估测雷达探测范围内的降水；所述的衰减订正方法建立的关系式k＝c*Zd应用于60公里范围内的降水；所述的X波段多普勒雷达的衰减订正方法包括以下步骤：

a)、X波段多普勒雷达垂直指向探测降水区，测量降水粒子的反射率因子和多普勒速度谱廓线；多普勒速度谱廓线还包括强度谱廓线和谱宽，以及根据速度谱廓线、强度谱廓线导出的二次产品，如回波顶高、回波底高、反射率因子、径向速度、谱宽、反射率因子剖面、径向速度剖面、谱宽剖面、组合反射率、垂直液态水含量、1小时累计含水量、3小时累计含水量、最大反射率、雨强因子，综合分析以上因子得出可能出现的降水量级；

b)、由降水粒子的多普勒速度谱廓线计算降水粒子的滴谱廓线；

c)、根据滴谱廓线计算米散射下的衰减系数廓线和降水粒子的反射率因子廓线；

d)、根据计算得到的衰减系数廓线和降水粒子的反射率因子廓线建立衰减系数k和反射率因子Z之间的关系，拟合确定系数c和d，从而确定衰减订正的关系式k＝c*Zd；

9)降水估测结束；

所述的基于X波段多普勒雷达降水估测方法的系统，包括X波段多普勒雷达、雨滴谱仪、监控主机和显示系统，所述X波段多普勒雷达包括外机和内机两个部分，所述的外机部分包括天线伺服分系统、发射分系统、接收分系统、监控标校分系统和配电系统；所述的内机部分包括信号处理器和用户分系统。

2.如权利要求1所述的基于X波段多普勒雷达的降水估测方法，其特征在于，所述监控主机设置有同步正交跳频信号盲源分离模块，所述步正交跳频信号盲源分离的信号处理方法包括：

步骤一，利用含有M个阵元的阵列天线接收来自多个同步正交跳频电台的跳频信号，对每一路接收信号进行采样，得到采样后的M路离散时域混合信号采集阵列天线节点间不同时间片的交互次数，根据得到的数据建立时间序列，通过三次指数平滑法来预测节点间下一个时间片的交互次数，将交互次数预测值与实际值的相对误差作为节点的直接信任值；直接信任值的具体计算步骤为：采集网络观测节点i与节点j之间的n个时间片的交互次数：选取一定时间间隔t作为一个观测时间片，以观测节点i和被测节点j在1个时间片内的交互次数作为观测指标，真实交互次数，记作y_t，依次记录n个时间片的y_n，并将其保存在节点i的通信记录表中；预测第n+1个时间片的交互次数：根据采集到的n个时间片的交互次数建立时间序列，采用三次指数平滑法预测下一个时间片n+1内节点i和j之间的交互次数，预测交互次数，记作计算公式如下：

${\hat{y}}_{n+1}=a_n+b_n+c_n$

预测系数a_n、b_n、c_n的取值可由如下公式计算得到：

$a_n=3{\hat{y}}_{n+1}^{\left(1\right)}-3{\hat{y}}_{n+1}^{\left(2\right)}+{\hat{y}}_{n+1}^{\left(3\right)}$

$b_n=\frac{\mathrm{\&alpha;}}{2{(1-\mathrm{\&alpha;})}^2}\&lsqb;(6-5\mathrm{\&alpha;}){\hat{y}}_{n+1}^{\left(1\right)}-2(5-4\mathrm{\&alpha;}){\hat{y}}_{n+1}^{\left(2\right)}+(4-3\mathrm{\&alpha;}){\hat{y}}_{n+1}^{\left(3\right)}\&rsqb;$

$c_n=\frac{{\mathrm{\&alpha;}}^2}{2{(1-\mathrm{\&alpha;})}^2}\&lsqb;{\hat{y}}_{n+1}^{\left(1\right)}-2{\hat{y}}_{n+1}^{\left(2\right)}+{\hat{y}}_{n+1}^{\left(3\right)}\&rsqb;$

其中：分别是一次、二次、三次指数平滑数，由如下公式计算得到：

${\hat{y}}_{n+1}^{\left(1\right)}=\mathrm{\&alpha;}\&times;y_n+(1-\mathrm{\&alpha;})\&times;{\hat{y}}_n^{\left(1\right)}$

${\hat{y}}_{n+1}^{\left(2\right)}=\mathrm{\&alpha;}\&times;{\hat{y}}_{n+1}^{\left(1\right)}+(1-\mathrm{\&alpha;})\&times;{\hat{y}}_n^{\left(2\right)}$

${\hat{y}}_{n+1}^{\left(3\right)}=\mathrm{\&alpha;}\&times;{\hat{y}}_{n+1}^{\left(2\right)}+(1-\mathrm{\&alpha;})\&times;{\hat{y}}_n^{\left(3\right)}$

是三次指数平滑法的初始值，其取值为

${\hat{y}}_0^{\left(1\right)}={\hat{y}}_0^{\left(2\right)}={\hat{y}}_0^{\left(3\right)}=\frac{y_1+y_2+y_3}{3}$

α是平滑系数(0＜α＜1)，体现信任的时间衰减特性，即离预测值越近的时间片的y_t权重越大，离预测值越远的时间片的y_t权重越小；如果数据波动较大，且长期趋势变化幅度较大，呈现明显迅速的上升或下降趋势时α应取较大值(0.6～0.8)，增加近期数据对预测结果的影响；当数据有波动，但长期趋势变化不大时，α可在0.1～0.4取值；如果数据波动平稳，α取0.05～0.20；

计算直接信任值：

节点j的直接信任值TD_ij为预测交互次数和真实交互次数y_n+1的相对误差，

采用多路径信任推荐方式而得到的计算式计算间接信任值；收集可信节点对节点j的直接信任值：节点i向所有满足TD_ik≤φ的可信关联节点询问其对节点j的直接信任值，其中φ为推荐节点的可信度阈值，根据可信度的要求精度，φ的取值范围为0～0.4；计算间接信任值：综合计算所收集到的信任值，得到节点j的间接信任值TR_ij，其中，Set(i)为观测节点i的关联节点中与j节点有过交互且其直接信任值满足TD_ik≤φ的节点集合；

由直接信任值和间接信任值整合计算得出综合信任值；综合信任值(T_ij)的计算公式如下：T_ij＝βTD_ij+(1-β)TR_ij，其中β(0≤β≤1)表示直接信任值的权重，当β＝0时，节点i和节点j没有直接交互关系，综合信任值的计算直接来自于间接信任值，判断较客观；当β＝1时，节点i对节点j的综合信任值全部来自于直接信任值，在这种情况下，判断较为主观，实际计算根据需要确定β的取值；

步骤二，对M路离散时域混合信号进行重叠加窗短时傅里叶变换，得到M个混合信号的时频域矩阵p＝0,1,…,P-1,q＝0,1,…,N_fft-1，其中P表示总的窗数，N_fft表示FFT变换长度；p，q)表示时频索引，具体的时频值为这里N_fft表示FFT变换的长度，p表示加窗次数，T_s表示采样间隔，f_s表示采样频率，C为整数，表示短时傅里叶变换加窗间隔的采样点数，C＜N_fft，且K_c＝N_fft/C为整数，也就是说采用的是重叠加窗的短时傅里叶变换；

步骤三，对步骤二中得到的跳频混合信号时频域矩阵进行预处理；

步骤四，利用聚类算法估计每一跳的跳变时刻以及各跳对应的归一化的混合矩阵列向量、跳频频率；在p(p＝0，1，2，…P-1)时刻，对表示的频率值进行聚类，得到的聚类中心个数表示p时刻存在的载频个数，个聚类中心则表示载频的大小，分别用表示；对每一采样时刻p(p＝0，1，2，…P-1)，利用聚类算法对进行聚类，同样可得到个聚类中心，用表示；对所有求均值并取整，得到源信号个数的估计即：

$\hat{N}=round\left(\frac{1}{p}\underset{p=0}{\overset{P-1}{\&Sigma;}}{\hat{N}}_p\right);$

找出的时刻，用p_h表示，对每一段连续取值的p_h求中值，用表示第l段相连p_h的中值，则表示第l个频率跳变时刻的估计；根据估计得到的p≠p_h以及第四步中估计得到的频率跳变时刻估计出每一跳对应的个混合矩阵列向量具体公式为：

${\hat{a}}_n\left(l\right)=\left\{\begin{array}{cc}\frac{1}{{\stackrel{\&OverBar;}{p}}_h\left(1\right)}\&CenterDot;\underset{p=1,p\&NotEqual;p_h}{\overset{{\stackrel{\&OverBar;}{p}}_h\left(1\right)}{\mathrm{\&Sigma;}}}{b}_{n,p}^0& l=1,\\ \frac{1}{{\stackrel{\&OverBar;}{p}}_h\left(l\right)-{\stackrel{\&OverBar;}{p}}_h\left(l-1\right)}\&CenterDot;\underset{p={\stackrel{\&OverBar;}{p}}_h\left(l-1\right)+1,p\&NotEqual;p_h}{\overset{{\stackrel{\&OverBar;}{p}}_h\left(l\right)}{\mathrm{\&Sigma;}}}{b}_{n,p}^0& l>1,\end{array}\right.,n=1,2,...,\hat{N};$

这里表示第l跳对应的个混合矩阵列向量估计值；估计每一跳对应的载频频率，用表示第l跳对应的个频率估计值，计算公式如下：

${\hat{f}}_{c,n}\left(l\right)=\left\{\begin{array}{cc}\frac{1}{{\stackrel{\&OverBar;}{p}}_h\left(1\right)}\&CenterDot;\underset{p=1,p\&NotEqual;p_h}{\overset{{\stackrel{\&OverBar;}{p}}_h\left(1\right)}{\mathrm{\&Sigma;}}}{f}_o^n\left(p\right)& l=1,\\ \frac{1}{{\stackrel{\&OverBar;}{p}}_h\left(l\right)-{\stackrel{\&OverBar;}{p}}_h\left(l-1\right)}\&CenterDot;\underset{p={\stackrel{\&OverBar;}{p}}_h\left(l-1\right)+1,p\&NotEqual;p_h}{\overset{{\stackrel{\&OverBar;}{p}}_h\left(l\right)}{\mathrm{\&Sigma;}}}{f}_o^n\left(p\right)& l>1,\end{array}\right.,n=1,2,...,\hat{N};$

步骤五，根据步骤四估计得到的归一化混合矩阵列向量估计时频域跳频源信号；

步骤六，对不同跳频点之间的时频域跳频源信号进行拼接；

步骤七，根据源信号时频域估计值，恢复时域跳频源信号；对每一采样时刻p(p＝0，1，2，…)的频域数据Y_n(p,q),q＝0,1,2,…,N_fft-1做N_fft点的IFFT变换，得到p采样时刻对应的时域跳频源信号，用y_n(p,q_t)(q_t＝0,1,2,…,N_fft-1)表示；对上述所有时刻得到的时域跳频源信号y_n(p,q_t)进行合并处理，得到最终的时域跳频源信号估计，具体公式如下：

$s_n\&lsqb;kC:(k+1)C-1\&rsqb;=\left\{\begin{array}{cc}\underset{m=0}{\overset{k}{\mathrm{\&Sigma;}}}{y}_n\&lsqb;m,(k-m)C:(k-m+1)C-1\&rsqb;& k<K_c\\ \underset{m=k-K_c+1}{\overset{k}{\mathrm{\&Sigma;}}}{y}_n\&lsqb;m,(k-m)C:(k-m+1)C-1\&rsqb;& k\&GreaterEqual;K_c\end{array}\right.,k=0,1,2,...$

这里K_c＝N_fft/C，C为短时傅里叶变换加窗间隔的采样点数，N_fft为FFT变换的长度。

3.如权利要求1所述的基于X波段多普勒雷达的降水估测方法，其特征在于，所述信号处理器设置有求解单元，所述求解单元的信号处理方法包括：

第一步，保证次级用户网络能长时间工作，需对次级用户的发射功率进行限制，保证次级用户网络的平均发射功率低于限定值：

E{α₀P₀+α₁P₁+β₀P₀+β₁P₁}≤P_av；

式中P_av是次级用户发射机SU-Tx的最大平均发射功率，这的平均是指信道衰减系数h_i,g_ss,g_sp随机变量的期望；

第二步，认知无线电网络的首要任务是保护主用户网络的服务质量，故对网络的干扰功率进行了限制；根据基于合作感知的频谱共享网络模型，知道干扰只在主用户PU处于忙状态时发生，所以平均干扰功率约束写成如下形式：

E{g_sp(β₀P₀+β₁P₁)}≤Q_av；

第三步，确保各个节点处的检测概率和网络的整体检测概率分别不低于各自的目标检测概率，关于检测概率的限制条件如下：

P_d≥P_th,P_di≥p_th,i＝1,2…k；

第四步，根据上述限制条件下，建立以最大化次级网络的平均吞吐量为目标函数的最优化问题：

$\begin{array}{cc}\underset{\left\{{\mathrm{\&tau;}}_s,\mathrm{\&epsiv;},\left\{{\mathrm{\&epsiv;}}_i\right\},P_0,P_1\right\}}{\max }& R\\ subjectto& \left(2\right),\left(3\right),\left(4\right),P_0\&GreaterEqual;0,P_1\&GreaterEqual;0\\ & 0\&le;{\mathrm{\&tau;}}_s\&le;T\end{array}---\left(\mathrm{Pr}oblem1\right)$

第五步，求解所建立的最优化问题，选择使得次级网络的吞吐量最大的合作感知的感知周期和次级用户的信号发射功率作为该频谱感知模型的感知参数；

具体包括以下的步骤：

1)，对不等式约束条件组(4)取等号，简化Problem 1为Problem 2；

$\begin{array}{cc}\underset{\left\{{\mathrm{\&tau;}}_s,P_0,P_1\right\}}{\max }& R\\ subjectto& \left(2\right),\left(3\right),P_0\&GreaterEqual;0,P_1\&GreaterEqual;0\\ & 0\&le;{\mathrm{\&tau;}}_s\&le;T\end{array}---\left(\mathrm{Pr}oblem2\right)$

2)，弱化对感知周期τ_s的求解，重点求解使平均吞吐量最大化的信号发射功率P₀，P₁；关于发射功率P₀和P₁的拉格朗日函数如下：

$\begin{array}{c}L\left(P_0,P_1,\mathrm{\&lambda;},\mathrm{\&mu;}\right)=E\left\{\frac{T-{\stackrel{\&OverBar;}{\mathrm{\&tau;}}}_s}{T}\&lsqb;{\mathrm{\&alpha;}}_0{r}_{00}+{\mathrm{\&alpha;}}_1{r}_{01}+{\mathrm{\&beta;}}_0{r}_{10}+{\mathrm{\&beta;}}_1{r}_{11}\&rsqb;\right\}\\ -\mathrm{\&lambda;}\&lsqb;E\left\{{\mathrm{\&alpha;}}_0{P}_0+{\mathrm{\&alpha;}}_1{P}_1+{\mathrm{\&beta;}}_0{P}_0+{\mathrm{\&beta;}}_1{P}_1\right\}-P_{av}\&rsqb;\\ -\mathrm{\&mu;}\&lsqb;E\left\{g_{sp}\left({\mathrm{\&beta;}}_0{P}_0+{\mathrm{\&beta;}}_1{P}_1\right)\right\}-Q_{av}\&rsqb;.\end{array}$

所以Problem 2的拉格朗日对偶优化问题为：

$\begin{array}{cc}\underset{\mathrm{\&lambda;}\&GreaterEqual;0,\mathrm{\&mu;}\&GreaterEqual;0}{\min imize}& g(\mathrm{\&lambda;},\mathrm{\&mu;})\end{array}---\left(\mathrm{Pr}oblem3\right)$

其中表示拉格朗日对偶函数；证明优化问题Problem 2与Problem 3的最优值差值为零，说明优化问题Problem 2与其拉格朗日对偶优化问题Problem 3之间是等价的，故只需求Problem 3的最优解即可；该问题是一个关于双变量P₀P₁的联合规划问题，为此将分解成两个子优化问题：

SP1：

SP2：

看出SP1和SP2分别是关于P₀P₁的无约束凸优化问题，此时运用拉格朗日函数以及KKT条件，便得到当检测到主用户PU处于闲状态时次级用户发射机SU-Tx的最优发射功率：

$P_0={\&lsqb;\frac{A_0+\sqrt{{\mathrm{\&Lambda;}}_0}}{2}\&rsqb;}^{+};$

其中：

${\mathrm{\&Lambda;}}_0=A_0^2-\frac{4}{g_{ss}}\left\{\frac{{\mathrm{\&sigma;}}_u^4+{\mathrm{\&sigma;}}_u^2{h}_k{P}_p}{g_{ss}}-\frac{{\log }_2\left(e\right)\&lsqb;{\mathrm{\&alpha;}}_0\left({\mathrm{\&sigma;}}_u^2+h_k{P}_p\right)+{\mathrm{\&beta;}}_0{\mathrm{\&sigma;}}_u^2\&rsqb;}{\mathrm{\&lambda;}\left({\mathrm{\&alpha;}}_0+{\mathrm{\&beta;}}_0\right)+{\mathrm{\&mu;\&beta;}}_0{g}_{sp}}\right\}$

当检测到主用户PU处于忙状态时，次级用户发射机SU-Tx的最优发射功率为：

$P_1={\&lsqb;\frac{A_1+\sqrt{{\mathrm{\&Lambda;}}_1}}{2}\&rsqb;}^{+};$

其中：

${\mathrm{\&Lambda;}}_1=A_1^2-\frac{4}{g_{ss}}\left\{\frac{{\mathrm{\&sigma;}}_u^4+{\mathrm{\&sigma;}}_u^2{h}_k{P}_p}{g_{ss}}-\frac{{\log }_2\left(e\right)\&lsqb;{\mathrm{\&alpha;}}_1\left({\mathrm{\&sigma;}}_u^2+h_k{P}_p\right)+{\mathrm{\&beta;}}_1{\mathrm{\&sigma;}}_u^2\&rsqb;}{\mathrm{\&lambda;}\left({\mathrm{\&alpha;}}_1+{\mathrm{\&beta;}}_1\right)+{\mathrm{\&mu;\&beta;}}_1{g}_{sp}}\right\}$

式中[x]⁺＝max{0,x}；λ≥0，μ≥0是式E{α₀P₀+α₁P₁+β₀P₀+β₁P₁}≤P_av和E{g_sp(β₀P₀+β₁P₁)}≤Q_av的拉格朗日乘子。