CN103678885A - 一种基于重力卫星的干旱指数构建及分析方法 - Google Patents

Abstract

本发明公开了一种基于重力卫星的干旱指数构建及分析方法。该方法包括：步骤一、重力卫星观测数据预处理；步骤二、区域平均陆地总水量变化过程提取；步骤三、区域陆地总水量空间分布计算；步骤四、区域陆地总水量标准化处理；步骤五、计算并分析基于重力卫星的干旱指数总水量异常指数时间序列及空间分布。本发明构建的干旱指数能有效反映陆地总水量的季节性和年际变化规律，科学揭示了干旱的成灾过程；考虑了干旱历时的致灾机制，能有效反映前期严重干旱对后期灾情的影响机理，突破了资料匮乏地区或无资料地区干旱监测与评估中存在的科学瓶颈问题。

Description

一种基于重力卫星的干旱指数构建及分析方法

技术领域

本发明涉及一种基于重力卫星的干旱指数构建与分析方法，尤其涉及一种基于GRCAE重力卫星数据反演的陆地总水量构建能够反映干旱历时影响的干旱指数的方法，属于重力场观测与水文学研究相结合的交叉学科技术领域。

背景技术

干旱是全世界范围内影响社会经济可持续发展的重大自然灾害之一，具有持续时间长、影响范围广、影响领域多等特点。干旱定义及强度分类没有统一的标准，不同干旱指标由于其构建目的、所依赖的水文输入数据及计算方法不同而各有其优缺点及适用性。

例如：

（1）、标准降雨指数（SPI）适用于研究不同时间尺度的干旱，计算稳定，对干旱反应较灵敏，但由于没有考虑蒸散发水分支出，故不能反映异常高温对干旱的加剧作用。

（2）、帕尔默干旱指数(PDSI)同时考虑了水分供给(降水)和水分支出(潜在蒸散发)，可以反映干旱的程度以及干旱历时。

（3）、水文干旱指数（HDI）利用地表径流和水库、河道水位等来构建干旱指标，由于缺乏流域土壤含水量、地下水埋深等实测数据，并不能很好地代表区域的整体情况。

因此，在实际应用中，上述三种指数都严重受到资料匮乏的限制和影响，其科学性、正确性、实时性大打折扣。

由于人财物能力限制，全世界很多地区土壤含水量、地下水埋深等实测水文资料都非常匮乏，严重制约了水资源的科学与合理利用，影响了社会经济的可持续发展。美国宇航局(NASA)和德国空间飞行中心(DLR)联合开发的GRACE重力卫星计划有效解决了大范围、长时段陆地水资源储量变化监测的科学难题，为研究气候变化条件下全球或区域陆地水资源储量变化提供了有力的数据支撑。

虽然，GRACE陆地水资源储量已被科学家应用于干旱监测评估研究。然而，目前基于GRACE陆地水资源储量的干旱指数国内外未见报道。因此，缺乏一个科学合理的量化指标来实现水资源和干旱的大范围、长时段实时监测和分析比较。

发明内容

本发明所要解决的技术问题是利用GRACE重力卫星数据，提出一种干旱指数构建和分析方法，能够从GRACE重力卫星观测的地球重力场模型中反演研究区陆地总水量变化过程，能快速实现广大资料匮乏或无资料地区的干旱监测，拓展重力场观测技术在陆地水资源及水旱灾害事件研究中的应用空间，填补无资料地区干旱监测的时空盲区，显著提高资料匮乏地区的干旱监测科学水平。

本发明为解决上述技术问题采用以下技术方案：

一种基于重力卫星的干旱指数构建及分析方法，具体包括如下步骤：

步骤1，对GRACE重力卫星观测数据进行预处理；

步骤2，从各月重力场模型中提取研究区平均陆地总水量变化过程；

步骤3，计算研究区陆地总水量空间分布；

步骤4，对陆地总水量进行归一化处理消除总水量季节性循环及地区差异的影响，具体为：

对研究时段内各月总水量扣除该月总水量同期平均值，得到该月总水量异常值，计算各月总水量异常对该月总水量变幅的比值，得到总水量异常百分比PTSA，计算公式如下：

${\mathrm{PTSA}}_{i,j}=\frac{{\mathrm{TSA}}_{i,j}-{\mathrm{MTSA}}_j}{{\mathrm{MTSA}}_j-\min {\mathrm{TSA}}_j}\×100,{\mathrm{TSA}}_{i,j}\≤{\mathrm{MTSA}}_j$

${\mathrm{PTSA}}_{i,j}=\frac{{\mathrm{TSA}}_{i,j}-{\mathrm{MTSA}}_j}{\max {\mathrm{TSA}}_j-M{\mathrm{TSA}}_j}\×100,{\mathrm{TSA}}_{i,j}>{\mathrm{MTSA}}_j$

式中：TSA_i,j为第i年第j月总水量异常，MTSA_j为各年第j月总水量平均值，minTSA_j和maxTSA_j分别为研究时段内第j月总水量最小值和最大值；

步骤5，构建考虑干旱历时的总水量异常指数，具体为：

首先，将极端干旱情况下总水量异常百分比累积量用一条直线表示：

$\underset{t=1}{\overset{j}{\mathrm{\Σ}}}{\mathrm{PTSA}}_t=-100j-100$

式中t为干旱的月数，j为参与计算的总月数；

其次，以此直线代表极端干旱的情况，将各月总水量异常指数TSAI表示为：

${\mathrm{TSAI}}_j=\frac{{\mathrm{\Σ}}_{t=1}^j{\mathrm{PTSA}}_t}{25j+25}$

之后，采用月增量形式计算干旱指数以考虑历时的影响，将各月总水量异常指数表示为：

${\mathrm{TSAI}}_j={\mathrm{TSAI}}_{j-1}+\frac{{\mathrm{PTSA}}_j}{50}-0.5{\mathrm{TSAI}}_{j-1}$

TSAI在-4到4之间变动，分别代表极端干旱至极端湿润的情况。

作为本发明的基于重力卫星的干旱指数构建方法进一步的优化方案，所述步骤1中对GRACE重力卫星观测数据进行预处理依次包括：

（1）数据的低阶项校正：使用卫星激光测距观测得到的C₂₀项对GRACE重力位模型的C₂₀项进行替换；

（2）去条带滤波：使用5阶多项式，对各月重力场模型8次以上各阶系数进行多项式拟合去除高阶项相关误差；

（3）高斯平滑滤波：使用平滑半径为200km的各项同性高斯平滑滤波器，去除GRACE重力场模型高阶项误差的影响。

作为本发明的基于重力卫星的干旱指数构建方法进一步的优化方案，所述步骤2是采用核函数法从各月重力场模型中提取研究区平均陆地总水量变化过程，其中基于核函数的研究区平均陆地总水量变化过程提取公式为：

$\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\σ}}}_{\mathrm{region}}=\frac{{\mathrm{a\ρ}}_E}{3{\mathrm{\Ω}}_{\mathrm{region}}}\underset{l=0}{\overset{l_{\max }}{\mathrm{\Σ}}}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}\frac{2l+1}{1+k_l}(W_{\mathrm{lm}}^c\mathrm{\Δ}{C}_{\mathrm{lm}}+W_{\mathrm{lm}}^s\mathrm{\Δ}{S}_{\mathrm{lm}})$

其中，为研究区平均陆地总水量变化，a为地球半径，ρ_E为地球平均密度，l和m分别为阶数和次数，l_max为重力场模型最高阶数，k_l为l阶荷勒夫数，表示地壳因表面荷载变化产生弹性形变进而引起的地球重力场变化，

和

分别为核函数球谐展开系数，ΔC_lm和ΔS_lm为月重力位系数异常值，Ω_region为区域面积。

作为本发明的基于重力卫星的干旱指数构建方法进一步的优化方案，步骤3所述计算研究区陆地总水量空间分布，其计算公式为：

式中Δh(θ,φ)为等效水高，a为地球半径，

和θ分别为计算点经度和余纬，

为归一化缔合勒让德多项式，ΔC_lm和ΔS_lm是月重力位系数异常值，l和m分别为阶数和次数，ρ_E和ρ_wat分别为地球平均密度和水的密度，l_max为重力场模型最高阶数，k_l为l阶荷勒夫数，表示地壳因表面荷载变化产生弹性形变进而引起的地球重力场变化。

本发明采用以上技术方案与现有技术相比，具有以下技术效果：

本发明运用GRACE观测的陆地总水量构建总水量异常指数(Total Storage Anomaly Index,TSAI)。基于总水量的干旱指数能够反映降雨、蒸散发等水文-气象过程对陆地水资源储量的综合作用，且克服了实测土壤含水量、地下水埋深资料稀缺和难以获取等科学瓶颈，充分拓展重力场观测技术在全球陆面水文过程及极端灾害事件研究中的应用空间，填补了无资料地区干旱监测的时空盲区，显著提高了资料匮乏地区的干旱监测科学水平。

本发明对输入数据要求少，计算简便快捷，有效解决了数据采集的困难，显著节约了数据采集的人、财、物成本，而且计算结果与目前全世界广泛使用的帕尔默干旱指数PDSI有较高的一致性。本发明适合于大尺度干旱灾害监测，突破了资料匮乏或无资料地区旱情监测评估的科技瓶颈，对抗旱减灾的生产实践具有较高推广应用价值。

附图说明

图1为本发明一种基于重力卫星的干旱指数构建及分析方法的流程图。

图2为GRCAE数据预处理效果图(单位：mm)。

图3为基于核函数提取研究区平均陆地总水量变化过程图。

图4为代表月份研究区陆地总水量空间分布图(单位：mm)。

图5A为研究区总水量异常百分比变化过程图。

图5B为代表月份研究区总水量异常百分比空间分布图。

图6A为研究区总水量异常指数、帕尔默干旱指数（PDSI）、30天标准降雨指数变化过程对比图；

图6B为代表月份总水量异常指数空间分布图。

具体实施方式

下面结合附图对本发明的技术方案做进一步的详细说明：

本发明提出一种基于重力卫星的干旱指数构建及分析方法，具体包括如下步骤：

步骤1，对GRACE重力卫星观测数据进行预处理；

步骤2，从各月重力场模型中提取研究区平均陆地总水量变化过程；

步骤3，计算研究区陆地总水量空间分布；

步骤4，对陆地总水量进行归一化处理消除总水量季节性循环及地区差异的影响；

步骤5，构建考虑干旱历时的总水量异常指数，分析从极端干旱至极端湿润的情况。

GRACE(Gravitv Recovery and Climate Experiment，重力恢复与气候实验)重力卫星计划由美国国家宇航局(NASA)和德国空间飞行中心(DLR)联合开发，该卫星于2002年3月17日在俄罗斯北部的普列谢茨克(Plesetsk)发射中心发射成功。

GRACE卫星采用低低卫—卫跟踪技术，即同时发射相距约220km的两颗在同一轨道上的低轨卫星(卫星几何尺寸：长3123mm，宽1942mm，高720mm)，两颗低轨卫星除去有星载GPS接收机准确确定其轨道位置外，还以微米级精度实时测量两个低轨卫星之间的距离及其变化率，以三轴加速度计测量非保守力，用S波段无线电波将每颗卫星的观测数据，包括重力相关测量和GPS掩星测量，传向地面监控站，表1给出了搭载仪器设备名称及其功能。

表1GRACE重力卫星主要载荷参数

仪器设备	功能
		K波段测距系统	测量两颗卫星之间的距离
S波段天线	从卫星至地面传输数据
		SuperSTAR加速度仪	精密测量非引力加速度
恒星照相仪	精密测量卫星相对于恒星的姿态
		超稳定振荡器	提供K波段测距系统所需频率
地球、太阳传感器	跟踪地球和太阳，精度较低，用于调节姿态
		质心调节仪	精密测量并调整质心的位置
GPS接收机	用于精密定轨和大气掩星研究
		激光反射镜	测量卫星至地面跟踪站距离
3-轴姿态稳定控制系统	稳定卫星姿态

对照图1，本发明基于2003年1月至2012年12月GRACE重力卫星数据，以揭示云南省2009～2010秋冬春连旱为例，阐明本发明的实施方式：

对照图2，如图2中（a）所示，未经预处理的GRACE重力卫星数据反演的陆地总水量信号被南-北方向条带噪声和高阶项噪声覆盖，很难识别有效的陆地总水量变化信号。采用卫星激光测距观测得到的C₂₀项对GRACE重力位模型的C₂₀项进行了替换；使用5阶多项式，对各月重力场模型8次以上各阶系数进行多项式拟合去除高阶项相关误差；并使用平滑半径为200km的各项同性高斯平滑滤波器，去除GRACE重力场模型高阶项误差的影响。如图2中（b）所示，经过预处理的GRACE重力场模型反演的陆地总水量噪声大幅降低，水文信号能够得到清晰反映。

根据研究区边界构建研究区核函数，基于核函数的区域平均陆地总水量提取公式为：

$\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\σ}}}_{\mathrm{region}}=\frac{{\mathrm{a\ρ}}_E}{3{\mathrm{\Ω}}_{\mathrm{region}}}\underset{l=0}{\overset{l_{\max }}{\mathrm{\Σ}}}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}\frac{2l+1}{1+k_l}(W_{\mathrm{lm}}^c\mathrm{\Δ}{C}_{\mathrm{lm}}+W_{\mathrm{lm}}^s\mathrm{\Δ}{S}_{\mathrm{lm}})$

其中，为区域或区域平均陆地总水量变化，a为地球半径，ρ_E为地球平均密度，l和m分别为阶数和次数，l_max为重力场模型最高阶数，k_l为l阶荷勒夫数，表示地壳因表面荷载变化产生弹性形变进而引起的地球重力场变化，和

为核函数球谐展开系数，ΔC_lm和ΔS_lm为月重力位系数异常值，Ω_region为区域面积。

对照图3，GRACE观测的云南地区总水量季节性变化规律明显，总水量最高值大致出现在7-9月份，最低值出现在2-4月份。总水量年际间变化特征显著，其中2009-2010年特大干旱在总水量过程线上有明显反映。2008年总水量为研究时段内最高，此后三年陆地总水量逐年下降，其中2010年春季总水量为整个研究时段内最低，表明GRACE观测的总水量可以有效反映西南特大干旱的形成过程。

对照图4，本例只给出代表月份研究区陆地总水量空间分布，计算区域陆地总水量空间分布，计算公式为：

式中Δh(θ,φ)为等效水高，a为地球半径，φ和θ分别为计算点经度和余纬，

为归一化缔合勒让德多项式，ΔC_lm和ΔS_lm就是月重力位系数异常值，l和m分别为阶数和次数，ρ_E和ρ_wat分别为地球平均密度和水的密度，k_l为l阶荷勒夫数，表示地壳因表面荷载变化产生弹性形变进而引起的地球重力场变化。

对区域陆地总水量标准化处理。对各月总水量扣除该月总水量同期平均值，得到该月总水量异常值，计算各月总水量异常对该月总水量变幅的比值，得总水量异常百分比(PercentageTotal Storage Anomaly，PTSA)，结果如图5A所示，计算公式如下：

${\mathrm{PTSA}}_{i,j}=\frac{{\mathrm{TSA}}_{i,j}-{\mathrm{MTSA}}_j}{{\mathrm{MTSA}}_j-\min {\mathrm{TSA}}_j}\×100,{\mathrm{TSA}}_{i,j}\≤{\mathrm{MTSA}}_j$

${\mathrm{PTSA}}_{i,j}=\frac{{\mathrm{TSA}}_{i,j}-{\mathrm{MTSA}}_j}{\max {\mathrm{TSA}}_j-M{\mathrm{TSA}}_j}\×100,{\mathrm{TSA}}_{i,j}>{\mathrm{MTSA}}_j$

式中：TSA_i,j为第i年第j月总水量异常，MTSA_j为各年第j月总水量平均值，minTSA_j和maxTSA_j分别为研究时段内第j月总水量最小和最大值。如图5A和5B所示，本例给出研究区总水量异常百分比变化过程及代表月份空间分布。

对照图6A，计算基于重力卫星的干旱指数总水量异常指数时间序列及空间分布。极端干旱情况下总水量异常百分比累积量可以用一条直线表示：

$\underset{t=1}{\overset{j}{\mathrm{\Σ}}}{\mathrm{PTSA}}_t=-100j-100;$

式中t为干旱的月数，j为参与计算的总月数。以此直线代表极端干旱(干旱等级为-4)的情况，总水量异常指数（Total Storage Anomaly Index，TSAI）可以表示为：

${\mathrm{TSAI}}_j=\frac{{\mathrm{\Σ}}_{t=1}^j{\mathrm{PTSA}}_t}{25j+25};$

采用月增量形式计算干旱指数以考虑历时的影响，各月总水量异常指数可以表示为：

${\mathrm{TSAI}}_j={\mathrm{TSAi}}_{j-1}+\frac{{\mathrm{PTSA}}_j}{50}-0.5{\mathrm{TSAI}}_{j-1};$

TSAI在-4到4之间变动，分别代表极端干旱和极端湿润的情况。

对照图3和图6A，虽然总水量在2010年2月达到最低值以后即开始恢复，但由于前期严重干旱的影响，干旱灾害并未得到缓解且一直持续到2010年5月达到极端干旱，表明考虑干旱历时影响的TSAI可以反映前期严重干旱对后期干旱程度的影响。

对照图4（a）、（b）和图6B（a）、（b），虽然2009年7月总水量较2009年4月丰沛，但TSAI显示的2009年7月干旱程度比2009年4月严重，表明基于GRCAE观测的总水量构建的干旱指数有效消除了总水量季节性循环变化对干旱评估的影响，相比于直接运用总水量评估干旱灾害更加直观。

对照图6A，本发明构建的TSAI指数与PDSI指数一致性较好。然而，PDSI指数计算需要降雨、径流、蒸散发、气温、风速等水文气象资料以及土壤蓄水容量等土壤特性资料，对资料要求高。而TSAI直接利用GRACE观测的陆地水储量构建，资料获取方便，计算便捷。因此，TSAI更适合于资料匮乏地区的干旱监测与评估。

对照图6A，SPI30相比于TSAI波动较为明显，这是因为SPI30只考虑降雨单因子导致的气象干旱灾害，对降雨变化敏感。但由于因子单一，不能反映干旱过程中各因子综合作用对陆地水量变化的影响。

以上实施方式仅用于说明本发明，而并非对本发明的限制，有关技术领域的普通技术人员，在不脱离本发明的精神和范围的情况下，还可以做出各种变化和变型，因此所有等同的技术方案也属于本发明的范畴，本发明的专利保护范围应由权利要求限定。

Claims (4)

1.一种基于重力卫星的干旱指数构建及分析方法，其特征在于，具体包括如下步骤：

步骤1，对GRACE重力卫星观测数据进行预处理；

步骤2，从各月重力场模型中提取研究区平均陆地总水量变化过程；

步骤3，计算研究区陆地总水量空间分布；

步骤4，对陆地总水量进行归一化处理消除总水量季节性循环及地区差异的影响，具体为：

对研究时段内各月总水量扣除该月总水量同期平均值，得到该月总水量异常值，计算各月总水量异常对该月总水量变幅的比值，得到总水量异常百分比PTSA，计算公式如下：

${\mathrm{PTSA}}_{i,j}=\frac{{\mathrm{TSA}}_{i,j}-{\mathrm{MTSA}}_j}{{\mathrm{MTSA}}_j-\min {\mathrm{TSA}}_j}\×100,{\mathrm{TSA}}_{i,j}\≤{\mathrm{MTSA}}_j$

${\mathrm{PTSA}}_{i,j}=\frac{{\mathrm{TSA}}_{i,j}-{\mathrm{MTSA}}_j}{\max {\mathrm{TSA}}_j-M{\mathrm{TSA}}_j}\×100,{\mathrm{TSA}}_{i,j}>{\mathrm{MTSA}}_j$

式中：TSA_i,j为第i年第j月总水量异常，MTSA_j为各年第j月总水量平均值，minTSA_j和maxTSA_j分别为研究时段内第j月总水量最小值和最大值；

步骤5，构建考虑干旱历时的总水量异常指数，具体为：

首先，将极端干旱情况下总水量异常百分比累积量用一条直线表示：

$\underset{t=1}{\overset{j}{\mathrm{\Σ}}}{\mathrm{PTSA}}_t=-100j-100$

式中t为干旱的月数，j为参与计算的总月数；

其次，以此直线代表极端干旱的情况，将各月总水量异常指数TSAI表示为：

${\mathrm{TSAI}}_j=\frac{{\mathrm{\Σ}}_{t=1}^j{\mathrm{PTSA}}_t}{25j+25}$

之后，采用月增量形式计算干旱指数以考虑历时的影响，将各月总水量异常指数表示为：

${\mathrm{TSAI}}_j={\mathrm{TSAi}}_{j-1}+\frac{{\mathrm{PTSA}}_j}{50}-0.5{\mathrm{TSAI}}_{j-1}$

TSAI在-4到4之间变动，分别代表极端干旱至极端湿润的情况。

2.根据权利要求1所述的基于重力卫星的干旱指数构建及分析方法，其特征在于，所述步骤1中对GRACE重力卫星观测数据进行预处理依次包括：

（1）数据的低阶项校正：使用卫星激光测距观测得到的C₂₀项对GRACE重力位模型的C₂₀项进行替换；

（2）去条带滤波：使用5阶多项式，对各月重力场模型8次以上各阶系数进行多项式拟合去除高阶项相关误差；

（3）高斯平滑滤波：使用平滑半径为200km的各项同性高斯平滑滤波器，去除GRACE重力场模型高阶项误差的影响。

3.根据权利要求1所述的基于重力卫星的干旱指数构建及分析方法，其特征在于，所述步骤2是采用核函数法从各月重力场模型中提取研究区平均陆地总水量变化过程，其中基于核函数的研究区平均陆地总水量变化过程提取公式为：

$\mathrm{\Δ}{\stackrel{\‾}{\mathrm{\σ}}}_{\mathrm{region}}=\frac{{\mathrm{a\ρ}}_E}{3{\mathrm{\Ω}}_{\mathrm{region}}}\underset{l=0}{\overset{l_{\max }}{\mathrm{\Σ}}}\underset{m=0}{\overset{l}{\mathrm{\Σ}}}\frac{2l+1}{1+k_l}(W_{\mathrm{lm}}^c\mathrm{\Δ}{C}_{\mathrm{lm}}+W_{\mathrm{lm}}^s\mathrm{\Δ}{S}_{\mathrm{lm}})$

其中，

为研究区平均陆地总水量变化，a为地球半径，ρ_E为地球平均密度，l和m分别为阶数和次数，l_max为重力场模型最高阶数，k_l为l阶荷勒夫数，表示地壳因表面荷载变化产生弹性形变进而引起的地球重力场变化，

和分别为核函数球谐展开系数，ΔC_lm和ΔS_lm为月重力位系数异常值，Ω_region为区域面积。

4.根据权利要求1所述的基于重力卫星的干旱指数构建及分析方法，其特征在于，步骤3所述计算研究区陆地总水量空间分布，其计算公式为：

式中Δh(θ,φ)为等效水高，a为地球半径，

和θ分别为计算点经度和余纬，

为归一化缔合勒让德多项式，ΔC_lm和ΔS_lm是月重力位系数异常值，l和m分别为阶数和次数，ρ_E和ρ_wat分别为地球平均密度和水的密度，l_max为重力场模型最高阶数，k_l为l阶荷勒夫数，表示地壳因表面荷载变化产生弹性形变进而引起的地球重力场变化。

Images