2. 河南省大气探测技术保障中心,郑州 450003
2. Henan Provincial Atmospheric Observation Technical Support Center, Zhengzhou 450003
风廓线雷达是利用大气湍流对电磁波的散射作用对大气风场等物理量进行探测的一种遥感设备。风廓线雷达观测网的建立弥补了常规高空探测站网空间密度和观测时次上的不足,在中小尺度灾害性天气的监测中发挥了重要作用(何平,2006;董丽萍等,2014;王令等,2014)。由于风廓线雷达回波信号非常微弱,具有明显的随机起伏涨落,且多普勒频移相对较小,因此,为保证风廓线雷达的定量测量精度,从随机涨落的回波信号中提取有用的气象信息,有必要对风廓线雷达系统进行标定。近年来,张沛源等(2001)和潘新民等(2002;2010a;2010b)介绍了新一代天气雷达回波强度标定原理,描述了利用测试信号,根据雷达气象方程计算出回波强度理论值和测量值的差值进行在线订正的方法。潘新民等(2010c)介绍了利用变化注入信号相位或频率对新一代天气雷达的速度进行标定的方法。这些方法对于风廓线雷达的回波强度和速度标定具有借鉴和指导意义,有利于减小风廓线雷达的测量误差,提高探测数据可靠性。
本文第一部分介绍了目前国内风廓线雷达回波强度的标定原理,借助天气雷达的标定方法,提出了在每个脉冲重复周期内,实时测量发射采样脉冲功率和回波功率在数字中频接收机的输出值的方法。采用该方法可以补偿发射机中固态功率放大器以及接收机中低噪声放大器、中频放大器等有源器件产生的增益漂移,改善回波强度的测量精度。本文第二部分介绍了风廓线雷达的速度标定原理,针对速度标定中经常出现的镜像谱分量现象建立了I/Q信号的数学模型,分析了I/Q信号幅相不平衡对速度测量的影响,并介绍了一种基于信号时域统计特性的订正方法,通过仿真实验证明该方法能够有效抑制镜像谱分量,对保证风廓线雷达的数据质量有一定的意义。
1 回波强度标定 1.1 标定原理与误差分析集中式风廓线雷达标定原理框图如图 1所示。为便于分析,定义参考平面1为天线端口,参考平面2为发射机输出端口,参考平面3为接收机低噪声放大器(LNA)输入端口,参考平面4为数字中频接收机输出端口。在实际标定过程中,回波功率的测量值往往是在数字中频接收机输出端而非天线端口获得的,因此需要把接收机增益和由于接收机有限带宽而造成的功率损失考虑在内。由雷达气象方程(Bringi et al, 2001;Doviak et al, 1993),在参考平面4测量的回波功率Pro为
$ {P_{ro}} = \frac{{C\eta }}{{r_0^2}} $ | (1) |
式中,r0为气象目标的距离,η为反射率,在风廓线雷达中定义为
$ \eta = 0.38C_n^2{\lambda ^{ - \frac{1}{3}}} $ | (2) |
式中,Cn2为大气折射率结构常数。式(1) 中雷达常数
$ C = \left( {\frac{c}{{1024{{\rm{ \mathsf{ π} }}^{\rm{2}}}\ln 2}}} \right)\left( {{\lambda ^2}{P_t}\tau } \right)\left( {G_0^2\theta \phi } \right)\left( {\frac{{{G_r}}}{{{l_r}}}} \right) $ | (3) |
由四项组成,其中第一项为常数项,c为真空中的光速;第二项为发射机参数,λ为发射波长,Pt为定义在参考平面1的发射脉冲功率,τ为发射脉冲宽度;第三项为实测的天线参数,G0为天线增益,θ、ϕ分别为天线H面和E面波束宽度;第四项为接收机参数,Gr定义为从参考平面1经环形器、T/R开关至参考平面4的接收机增益,lr为接收机有限带宽损失因子。
在式(3) 中,λ、τ、G0、θ、ϕ、lr在雷达工作中基本保持不变,因此影响回波强度测量精度的因素主要是Pt和Gr。风廓线雷达采用固态发射机和数字中频接收机体制,受温度变化影响较大。当环境温度、电源电压或电路工作参数发生变化时,有源微波器件如固态功率放大器和低噪声放大器等容易产生增益漂移,表现为Pt和Gr的变化引起回波功率Pro的变化,最终导致反射率η的测量误差,而这类误差往往难以在雷达产品数据中发现。
传统的回波强度标定方法是从参考平面3注入固定功率的连续波(CW)测试信号,根据实测发射功率和固定距离库计算出回波强度理论值和测量值的差值进行在线订正,即用前一个扫描周期标定出的回波强度测量误差订正后一个扫描周期的回波强度测量值。对于业务中采用五波束扫描的风廓线雷达,扫描周期为5 min,如果在扫描周期内,CW测试信号功率和接收机增益发生变化,则有可能对回波强度进行错误订正,导致回波强度测量误差增大。
1.2 提高回波强度标定精度的方法 1.2.1 补偿发射功率变化对于发射功率变化,可以采用硬件和软件两种方法补偿。硬件上,在固态功率放大器中采用增益-温度补偿技术,当环境温度变化时,调整固态功率放大器的输出电压和电流,减小增益漂移,保持发射功率电平稳定在±1 dB之内。软件上,可以通过在线监测发射功率,当监测到发射功率有变化时,在线修正发射功率的测量值,保证回波强度的测量精度。
1.2.2 补偿接收机增益变化由于目前国内的风廓线雷达接收机中LNA端口注入的测试信号功率尚无实时监测功能,因此一旦测试信号功率或接收机增益发生变化,将使得回波强度测量值得到错误校正,直接导致回波强度测量误差。为解决这一问题,可以利用收发通道内的实时标定环路来补偿接收机增益变化,方法如下所述。在每个脉冲重复周期内,将发射机输出信号经过耦合、衰减,从参考平面3注入接收机中,将该采样信号按照时间顺序加以记录并保存。如图 1所示,在参考平面4,采样脉冲的功率测量值S1可以表示为
$ {S_1} = {P_t} + {L_{21}} - {L_{23}} + {G_{r1}} $ | (4) |
而后向散射的实际大气回波功率的测量值S2可以表示为
$ {S_2} = {P_r} - {L_{13}} + {G_{r1}} $ | (5) |
式(4)、式(5) 中各参数均已换算成dB单位,发射功率Pt和回波功率Pr都是在参考平面1的测量值,Gr1为从参考平面3到参考平面4的增益,用Lij表示从参考平面i到参考平面j之间的损耗值(Lij为正数),则L21为从参考平面2到参考平面1的损耗,L23为从参考平面2经耦合器、衰减器到参考平面3的损耗,L13为从参考平面1经环形器、T/R开关到参考平面3的损耗。目前国内风廓线雷达的脉冲重复周期范围为20~100 μs,在如此短暂的脉冲重复周期内,可以认为Gr1基本不变。由式(1)、(3) 可知,
$ \frac{{{P_r}}}{{{P_t}}} = \left( {\frac{c}{{1024{{\rm{ \mathsf{ π} }}^2}{\rm{ln}}2}}} \right)\left( {{\lambda ^2}\tau } \right)\left( {G_0^2\theta \phi } \right)\frac{\eta }{{r_0^2}} = \frac{{{C_1}\eta }}{{r_0^2}} $ | (6) |
C1为不含发射功率的雷达常数,由式(6),η与Pr/Pt成正比。若将式(4) 与式(5) 相减,可得
$ {P_r} - {P_t} = {S_2} - {S_1} + {L_{13}} + {L_{21}} - {L_{23}} $ | (7) |
式(7) 中各参数均采用dB单位。
不难发现,Gr1从式(4)、式(5) 中对消掉,将式(7) 代入式(6),即可得到
$ \begin{array}{l} \eta \left({{\rm{dB}}} \right) = \left({{S_2} - {S_1} + {L_{13}} + {L_{21}} - {L_{23}}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;20{\rm{lo}}{{\rm{g}}_{10}}{r_0} - 10{\rm{lo}}{{\rm{g}}_{10}}{C_1} \end{array} $ |
由于式(8) 中已不含发射功率和接收机增益项,因此,采用该方法可以补偿发射机中固态功率放大器以及接收机中低噪声放大器、中频放大器等有源器件产生的增益漂移,保证了回波强度的测量精度。
2 速度标定 2.1 标定原理风廓线雷达速度标定主要采用变化注入信号相位(简称“移相法”)或频率(简称“频偏法”)的方法。当采用移相法时,将图 1中内置信号源输出的CW测试信号经移相器变化相位Δφ后从参考平面3注入接收机,需要注意的是,由于风廓线雷达采用时域相干积累技术,因此采样时间Ts=Ncoh·PRT,Ncoh为时域相干积累次数,PRT为脉冲重复周期,则速度的理论值为
$ V = \frac{\lambda }{{4{\rm{ \mathsf{ π} }}}} \cdot \frac{{\Delta \varphi }}{{{T_s}}} = \frac{\lambda }{{4{\rm{ \mathsf{ π} }}}} \cdot \frac{{\Delta \varphi }}{{{N_{{\rm{coh}}}} \cdot PRT}} $ | (9) |
当采用频偏法时,利用机外信号源输出频率为fc+fd的测试信号,从参考平面3注入接收机,fc为雷达工作频率,fd为频偏量,则速度的理论值为V=(λ·fd)/2。
2.2 速度标定中的镜像谱现象分析在风廓线雷达速度标定时,功率谱中经常会出现与真实速度关于零频对称的镜像谱分量,如图 2所示,真实速度为v=4.04 m·s-1,但在v=-4.04 m·s-1处出现了镜像谱分量,下面就对该镜像谱分量产生的原因做一下理论分析。
理想情况下,数字中频接收机输出的等幅正交I/Q信号(Skolnik, 1990)分别为
$ {I_0}\left( n \right) = {\rm{cos}}\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right) $ | (10) |
$ {Q_0}\left( n \right) = {\rm{sin}}\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right) $ | (11) |
式中,fd代表大气运动的多普勒频移。为简化分析,将I/Q信号的幅度作归一化处理,且相位为零;并以I通道为基准,令幅度和相位不平衡均出现在Q通道中。实际情况下,数字中频接收机输出的I/Q信号分别为
$ {I_1}\left( n \right) = {\rm{cos}}\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right) + {\beta _{\rm{I}}} $ | (12) |
$ {Q_1}\left( n \right) = \alpha {\rm{sin}}\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t + \varphi } \right) + {\beta _{\rm{Q}}} $ | (13) |
式中,α、β分别代表正交信号Q1(n)相对于I1(n)的幅度和相位误差(对于幅度相位平衡的I、Q信号,有α=1,φ=0),βI和βQ分别为I、Q通道的直流偏移量。合成的复电压信号为
$ \begin{array}{*{20}{l}} {x\left( n \right) = {I_1}\left( n \right) + j{Q_1}\left( n \right) = {\rm{cos}}\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right) + }\\ {\;\;\;\;\;\;\;\;\;j\alpha {\rm{sin}}\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t + \varphi } \right) + {\beta _{\rm{I}}} + j{\beta _{\rm{Q}}}} \end{array} $ | (14) |
由Euler公式,有
$ \begin{array}{*{20}{l}} {{\rm{cos}}\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right) = \frac{{{{\rm{e}}^{j\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right)}} + {{\rm{e}}^{ - j\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right)}}}}{2},}\\ {{\rm{sin}}\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right) = \frac{{{{\rm{e}}^{j\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right)}} + {{\rm{e}}^{ - j\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right)}}}}{{2j}},} \end{array} $ |
则
$ \begin{array}{*{20}{l}} {x\left( n \right) = \frac{{1 + \alpha {\rm{cos}}\varphi + j\alpha {\rm{sin}}\varphi }}{2} \cdot {{\rm{e}}^{j\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right)}} + }\\ {\;\;\;\;\frac{{1 - \alpha {\rm{cos}}\varphi + j\alpha {\rm{sin}}\varphi }}{2} \cdot {{\rm{e}}^{ - j(2{\rm{ \mathsf{ π} }}{f_d}n\Delta t)}} + {\beta _{\rm{I}}} + j{\beta _{\rm{Q}}}} \end{array} $ | (15) |
由式(15) 可知,大气湍流回波经过幅度相位不平衡的I、Q通道后,输出的信号功率谱中除了真实的谱分量
直流偏移量为
$ \begin{array}{*{20}{l}} {{I_2}\left( n \right) = {\rm{cos}}\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t} \right),}\\ {{Q_2}\left( n \right) = \alpha {\rm{sin}}\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t + \varphi } \right)} \end{array} $ |
I2(n)、Q2(n)与I0(n)、Q0(n)的关系可以通过矩阵C来描述(Ellingson, 2003),即
$ \begin{array}{l} \left[ \begin{array}{l} {I_2}\left(n \right)\\ {Q_2}\left(n \right) \end{array} \right] = \mathit{\boldsymbol{C}}\left[ \begin{array}{l} {I_0}\left(n \right)\\ {Q_0}\left(n \right) \end{array} \right] = \\ \;\;\;\left[ \begin{array}{l} \;\;\;1\;\;\;\;\;\;\;\;\;\;\;\;\;\;0\\ \alpha {\rm{sin}}\varphi \;\;\;\alpha {\rm{cos}}\varphi \end{array} \right]\left[ \begin{array}{l} {I_0}\left(n \right)\\ {Q_0}\left(n \right) \end{array} \right] \end{array} $ | (16) |
式(16) 等式两边同乘以C-1,得
$ \begin{array}{l} \left[ \begin{array}{l} {I_0}\left(n \right)\\ {Q_0}\left(n \right) \end{array} \right] = {\mathit{\boldsymbol{C}}^{ - 1}}\left[ \begin{array}{l} {I_2}\left(n \right)\\ {Q_2}\left(n \right) \end{array} \right] = \\ \;\;\;\left[ \begin{array}{l} \;\;\;\;1\;\;\;\;\;\;\;\;\;\;0\\ - {\rm{tg}}\varphi \;\;\;\frac{1}{{\alpha {\rm{cos}}\varphi }} \end{array} \right]\left[ \begin{array}{l} {I_2}\left(n \right)\\ {Q_2}\left(n \right) \end{array} \right] \end{array} $ | (17) |
只要找到α和φ,就可以获得订正矩阵C-1中的未知元素。由
$ \begin{array}{l} \frac{1}{N}\sum\limits_{n = 1}^N {Q_2^2\left( n \right) = \frac{1}{N}} \sum\limits_{n = 1}^N {{\alpha ^2}{\rm{si}}{{\rm{n}}^2}\left( {2{\rm{ \mathsf{ π} }}{f_d}n{\rm{ }}\Delta {\rm{ }}t + \varphi } \right)} = \\ \frac{{{\alpha ^2}}}{2} \cdot \frac{1}{N}\sum\limits_{n = 1}^N {\left\{ {1 - {\rm{cos}}\left[ {2\left( {2{\rm{ \mathsf{ π} }}{f_d}n\Delta t + \varphi } \right)} \right]} \right\}} = \frac{{{\alpha ^2}}}{2} \end{array} $ | (18) |
可得
$ \alpha =\sqrt{\frac{2}{N}\sum\limits_{n=1}^{N}{Q_{2}^{2}\left(n \right)}} $ | (19) |
同理,由
$ \frac{1}{N}\sum\limits_{n=1}^{N}{\left[ {{I}_{2}}\left(n \right){{Q}_{2}}\left(n \right) \right]}=\frac{\alpha \text{sin}\varphi }{2} $ | (20) |
可得
$ \varphi =\text{arcsin}\left\{ \frac{2}{N\alpha }\sum\limits_{n=1}^{N}{\left[ {{I}_{2}}\left(n \right)Q_2\left(n \right) \right]} \right\} $ | (21) |
将式(19)、式(21) 代入式(17),可得订正矩阵C-1中的未知元素。在信号处理器中利用订正矩阵C-1对I/Q数据进行订正,即可抑制甚至消除速度标定中的镜像谱分量。
2.4 仿真实验在MATLAB软件中设置仿真参数如下:风廓线雷达工作频率1290 MHz,采样频率fs=1000 Hz,FFT点数256,I/Q数据量为4096,多普勒频率100 Hz,I2(n)=cos(2π·100·nΔt),Q2(n)=0.9sin(2π·100·nΔt+0.5),则多普勒速度理论值V=11.6 m·s-1。如图 3所示,订正前在-11.6 m·s-1处存在镜像谱分量,订正后镜像谱分量消失,证明该订正方法是有效的。
本文研究了风廓线雷达标定系统中的两个问题。首先,介绍了风廓线雷达回波强度的标定原理,根据雷达气象方程,分析了发射功率变化和接收机增益变化对回波强度测量精度的影响,并借助天气雷达的标定方法,提出了通过实时测量每个脉冲重复周期内发射采样脉冲功率和回波功率在数字中频接收机的输出值,改善回波强度测量精度的方法。其次,介绍了风廓线雷达速度标定原理,分析了I/Q信号的幅度、相位不平衡对速度测量的影响,提出了一种基于信号时域统计特性的订正方法,并利用I/Q信号数学模型进行了仿真实验。结果表明,该方法可以有效抑制I/Q信号幅相不平衡产生的镜像谱分量。由于目前试验条件限制,实际风廓线雷达精确标定结果将在下一步工作中完善。
董丽萍, 吴蕾, 王令, 等, 2014. 风廓线雷达组网资料初步对比分析[J]. 气象, 40(9): 1145-1151. DOI:10.7519/j.issn.1000-0526.2014.09.012 |
何平, 2006. 相控阵风廓线雷达[M]. 北京: 气象出版社.
|
潘新民, 柴秀梅, 崔柄俭, 等, 2010a. CINRAD/SB雷达回波强度定标调校方法[J]. 应用气象学报, 21(6): 739-746. |
潘新民, 柴秀梅, 黄跃青, 等, 2010b. CINRAD/SA & SB回波强度定标故障的诊断分析和解决方法[J]. 气象, 36(12): 122-127. |
潘新民, 柴秀梅, 徐俊领, 等, 2010c. 新一代天气雷达测速定标精度检查方法[J]. 气象科技, 38(2): 214-221. |
潘新民, 汤志亚, 2002. 天气雷达接收功率标定的检验方法探讨[J]. 气象, 28(4): 34-37. DOI:10.7519/j.issn.1000-0526.2002.04.007 |
王令, 王国荣, 古月, 等, 2014. 风廓线雷达垂直径向速度应用初探[J]. 气象, 40(3): 290-296. DOI:10.7519/j.issn.1000-0526.2014.03.004 |
张沛源, 周海光, 梁海河, 等, 2001. 数字化天气雷达定标中应注意的一些问题[J]. 气象, 27(6): 27-32. DOI:10.7519/j.issn.1000-0526.2001.06.005 |
Bringi V N, Chandrasekar V, 2001. Polarimetric Doppler Weather Radar:Principles and Applications[M].
Cambridge: Cambridge University Press, 294-316.
|
Doviak R J, Zrnic D S, 1993. Doppler Radar and Weather Observations[M].
San Diego: Academic Press, 145-158.
|
Ellingson S W. 2003. Correcting I-Q Imbalance in Direct Conversion Receivers.
|
Skolnik M I, 1990. Radar Handbook[M].
Boston: McGraw-Hill Press.
|
Vega M A, Chandrasekar V, Nguyen C, et al, 2012. Calibration of the NASA Dual-Frequency, Dual-Polarized, Doppler Radar. 2012 IEEE International Geoscience and Remote Sensing Symposium, Munich[J]. IEEE Geoscience and Remote Sensing Societ: 4625-4628. |