全文HTML
--> --> -->除了近地面传播的视距微波链路以外, 地面还广泛覆盖了卫星和地面站之间传播的星地链路[13]. 高频星地链路在大气中传播时同样会受到降水的衰减, 因此同样也可以用于降雨测量. 李黄[14]首先提出利用Ku波段卫星通信雨衰进行降雨探测的设想, 并开展初步实验. Arslan等[15]利用Ku波段地球同步轨道(GEO)广播卫星信号测量降雨, 其结果与雷达数据相比具有较高的一致性. Francois等[16]搭建4条GEO卫星链路并联合四维变分同化的方法对近地面区域降雨进行反演, 证明星地链路探测区域降雨的可行性. Adirosi等[17]使用雨量筒和激光雨滴谱仪对星地链路测雨结果进行定标, 提高了降雨探测的准确性.
现有研究大多针对GEO卫星, 卫星和地面站的位置相对固定, 仅能够探测路径平均降雨. 此外, 未来还将发射多个低轨(LEO)或中轨(MEO)的宽带通信卫星星座[13]. 与GEO卫星不同, 在LEO和MEO卫星与地面站通信过程中, 随着卫星的移动, 星地链路像X光扫描人体一样对降雨区进行“扫描”, 并且由于星地链路倾斜穿过整层降雨区域, 其衰减中还包含各个高度上的降雨信息. 利用计算机断层成像技术等方法就可以实现对二维垂直降雨场的重建. 因此可以成为一种测量降雨垂直分布的新方法. 星地链路在搭建后可自主运行, 实测数据可通过卫星回传, 并且卫星星座覆盖全球, 因此在高原、山区、海岛等传统地面观测资料缺失的地区使用具有独特优势. 本文在建立星地链路大气衰减定量模型的基础上, 提出并研究基于联合代数重建技术的星地链路反演二维垂直降雨场方法. 利用实测降雨资料开展仿真分析, 验证星地链路反演垂直降雨场的可行性和有效性, 为发展星地链路测雨新方法奠定技术基础.
2.1.星地链路几何结构
卫星从位置A移动到位置C时, 地球站持续接收信号过程如图1所示, 其中, O代表地心, E为卫星轨道正下方地球站位置, h和R分别表示轨道高度和地球平均半径,

Figure1. Geometry of earth-space link.


假设卫星从位置A移动到位置B时地球站位于降雨区域内, 星地链路探测垂直降雨场的过程如图2所示, 图中hrain表示雨区相对于地球站的垂直高度, Lr和Lo分别代表星地链路在雨区内和雨区外的传播距离, Lo远大于Lr, 则星地链路的水平探测尺度W为

Figure2. The detection of vertical rainfall field with earth-space link.
2
2.2.星地链路信号衰减
在星地通信系统中, 地球站接收功率Pr (dB)可以表示为[18]自由空间衰减、闪烁衰减、气体衰减和云致衰减的计算可参考ITU-R P.618建议书[18], 在南京地区接收AsiaSat5卫星Ku波段信号, 电磁波为垂直极化, 卫星轨道高度为36000 km, 天线仰角为43°, 统计高层云底高为3.55 km, 云顶高为4.77 km[19], 则在12—18 GHz范围内, 星地链路上的自由空间衰减、闪烁衰减、整层大气衰减、云致衰减以及10 mm/h和20 mm/h雨强时的降雨衰减如图3所示. 整体而言, 气体衰减最小在10–1—1 dB范围内, 闪烁衰减和云致衰减次之在1 dB左右, 降雨引起的衰减明显, 并且随着降雨强度和频率的增加而增加, 能够达到10 dB以上, 高于气体衰减、闪烁衰减及云致衰减而不会被其淹没. 虽然自由空间衰减超过102 dB, 但由于其仅是电波频率和传播距离的函数[18]并且当电波频率偏移200 MHz, 卫星轨道高度变化1000 km时自由空间衰减的相对变化量不足0.2%, 因此在已知卫星工作频率和轨道的条件下能够精确计算自由空间衰减, 可以将其从总衰减中剔除, 消除对降雨衰减测量的影响.

Figure3. Scintillation attenuation, cloud attenuation, gas attenuation and rain attenuation with 10 mm/h and 20 mm/h rain rate in earth space link.
3.1.垂直降雨场反演模型
星地链路在某一条路径上的降雨总衰减为为实现垂直降雨场的反演, 首先将探测区域划分成大小相等的方形网格, 假设每个网格内衰减系数均匀, 降雨场网格化示意图如图4所示. 其中


Figure4. The sketch of calculation of earth-space link rain attenuation.




雨衰系数γR与雨强R的关系满足幂律关系[20]:
由ITU-R P.838-3建议书[20]发现, 当星地链路频率为垂直极化17 GHz (Ku波段)时, 幂律系数随仰角的变化很小(图5), α和β均在0.04范围内变化, 两者的标准差(standard deviation, STD)和标准差系数(coefficient of variation, CV)分别为0.0012, 1.74%和0.0136, 1.32%, 并且β接近于1, 雨强和衰减系数近似呈线性关系. 因此, 本文在反演降雨强度时幂律系数选择为



Figure5. The change of power law coefficients with elevation.
2
3.2.基于联合代数重建技术(SART)的降雨场反演算法
SART作为一种CT成像技术中巨型稀疏矩阵求解的经典算法, 能够有效地解决方程组超定或欠定问题[21], 是Kaczmarz[22]的迭代重建技术 (algebraic reconstruction technique, ART)的进一步发展. SART算法的优势在于每次迭代过程中所有的方程都被用于像素解的更新, 实现噪声的平滑处理, 比ART算法的稳定性更高[21]. 基于SART算法, 垂直降水场衰减系数向量γ的迭代过程为[23]





由于衰减系数取值为正数, 所以需在迭代过程中要加入非负约束. 定义投影函数





2
3.3.反演算法流程
反演算法流程图如图6所示, 具体的反演过程为
Figure6. The flow chart of vertical rainfall field inversion method.
1)设置最大迭代次数Kmax和误差阈值

2)初始化γ0(很小的随机向量或零向量);
3)由(22)式和(23)式计算矩阵Dc, Dr;
4)由(25)式迭代求解衰减系数;
5) k = k + 1;
6)


7)由(18)式计算降雨强度.
4.1.初始场参数设置
利用南京地区2014年6月30日、7月2日和10月28日的微雨雷达(micro rain radar, MRR)实测资料, 将时间尺度转为水平距离尺度, 构建3类二维垂直降雨场I, II和III, 如图7所示. 可以清楚地看到降雨垂直方向上分布的不均匀性. 其中, 降雨场I中5 km高度上出现3个强降雨中心, 并且在水平5—10 km范围上空出现连续性垂直降雨带, 最大降雨强度为8 mm/h; 在降雨场II中5 km高度上出现1个强降雨中心, 最大降雨强度为18 mm/h; 在降雨场III中3—4 km高度范围内出现连续性水平降雨区, 最大降雨强度为40 mm/h.
Figure7. Vertical rainfall field derived from the data of MRR: (a) rainfall field I; (b) rainfall field II; (c) rainfall field III.
将上述垂直降雨场划分成1个31 × 31方形组成的网格, 其垂直分辨率为200 m/格, 水平分辨率为1 km/格, 并由位于(–10 km, 0 km)的地球站1, (64 km, 0 km)的地球站2和(15 km, 0 km)的地球站3进行探测, 如图8所示, 图中θ1, θ2和θ3分别表示3个地球站天线的仰角, Δθ1, Δθ2和Δθ3分别表示3个天线相邻两次探测的间隔仰角. 参考Ku工作波段AsiaSat5卫星的EIRP = 56 dB, 直径为1.5 m抛物面接收天线增益Gr = 49 dB, 则增益常数C = 105 dB. 星地链路具体参数如表1所示.
地球站 | 频率/ GHz | 极化方式 | Δθ/(°) | θMin/(°) | C/dB |
1 | 17 | 垂直 | 0.1 | 0.091 | 105 |
2 | 17 | 垂直 | 0.1 | 0.065 | 105 |
3 | 17 | 垂直 | 0.1 | 1.00 | 105 |
表1星地链路参数
Table1.Parameters of earth-space link.

Figure8. Sketch of grid for rainfall field and detection of earth-space links.
在MRR实测资料的基础上开展仿真实验, 由(17)式得到每个网格内的衰减系数, 根据表1中星地链路参数得到探测过程中的降雨衰减并建立线性方程组(15), 然后基于带有非负限制的SART反演算法对方程组求衰减系数解, 最后结合(18)式由衰减系数得到降雨强度, 从而重构垂直降雨场.
2
4.2.反演结果分析
首先, 用地球站1的探测数据来反演垂直降雨场, 结果如图9(a)—(c)所示. 可以看到, 反演场的降雨强度及其分布与真实场相比存在明显的差异. 其中, 降雨场I, II和III与其反演场的相关系数、平均偏差分别为0.556, 0.630 mm/h; 0.504, 0.865 mm/h和0.364, 2.553 mm/h. 3个反演场欧式距离和熵随迭代次数的变化如图10(a)—(c)和图11(a)—(c)所示. 从图10(a)—(c)可以看到, 虽然整个迭代过程的反演场逐渐趋向于真实场, 但经过500次迭代后的欧式距离分别为1.046, 1.429和4.103 mm/h, 与真实场的差距较大. 图11(a)—(c)中, 黑色虚线代表实际降雨场的熵, 在整个反演过程中反演场熵的变化并不明显, 表明反演效果随迭代次数未发生显著的变化, 500次迭代后反演场熵和真实场熵的相对误差分别为11.35%, 3.65%和3.25%. 改用地球站2或地球站3的探测数据后也出现了类似的情况. 之后又尝试减小θmin和Δθ以得到更多的探测数据, 但反演效果并没有明显的改善. 实验结果表明, 仅利用1个地球站难以实现垂直降雨场的反演. 这是因为在探测过程中出现了多条相邻星地链路经过同一组网格的现象, 导致方程组(15)的有效方程数量即矩阵L的秩减少, 即使设置Δθ = 0.05°使矩阵L的维度达到1441 × 962, 但有效方程的数量也仅为572, 同时由于降雨场网格数量庞大, 使得SART算法在迭代过程中易趋向局部最优解.
Figure9. Inversed vertical rainfall field: (a) Inversed field I with one earth station; (b) inversed field II with one earth station; (c) inversed field III with one earth station; (d) inversed field I with two earth stations; (e) inversed field II with two earth stations; (f) inversed field III with two earth stations; (g) inversed field I with three earth stations; (h) inversed field II with three earth stations; (i) inversed field III with three earth stations.

Figure10. The change of Euclidean distance with iterations: (a) Inversed field I with one earth station; (b) inversed field II with one earth station; (c) inversed field III with one earth station;(d) inversed field I with two earth stations; (e) inversed field II with two earth stations; (f) inversed field III with two earth stations;(g) inversed field I with three earth stations; (h) inversed field II with three earth stations; (i) inversed field III with three earth stations.

Figure11. The change of entropy with iterations: (a) inversed field I with one earth station; (b) inversed field II with one earth station; (c) inversed field III with one earth station; (d) inversed field I with two earth stations; (e) inversed field II with two earth stations; (f) inversed field III with two earth stations; (g) inversed field I with three earth stations; (h) inversed field II with three earth stations; (i) inversed field III with three earth stations.
基于第1次的实验结果, 本文联合两个地球站1和地球站2对降雨场进行探测, 反演结果示于图9(d)—(f). 可知, 垂直降雨场的整体反演效果较好, 能够清楚地看到降雨场I中5 km高度上的3个强降雨中心和连续性垂直降雨带, 降雨场II中5 km上的1个强降雨中心和降雨场III中的连续性水平降雨区, 并且准确得到每个降雨场的最大降雨强度. 但是反演场I中连续性垂直降雨带内部的降雨分布、反演场II中3 km高度上的降雨中心以及反演场III中3 km高度以下的弱降水中心仍与真实场存在偏差. 其中降雨场I, II, III与其反演场相关系数、平均偏差分别为0.980, 0.122 mm/h; 0.989, 0.159 mm/h和0.982, 0.537 mm/h. 欧式距离和熵的变化情况如图10(d)—(f)和图11(d)—(f)所示, 经过500次迭代后3个反演场的欧式距离、熵的相对误差分别为0.246 mm/h、1.53%, 0.235 mm/h、0.061%和0.812 mm/h、0.23%. 虽然实验结果与真实场仍存在部分微小偏差, 但两个地球站已基本实现基于星地链路的垂直降雨场反演. 与单地球站相比, 在本次实验中从两个不同的位置对降雨场进行探测, 使有效方程的数量增加到944个接近与降雨场的网格数量, SART算法在迭代过程中能够更加准确地趋向全局最优解.
在前两次实验的基础上, 本文联合3个地球站对降雨场进行探测, 反演结果如图9(g)—(i)所示, 与基于双地球站的反演结果相比, 前者存在的偏差在本次实验结果中均已消失, 反演场更加接近于真实场, 其中降雨场I, II, III与其反演场的相关系数接近1, 平均偏差分别为4.22 × 10–12, 2.65 × 10–12和3.64 × 10–12 mm/h. 从图10(g)—(i)可以看到欧式距离均小于0.01 mm/h, 图11(g)—(i)中反演场熵的相对误差均在0.01%以下. 实验结果表明, 基于3个地球站已实现垂直降雨场的准确反演. 同时, 在本次实验中, 有效方程的数量为961, 与降雨场网格数量相同, 在反演过程中SART算法能精准地找到全局最优解.
1)对于非降雨因素引起的衰减, 气体衰减最小在10–1—1 dB范围内, 闪烁衰减和云致衰减次之在1 dB左右, 明显小于降雨衰减而不会将其淹没. 自由空间衰减最严重为102 dB高于其他衰减的1—3个数量级, 但由于其仅是电波频率和传播距离的函数, 并且当电波频率变化为200 MHz, 卫星轨道变化为1000 km时自由空间衰减的相对变化量不足0.2%, 因此在已知卫星工作频率和轨道的条件下可以对自由空间衰减进行精确计算, 以将其从总衰减中剔除从而降低对降雨衰减测量的影响.
2)在卫星信号频率为17 GHz, 极化方式为垂直极化的条件下, 仅利用1个地球站难以实现对垂直降雨场的反演, 3个反演场与真实场之间的相关系数分别为0.556, 0.504和0.364, 平均偏差分别为0.630, 0.865和2.522 mm/h. 这是因为在探测过程中出现了多条相邻星地链路经过同一组网格的现象, 导致方程组的有效方程数量即矩阵L的秩减少, 同时由于降雨场网格数量庞大, 使得SART算法在迭代过程中易趋向局部最优解.
3)地球站数量增加到两个的情况下, 反演效果得到显著提升, 此时反演结果与真实降雨场间的相关系数分别为0.980, 0.989和0.982, 平均偏差分别为 0.122, 0.159和0.537 mm/h, 基本实现了基于星地链路的垂直降雨场反演; 当地球站数量上升到3个时, 反演结果更加精确, 此时相关系数均接近1, 真实降雨场被准确反演.
通过开展数值仿真实验证明了基于星地链路的二维垂直降雨场反演方法的可行性, 下一步在利用实际链路探测降雨时, 还需在以下方面做进一步研究: 1)考虑实际链路的信号噪声, 选用合适的信号去噪方法确保降雨衰减信息的真实有效; 2)根据实际链路的密度和结构, 选择合适的降雨场网格格距以提升反演效果.