摘要：2023年7月底，受台风“杜苏芮”“卡努”及地形等因素综合影响，北京及周边发生罕见的极端暴雨，造成严重的人员伤亡与经济损失；大气可降水量（PWV）作为影响降雨的主要因素之一，探究暴雨发生及发展过程与PWV的关系，对于进一步建立暴雨预警模型有重大意义。本文选取U

本文内容来源于《测绘学报》2024年第1期（审图号 GS京(2025)0114号）

基于GNSS观测的2023北京特大暴雨分析

1.中国矿业大学（北京）地球科学与测绘工程学院，北京 100083

2.

3.

4.

基金项目

作者简介

第一作者：杨飞（1991—），男，博士，副教授，研究方向为GNSS高精度数据处理及GNSS气象学。 E-mail：

通讯作者: 李志才 E-mail：zcli@cumtb.edu.cn

摘要

2023年7月底，受台风“杜苏芮”“卡努”及地形等因素综合影响，北京及周边发生罕见的极端暴雨，造成严重的人员伤亡与经济损失；大气可降水量（PWV）作为影响降雨的主要因素之一，探究暴雨发生及发展过程与PWV的关系，对于进一步建立暴雨预警模型有重大意义。本文选取UTC 2023年7月25日至8月1日北京及周边地区34个GNSS站点、34个气象站、1个探空站点和ERA5数据，采用GAMIT 10.71反演高精度GNSS-PWV；提出了一种改进的插值算法并获取了本次极端暴雨期间北京及周边地区高时空分辨率的PWV空间数据；利用探空与ERA5数据，从多个角度对反演的PWV精度进行了详细评估；最后结合气象站降雨量数据，从时间和空间的角度分析了PWV变化与极端降雨之间的关系，并首次探讨了对流层延迟梯度和降雨走势的关系。结果表明，GNSS-PWV与RS-PWV的相关系数达到0.99，RMSE和bias约为0.52 mm和-0.52 mm；GNSS-PWV与ERA5-PWV的RMSE小于6 mm，bias范围为-4~1.5 mm；GNSS-PWV格网数据与ERA5反演的格网结果RMSE和bias均值分别为4 mm和1 mm。通过分析PWV与并址气象站降水观测的时间序列，发现PWV在此次降雨前极速上升，在降水过程中仍累积上升，降水结束后无法立刻消散的特性，这与此次降水接连受台风“杜苏芮”和“卡努”影响相关；其次，各站点处对流层延迟梯度展现出指向东北方向的一致性，这与PWV高值在空间上呈现由西南向东北的输送趋势一致，与实际降水路线相符。

关键词

GNSS气象学大气可降水量降雨量

本文引用格式

杨飞, 汪莹莹, 李志才, 余博尧, 武军郦, 曹云昌, 张澍.

YANG Fei, WANG Yingying, LI Zhicai, YU Boyao, WU Junli, CAO Yunchang, ZHANG Shu.

阅读全文

2023年7月底至8月初，京津冀地区发生极端罕见的特大暴雨事件，本次暴雨累计降雨量大、持续时间长、影响范围广。以北京为例，7月29日20时至8月2日7时，北京市平均降雨量达到331 mm，这83 h内降雨量已达到常年全年降雨量的60%，受灾严重的门头沟区平均降雨量为538.1 mm，房山区平均降雨量为598.7 mm。本次强降水导致北京房山、门头沟和河北涿州等地出现内涝，北京、河北多地出现山洪、塌方等灾害。此次暴雨是在台风“杜苏芮”残余环流与副热带高压、台风“卡努”及地形综合作用下形成的，台风“杜苏芮”北上时携带的大量水汽显著影响此次降雨的发生规模。因此，作为水汽含量重要的表现形式[1]，大气可降水量（precipitable water vapor，PWV）可用于分析此次降水的发生过程。

常见的水汽探测的方法有探空、卫星遥感等，但都难以获取降雨期间准确连续的水汽信息。探空往往受限于站点的数量和发射频率，而卫星遥感则受到轨道重返周期及云层的影响[2]。而基于GNSS观测进行大气水汽反演，具有低成本、高精度、高时空分辨的优势[3]。近年来，相关学者尝试开展了GNSS水汽与降水在时间及空间尺度的相关研究。文献[4]基于印度尼西亚的观测验证了PWV与降雨量的正相关关系；文献[5]开展了我国中南半岛地区2006—2016年GNSS水汽产品及降水数据分析，并发现PWV和降水的相关性呈现随纬度减小而降低的特点；文献[6]探究了GNSS-PWV的时空特征与梅雨过程的关系，并发现PWV在降雨过程中一直处于高值状态、降雨前后会出现PWV的迅速上升和逐渐下降过程；文献[7]利用GNSS观测开展了河南“21·7”极端暴雨过程的PWV与降雨过程研究。此外，基于降雨发生前会伴随PWV累积的研究结论，相关学者尝试建立了基于GNSS-PWV的降雨预测模型。文献[8]利用PWV上升期的增量及增率建立降水预报模型，可提前2~6 h预报出80%~90%的降水事件和90%以上的暴雨事件；文献[9]提出一种非线性自回归神经网络的48 h降雨预报模型，实测降雨时间序列与模型预测降雨时间序列的相关系数可达0.90；文献[10]基于改进的BPNN法，利用PWV及多个气象参数建立的降雨预报模型，其正确率可达96%；文献[11]基于NARX神经网络，提出融合GNSS-PWV与气象参数的短时强降水预报方法，提高了对强降水事件的检测精度；文献[12]利用PWV绝对值、变化量及其变化率建立了短期降水预报模型，其预测正确率可以达到80%左右；文献[13]基于机器学习方法建立了一种线性-非线性相结合的短期降水预报方法，可有效降低降水预报的误报率。

本文针对此次北京及周边地区的极端暴雨，利用该地区GNSS观测数据反演得到各站点PWV，基于一种改进的插值算法获取研究区0.25°×0.25°的PWV格网数据，并采用探空、气象再分析资料进行GNSS-PWV精度的多角度评估；进一步，结合气象站提供的暴雨发生过程的降雨信息，分析了降雨量与GNSS-PWV、梯度在时间和空间上的关系，可为建立极端降雨灾害预报系统提供科学的理论基础和数据支撑，进一步减少极端天气带来的社会损失。

北京及周边地区属于典型的温带半湿润半干旱季风气候，冬季寒冷而干燥，夏季炎热而潮湿，降雨受季节影响明显并集中在7月下旬。该地区位于华北平原，北靠燕山山脉，西部紧邻太行山，东邻渤海，西北和北面地形较高，南面和东面地形较为平坦，呈现出西北高东南低的地形特点，如图1（a）高程信息所示。7月29日，台风“杜苏芮”残余环流北上至京津冀地区，携带的大量水汽在此聚集；同时西太平洋的副热带高压造成的强烈气压梯度显著增强了该区域东南风力；在此东南风的作用下，台风“卡努”携带的大量水汽也被输送至华北平原，大量的水汽在此区域聚集，并在太行山东侧抬升导致降雨。因此，位于太行山东侧的石家庄、保定和北京西南部的房山、门头沟成为本次降水过程的核心影响区域，如图1（b）所示，降雨量高达600 mm。

图1 北京及周边地区的GNSS站点、探空站及气象站分布和总降水量

Fig. 1 Distribution of GNSS stations, radiosonde stations and meteorological stations in Beijing and its surrounding areas and total precipitation

本文使用的GNSS站分布如图1（a）中蓝色圆圈所示，14个测站分布在北京市内并集中于北京市中南部，11个测站位于河北省保定市，7个站点分布在河北省石家庄市。上述GNSS数据来自北京迅腾智慧科技股份有限公司的CORS，时间范围为2023年7月25日至2023年8月1日，反演了1 h分辨率的GNSS-PWV。本文选择与GNSS测站共址的地面气象站用于降水分析，其位置分布如图1（a）中红色三角所示。上述降水资料来自中国气象局的逐小时累计降雨量观测，时段为2023年7月25日至8月1日。

此外，本文还选用了研究区域涉及的探空数据和ERA5再分析数据。作为获取高精度PWV的测量方式，探空站通过发射搭载无线电探空仪的探空气球，来测量不同高度处的气象数据，以获取有关大气垂直分布的信息[14]。在本文中，使用的探空站编号为54 511，具体位置如图1（a），其数据来自美国怀俄明大学网站（

对地面上空每个气压层的可降雨量进行数值积分，探空和ERA5的PWV计算过程[17]如下

(1)

式中，g表示当地的重力加速度，通常取值为9.8 m/s2；PT表示对流层顶气压；PG表示地面层气压，单位为hPa；q表示比湿，单位为g/kg，其计算过程如下

(2)

式中，e表示水汽压，其计算公式如下

(3)

式中，a、b的数值由温度决定，当温度低于40℃时，a=21.87，b=265.49，当温度高于0℃时，a=17.26，b=237.79，当温度介于0℃与40℃，则通过内插计算其数值[18]。

引入构成长基线的IGS测站（如JFNG、LHAZ、URUM），采用GAMIT10.71进行双差解算，能够获取天顶总延迟（zenith total delay，ZTD）及东西向和南北向的对流层梯度因子。其中，ZTD由天顶静力学延迟（zenith hydrostatic delay，ZHD）和天顶湿延迟（zenith wet delay，ZWD）组成[19]，如式（4）所示

(4)

式中，ZHD可基于Saastamoinen模型准确估计[20]，计算过程如下

(5)

式中，P0是指地表气压，单位为hPa；H和φ分别表示测站处的大地高和纬度，单位分别为弧度和km。将上述模型计算的ZHD从其对应的ZTD中减去，便可对剩余的ZWD进行PWV转换[21]，计算过程如下

(6)

(7)

式中，Π和Tm分别代表转换因子和测站上空的加权平均温度；Rw表示水汽气体常数，其值为461.495 J×kg-1 k-1；ρw代表液态水的密度；k'2和k3指折射率常数，其值分别为（22.13±2.20）K/hPa、（3.739±0.012）×105 K2/hPa。

为更加准确地了解研究区PWV的空间分布变化，基于“移去-内插-恢复”思想[22]，提出一种改进的插值算法来获取研究区内0.25°×0.25°的PWV空间格网数据。其具体步骤如下。

（1）需要将GNSS-PWV归算至统一的高程面，具体公式如下[23]

(8)

式中，PWV0和H0分别表示GNSS-PWV和测站处的高程；Hh和PWVh表示归算的目标高度和目标高度处的PWV。

（2）根据统一高度处GNSS各测站经纬度坐标和时间信息，利用GPT3模型分别计算出各测站处和各0.25°格网时处的水汽压、水汽压递减率等气象参数，并基于式（9）计算ZWD估值[24]

(9)

（3）利用式（6）和式（7）即可得到PWV，记为GPT3-PWV。在各测站计算GPT3-PWV与GNSS-PWV的差值，记为PWV-残差，并对其进行克里金内插，以得到研究区内统一高度处0.25°的PWV-残差格网分布，进而可恢复各格网PWV信息。

（4）再次利用式（8），结合美国国家航天航空局NASA发布的DEM产品，将每个格网点的PWV数值归算至其对应的真实地表。至此，可以得到研究范围内小时分辨率的GNSS空间格网数据。

将探空计算得到的PWV作为参考值，对GNSS-PWV精度进行评估。已有研究结果表明，当探空站与GNSS站水平距离小于50 km时，可认为二者共址[25]。而本文研究中有9个GNSS测站满足上述要求，与54 511探空站共址。为此，本文使用两种方法进行PWV精度对比：①匹配方法1，对9个共址的GNSS站点PWV进行高程改正，使其位于统一的高程面，然后使用二维的克里金插值法，得到探空站位置处的PWV，记为GNSS-PWV1；②匹配方法2，直接选择与探空站最近的GNSS站，进行高程改正得到相应PWV，并记作GNSS-PWV2。

图2展示了在台风北上导致的北京等地暴雨过程中，探空PWV与两种方法获得的GNSS-PWV的时序分布。由图2可以看出，两种方法得到的PWV具有高度的一致性，且都与探空PWV呈现相同的变化趋势：即在暴雨前累积上升，暴雨结束后呈下降状态；探空站和GNSS的PWV峰值皆出现在7月31日12时，三者最大差值为6 mm出现在7月26日0时，最小差异为0.6 mm出现7月28日10时。

图2 RS-PWV与GNSS-PWV时序

Fig. 2 Sequence of RS-PWV and GNSS-PWV

图3展示了两种GNSS-PWV与RS-PWV的相关性。图3（a）中，RS-PWV与GNSS-PWV的相关系数（r）、bias和RMSE分别为0.99、-0.52 mm、0.52 mm；图3（b）中相关指标分别为0.99、-0.40 mm、0.40 mm。可以看到，将探空计算的PWV作为参考值，所选取的两种比较都表明GNSS反演的PWV具有较高的精度。

图3 RS-PWV与GNSS-PWV的相关性

Fig. 3 Correlation between RS-PWV and GNSS-PWV

本文利用ERA5计算了各个GNSS站点的PWV，并与GNSS反演的PWV进行了比较。图4展示了34个GNSS站点的bias和RMSE分布。可以看到，所有测站PWV的bias都小于1.5 mm，而RMSE在2.5~5.5 mm之间。受暴雨影响，北京市南部、保定市东北部、石家庄市中部bias、RMSE值稍高，但此精度可满足后续对GNSS-PWV作时空特征分析[26]。

图4 GNSS-PWV与ERA5-PWV的bias和RMSE分布

Fig. 4 Distribution of bias and RMSE in GNSS-PWV and ERA5-PWV

此外，基于改进内插算法得到的GNSS格网PWV也同ERA5格网PWV进行了比较，图5展示了其bias和RMSE分布。可以看到，研究区bias处于-2.5~4 mm，尤其在GNSS测站密集区域其bias小于1 mm；研究区RMSE处于3~7 mm，且超过90%区域小于4 mm。研究区bias和RMSE的均值分别为-0.75 mm和3.5 mm，反映出GNSS反演PWV的高精度。

图5 GNSS-PWV空间格网数据与ERA5-PWV的bias和RMSE

Fig. 5 GNSS-PWV spatial grid data and ERA5-PWV bias and RMSE

本文选取了受灾严重区域的4个GNSS测站及其共址气象站进行极端暴雨过程PWV与降雨的时间变化分析。其中，GNSS站点分别为石家庄市sjps、保定市bddx、北京市sjgb和gyyc。图6展示了4个测站7月25日0时至8月1日23时期间PWV变化与小时降雨量的时序关系，虚线描绘出了集中降雨的时间。可以看到，4个站点的PWV均呈现先增后减的趋势，即PWV在暴雨发生前及发生期间处于波动上升状态，各测站降雨量达到最大时，对应时刻的PWV皆高于60 mm；在强降水结束后，水汽仍维持在较高水平，无法被立刻释放至降水前的数值大小。

图6 GNSS-PWV和降雨量的时序关系

Fig. 6 Temporal relationship between GNSS-PWV and rainfall

其中，sjps站在29日21时降雨量达到最大40 mm，其他持续降雨发生时的雨量约在20 mm。在29日至31日集中降雨期间，PWV在最大降雨发生后仍一直处于较高水平，表明水汽仍在积累；31日后PWV有小幅度降低，但仍维持在55 mm的高位，预示降水还将发生。本次降水持续时间长，从28日16时开始降水至29日凌晨3时降雨强度增大时，PWV增加量达到了15 mm，平均每小时增长1.5 mm；相应在29日0时至2时PWV变化率最大，达到4 mm/h。

bddx站在30日7时降雨量达到最大43 mm，29日4时降雨量为24 mm，其余时间降雨量较少处于10 mm以下。从28日11时至29日发生降水前，PWV持续上升并累积了10 mm；29日0时至3时，PWV在短时间内急剧上升，变化率达到2 mm/h；在29日4时至30日之间，降雨量虽保持较低水平，但水汽仍持续积累甚至高达70 mm；这为30日开始的高强度暴雨提供了丰富的水汽积累，反映了水汽含量和降水的前置关系。

sjgb站存在两次强降雨，分别出现在30日10时和31日3时，降雨量为42 mm和51 mm；在26日9时出现低强度降雨时，PWV提前8 h从30 mm增加至49 mm，降雨结束后，PWV逐渐降至40 mm；29日7时开始，降雨量逐渐增加，PWV同时从43 mm增加至60 mm以上；31日3时后，降雨量逐渐降至10 mm及以下，PWV也随之降至约50 mm。

gyyc站暴雨开始时间稍晚于其他站点，28日8时，PWV开始呈现累积上升状，29日21时降雨量达到第1个峰值23.9 mm，期间PWV增加了20 mm；30日13时至21时降雨量降至5 mm以下，而后到31日2时降雨量迅速增加达到最大值50.9 mm，这段时间内，PWV持续上升至76 mm。强暴雨结束后，PWV没有立刻消散；从31日12时开始，PWV逐渐降低，到20时，PWV从78 mm降至68 mm；而后少量降水期间，PWV虽有小幅度上升，但整体趋势保持下降。

为更直观地分析暴雨前PWV的变化量，分别计算了各站点暴雨前（7月25日至29日）和暴雨期间（7月29日至8月1日）PWV的均值，以及二者的变化率，结果如图7所示，其中，蓝色代表暴雨发生前PWV均值，红色代表暴雨发生期间PWV均值，红色折线代表二者的变化率。降雨前，各站点PWV均值约为50 mm，降雨期间，PWV均值约为70 mm，且各站点变化率皆高于30%，尤其在北京市的bjyq站点，增幅达到49%，表明暴雨的发生伴随PWV大量聚集的现象。

图7 暴雨发生前及发生期间GNSS-PWV均值变化及其变化率

Fig. 7 Variation and rate of GNSS-PWV mean value before and during rainstorm

为进一步探究暴雨过程PWV与降雨量的空间变化关系，选取降雨期间4个时刻的降雨量、PWV及其2 h前PWV进行分析并绘制空间分布如图8。可以看到，7月29日8时，北京地区的PWV为55~60 mm且处于整体上升期；石家庄及保定南部地区PWV相对较高，达到65 mm左右，该区域有降水产生，降雨量都小于10 mm。7月30日9时，PWV整体升高，尤其是南部地区普遍达到了65 mm以上；北京地区PWV也有大幅度上升，2 h后整个范围内有降水产生且降雨量在10~20 mm。31日10时，北京及其周边的廊坊地区，PWV处于最高值，普遍维持在70 mm左右，而石家庄地区的PWV相比前两日有所降低。由图8（a）—（c）可以明显看出，在此次强暴雨中，水汽呈现由西南向东北逐渐转移的趋势；图8（h）表现出石家庄及保定地区降雨强度逐渐减弱，强降雨主要集中在北京，这种情况和水汽由西南向东北的输送趋势相符。至8月1日7时，强降雨主要集中在北京东南部，对应PWV处于较高值，约在60 mm，其余地区降水低于10 mm。由图8（g）—（l）可以看出，虽然部分地区降雨已经结束，水汽含量有所较低，但对比降雨初期仍处于较高水平，在55~65 mm之间，并未立刻释放。

图8 降雨量GNSS-PWV空间变化

Fig. 8 GNSS-PWV spatial variation of rainfall

采用空间插值思想生成的北京及其周边地区水汽分布图，能更加直观地显示出大范围的水汽分布情况，并能体现出水汽的输送方向。从空间上来看，有强降雨发生地区的PWV比周边区域高出10~20 mm，降雨结束后PWV有所降低；PWV较高的区域与有降水的区域高度重合，证明了两者间的相关性。

图8反映出此次强降雨中心由河北西南部向北京北部转移，GNSS-PWV高值也呈现相同的趋势。为更清晰地判断对流层水汽与降雨走势的关系，本文进一步绘制了7月29日至8月1日每隔6 h的对流层延迟梯度变化图，如图9所示。可以看出，在每个时刻，几乎所有GNSS站点处的对流层梯度都指向东北方向，表明东北方向比西南方向有相应更大的水汽含量；7月29日，位于河北石家庄市GNSS站点对流层延迟梯度大多指向北方向，部分指向西南方向，至8月1日，各站点对流层延迟梯度的方向逐渐趋于一致，指向北京的西南部；由图8（k）可以看出，此时北京西南部降雨量仍维持在30 mm，其余站点降雨量皆在10 mm以下，这表明对流层延迟梯度的指向和降雨中心的变化具有一致性。

图9 对流层延迟梯度变化

Fig. 9 Variation of the tropospheric delay gradient

针对2023年北京及周边区域的特大暴雨，本文利用GNSS观测反演了UTC 2023年7月25日至8月1日期间逐小时的GNSS-PWV，并通过一种改进的插值算法生成了研究区内0.25°的GNSS-PWV格网数据；结合气象站降雨量数据，从时间与空间的角度探究了GNSS-PWV、梯度变化与此次极端暴雨的关系，结果表明：

（1）GNSS-PWV与RS-PWV相关性高达0.99，RMSE为0.52 mm；GNSS-PWV与ERA5-PWV在各站点处的RMSE范围为2.5~5.5 mm；对于GNSS格网数据，超过90%地区RMSE在5 mm以下，这些结果皆表明GNSS-PWV精度满足降雨分析的需求。

（2）GNSS测站PWV与共址气象站实测降雨量的时间序列分析表明，本次暴雨前后PWV整体呈现先增后减的趋势；降雨发生前12~14 h内，PWV开始迅速增加，且在暴雨发生期间，PWV仍呈现累积上升的状态，表明水汽在持续积累并导致此次暴雨持续时间长；从暴雨发生前夕至降雨量达到顶峰，PWV相差约20 mm。此外，若在降雨产生后，PWV并未降低或只减少了较小的值，代表大气中水汽释放不完全，还将出现降雨现象。

（3）PWV格网数据与降雨量空间分布分析表明，水汽高值与降雨的空间变化趋势一致，都为“西南至东北”；且强降雨发生区域，PWV数值也较大，降雨量减少后，PWV数值也会相应降低。

（4）此次暴雨过程中，GNSS站点的对流层延迟梯度指向东北方向，表明东北方向比西南方向有更大的水汽含量；且随时间变化，这种方向性更加明显，表明对流层延迟梯度的指向性与强降雨中心的转移具有一致性。

来源：测绘学报