本发明公开了一种基于卫星降水数据的降雨侵蚀力计算方法，所述方法包括以下步骤：步骤1、获取数据；步骤2、计算降雨侵蚀力；步骤3、计算校正系数；步骤4、结果校正。本发明所述方法在卫星降水数据的基础上，结合不同类型子气候带的校正系数，对降雨侵蚀力数据进行校正，减小了由卫星降水数据精度问题所导致的降雨侵蚀力计算偏差；提高了降雨侵蚀力计算结果的准确性。
1.一种基于卫星降水数据的降雨侵蚀力计算方法，其特征在于，所述方法包括以下步骤：步骤1、获取数据：获取卫星降水数据、站点插值降雨侵蚀力数据以及气候带数据；
步骤2、计算降雨侵蚀力：根据获取的卫星降水数据，对单场降雨事件是否为侵蚀性降雨事件进行判定，若判定为侵蚀性降雨事件，则对侵蚀性降雨事件进行降雨侵蚀力计算，进而计算得到年平均降雨侵蚀力数据；
步骤3、计算校正系数：将步骤2中得到的年平均降雨侵蚀力数据与步骤1中获取的站点插值降雨侵蚀力数据在气候带数据上进行一元线性回归分析，通过计算得到校正系数，并对校正系数进行可靠性验证；
步骤4、结果校正：根据步骤3中得到的校正系数对步骤2中得到的年平均降雨侵蚀力数据进行校正，得到校正后的降雨侵蚀力数据；
其中，R表示降雨侵蚀力，单位为MJ mm h ha ；E表示降雨总动能，单位为MJ ha ；I30表‑1示30min内最大降雨强度，单位为mm h ；n表示侵蚀性降雨事件次数；
其中，k表示将单场侵蚀性降雨事件划分为k个片段；vr为第r个片段的降雨量；er为第r‑1 ‑1个片段的降雨动能，单位为MJ ha mm ；
根据RUSLE2中的降雨能量计算公式，则er计算公式为：‑1
其中，ir表示第r个片段降雨强度的大小，单位为mm h ；e为自然常数；
其中，Rm为年平均降雨侵蚀力，单位为MJ mm h ha yr ；m为年份，m≥2。
2.根据权利要求1所述的一种基于卫星降水数据的降雨侵蚀力计算方法，其特征在于，步骤1中所述气候带数据包括多个不同类型的子气候带数据，分别为：热带气候带数据；半干旱草原气候带数据；沙漠气候带数据；海洋性气候带数据；地中海气候带数据；亚热带湿润气候带数据；温带大陆性湿润气候带数据；温带季风气候带数据。
3.根据权利要求1所述的一种基于卫星降水数据的降雨侵蚀力计算方法，其特征在于，步骤2中所述对单场降雨事件是否为侵蚀性降雨事件进行判定的具体过程为：按下列条件进行判定：(1)单场降雨总量P总超过12.7mm；
(2)单场降雨任意30min内降雨量P30超过6.35mm；
(3)如果连续降雨在6h内降雨量P6h不超过1.27mm，则将该场连续降雨事件从第6h处划分为两场降雨事件，再分别利用条件(1)和(2)进行判断；
符合上述条件之一的降雨事件则判定为侵蚀性降雨事件。
4.根据权利要求2所述的一种基于卫星降水数据的降雨侵蚀力计算方法，其特征在于，步骤3中所述校正系数的具体计算过程为：读取步骤2中得到的年平均降雨侵蚀力数据以及步骤1中获取的站点插值降雨侵蚀力数据，根据不同类型的子气候带数据对计算得到的年平均降雨侵蚀力数据以及站点插值降雨侵蚀力数据分别进行裁剪；对同一子气候带类型的不同降雨侵蚀力数据进行一元线性回归分析，选取年平均降雨侵蚀力数据作为自变量x，站点插值降雨侵蚀力数据作为因变量y，进行回归分析，其一元线性回归方程为：y＝θ0+θ1x (5)
根据最小二乘法计算求得θ0与θ1，计算过程为：设有q组样本点：(xa，ya)，a＝1，…，q；
线性模型预测值 与真实值ya之间的差值,称为预测误差 其计算公式为：通过下述公式(8)来调整θ0与θ1的取值，使得预测误差 最小：最终得到该子气候带类型的校正系数θ1；
其中，SSR为回归平方和；SST为总平方和；yi为真实值； 为预测值； 为y平均值；
5.根据权利要求2所述的一种基于卫星降水数据的降雨侵蚀力计算方法，其特征在于，步骤4中所述根据步骤3中得到的校正系数对步骤2中得到的年平均降雨侵蚀力数据进行校正的具体过程为：在计算得到的年平均降雨侵蚀力数据基础上，按照不同子气候带类型将每个点位上的数值乘以对应的校正系数，得到不同子气候带上校正后的降雨侵蚀力数据结果，然后将不同子气候带降雨侵蚀力校正结果进行拼接，最后得到校正后的全球范围内降雨侵蚀力计算结果。
一种基于卫星降水数据的降雨侵蚀力计算方法
技术领域
[0001] 本发明属于降雨侵蚀技术领域，具体涉及一种基于卫星降水数据的降雨侵蚀力计算方法。
背景技术
[0002] 在全球气候变化背景下，全球范围内的水力侵蚀呈现加剧趋势，准确刻画降雨侵蚀力分布情况，有助于了解不同地区受到水力侵蚀影响的程度，为各国土地利用规划和管理提供参考。获取精准的降水数据是计算降雨侵蚀力的前提，当前对降雨量的测量，主要由分布在全球各地的气象站、雷达和卫星等实现。基于大量气象站观测的降水数据，具有高精度等特点，借助空间插值等方法，可以较为精确的描绘中小尺度范围内的降雨侵蚀力分布情况。但是，在气象站分布稀疏的地区，依靠空间插值方法进行降雨量估算将降低计算结果的准确性。而气象站的分布受限于各国的自然地理条件以及经济情况，难以实现全球覆盖。因此，传统站点插值得到的降雨侵蚀力数据受限于站点分布密度，在大尺度上无法准确反映出降雨侵蚀力分布情况。随着科技的发展，卫星降水数据为更大尺度上的降雨监测提供了可能，该数据具有更大的空间范围和实施观测等优点，弥补了站点数据的不足，但受限于卫星估算降雨量的技术方法，已有研究表明卫星降水数据对于强降雨估算值会偏小，进一步导致整体上计算得到的降雨侵蚀力值偏小。并且以往研究中对于侵蚀性降雨事件的评定标准不一，使得降雨侵蚀力计算结果误差较大。因此如何提高降雨侵蚀力计算结果的准确性是当前需要解决的技术问题。
发明内容
[0003] 本发明的目的在于，提出一种基于卫星降水数据的降雨侵蚀力计算方法，以解决上述技术问题。
[0004] 本发明是通过以下技术方案实现的：
[0005] 本发明提出一种基于卫星降水数据的降雨侵蚀力计算方法，所述方法包括以下步骤：
[0006] 步骤1、获取数据：获取卫星降水数据、站点插值降雨侵蚀力数据以及气候带数据；
[0007] 步骤2、计算降雨侵蚀力：根据获取的卫星降水数据，对单场降雨事件是否为侵蚀性降雨事件进行判定，若判定为侵蚀性降雨事件，则对侵蚀性降雨事件进行降雨侵蚀力计算，进而计算得到年平均降雨侵蚀力数据；
[0008] 步骤3、计算校正系数：将步骤2中得到的年平均降雨侵蚀力数据与步骤1中获取的站点插值降雨侵蚀力数据在气候带数据上进行一元线性回归分析，通过计算得到校正系数，并对校正系数进行可靠性验证；
[0009] 步骤4、结果校正：根据步骤3中得到的校正系数对步骤2中得到的年平均降雨侵蚀力数据进行校正，得到校正后的降雨侵蚀力数据。
[0010] 进一步的是，步骤1中所述气候带数据包括多个不同类型的子气候带数据，分别为：热带气候带数据；半干旱草原气候带数据；沙漠气候带数据；海洋性气候带数据；地中海气候带数据；亚热带湿润气候带数据；温带大陆性湿润气候带数据；温带季风气候带数据。
[0011] 进一步的是，步骤2中所述对单场降雨事件是否为侵蚀性降雨事件进行判定的具体过程为：按下列条件进行判定：
[0012] (1)单场降雨总量P总超过12.7mm；
[0013] (2)单场降雨任意30min内降雨量P30超过6.35mm；
[0014] (3)如果连续降雨在6h内降雨量P6h不超过1.27mm，则将该场连续降雨事件从第6h处划分为两场降雨事件，再分别利用条件(1)和(2)进行判断；
[0015] 符合上述条件之一的降雨事件则判定为侵蚀性降雨事件。
[0016] 进一步的是，步骤2中所述降雨侵蚀力的计算公式为：
[0017]
[0018] 其中，R表示降雨侵蚀力，单位为MJ mm h‑1ha‑1；E表示降雨总动能，单位为MJ ha‑1；‑1
I30表示30min内最大降雨强度，单位为mm h ；n表示侵蚀性降雨事件次数；
[0019] 降雨总动能E的计算公式为：
[0020]
[0021] 其中，k表示将单场侵蚀性降雨事件划分为k个片段；vr为第r个片段的降雨量；er为‑1 ‑1第r个片段的降雨动能，单位为MJ ha mm ；
[0022] 根据RUSLE2中的降雨能量计算公式，则er计算公式为：
[0023]
[0024] 其中，ir表示第r个片段降雨强度的大小，单位为mm h‑1；e为自然常数；所述年平均降雨侵蚀力的计算公式为：
[0025]
[0026] 其中，Rm为年平均降雨侵蚀力，单位为MJ mm h‑1ha‑1yr‑1；m为年份，m≥2。
[0027] 进一步的是，步骤3中所述校正系数的具体计算过程为：
[0028] 读取步骤2中得到的年平均降雨侵蚀力数据以及步骤1中获取的站点插值降雨侵蚀力数据，根据不同类型的子气候带数据对计算得到的年平均降雨侵蚀力数据以及站点插值降雨侵蚀力数据分别进行裁剪；对同一子气候带类型的不同降雨侵蚀力数据进行一元线性回归分析，选取年平均降雨侵蚀力数据作为自变量x，站点插值降雨侵蚀力数据作为因变量y，进行回归分析，其一元线性回归方程为：
[0029] y＝θ0+θ1x (5)
[0030] 其中，θ0为截距项；θ1为校正系数；
[0031] 根据最小二乘法计算求得θ0与θ1，计算过程为：
[0032] 设有q组样本点：(xa，ya)，a＝1，…，q；
[0033] 对每个xa都有一个线性模型预测值 公式为：
[0034]
[0035] 线性模型预测值 与真实值ya之间的差值,称为预测误差 其计算公式为：
[0036]
[0037] 通过下述公式(8)来调整θ0与θ1的取值，使得预测误差 最小：
[0038]
[0039] 最终得到该子气候带类型的校正系数θ1；
[0040] 所述对校正系数进行可靠性验证的具体过程为：
[0041] 计算一元线性回归方程的决定系数R2：
[0042]
[0043] 其中，SSR为回归平方和；SST为总平方和；yi为真实值； 为预测值； 为y平均值；
[0044] 当R2大于0.5时，表明校正系数通过验证。
[0045] 进一步的是，步骤4中所述根据步骤3中得到的校正系数对步骤2中得到的年平均降雨侵蚀力数据进行校正的具体过程为：在计算得到的年平均降雨侵蚀力数据基础上，按照不同子气候带类型将每个点位上的数值乘以对应的校正系数，得到不同子气候带上校正后的降雨侵蚀力数据结果，然后将不同子气候带降雨侵蚀力校正结果进行拼接，最后得到校正后的全球范围内降雨侵蚀力计算结果。
[0046] 本发明的有益效果是：本发明所述方法在卫星降水数据的基础上，结合不同类型子气候带的校正系数，对降雨侵蚀力数据进行校正，减小了由卫星降水数据精度问题所导致的降雨侵蚀力计算偏差；提高了降雨侵蚀力计算结果的准确性。
[0047] 下面结合附图和具体实施方式对本发明作进一步详细说明。
附图说明
[0048] 图1为本发明所述方法流程示意图；
[0049] 图2为实施例一中在不同子气候带类型上的一元线性回归分析结果示意图。
具体实施方式
[0050] 本发明提出一种基于卫星降水数据的降雨侵蚀力计算方法，如图1所示，包括以下步骤：
[0051] 步骤1、获取数据：获取卫星降水数据、站点插值降雨侵蚀力数据以及气候带数据。
[0052] 卫星降水数据为通过在美国宇航局(NASA)下载得到的GPM卫星的Level 3(Final Run)降水数据；站点插值降雨侵蚀力数据为通过在欧洲土壤数据中心(ESDAC)下载的全球范围的站点插值降雨侵蚀力数据；气候带数据为通过在国家地球系统科学数据中心下载的‑Geiger气候带数据；气候带数据包括多个不同类型的子气候带数据，分别为：热带气候带数据；半干旱草原气候带数据；沙漠气候带数据；海洋性气候带数据；地中海气候带数据；亚热带湿润气候带数据；温带大陆性湿润气候带数据；温带季风气候带数据。
[0053] 步骤2、计算降雨侵蚀力：根据获取的卫星降水数据，对单场降雨事件是否为侵蚀性降雨事件进行判定，若判定为侵蚀性降雨事件，则对侵蚀性降雨事件进行降雨侵蚀力计算，进而计算得到年平均降雨侵蚀力数据。
[0054] 首先根据获取的卫星降水数据，计算得到单场降雨总量P总以及单场降雨任意30min内降雨量P30，然后根据下列条件判定单场降雨事件是否为侵蚀性降雨事件：
[0055] (1)单场降雨总量P总超过12.7mm；
[0056] (2)单场降雨任意30min内降雨量P30超过6.35mm；
[0057] (3)如果连续降雨在6h内降雨量P6h不超过1.27mm，则将该场连续降雨事件从第6h处划分为两场降雨事件，再分别利用条件(1)和(2)进行判断。
[0058] 符合上述条件之一的降雨事件则判定为侵蚀性降雨事件，对侵蚀性降雨事件进行降雨侵蚀力计算。
[0059] 降雨侵蚀力R的计算公式为：
[0060]
[0061] 其中，R表示降雨侵蚀力，单位为MJ mm h‑1ha‑1；E表示降雨总动能，单位为MJ ha‑1；‑1
I30表示30min内最大降雨强度，单位为mm h ；n表示侵蚀性降雨事件次数；I30通过逐个比对该场侵蚀性降雨中每个30min的降雨强度，筛选得到的30min内最大降雨强度。
[0062] 降雨总动能E的计算公式为：
[0063]
[0064] 其中，k表示将单场侵蚀性降雨事件划分为k个片段；vr为第r个片段的降雨量；er为‑1 ‑1第r个片段的降雨动能，单位为MJ ha mm 。
[0065] 根据RUSLE2中的降雨能量计算公式，则er计算公式为：
[0066]
[0067] 其中，ir表示第r个片段降雨强度的大小，单位为mm h‑1；e为自然常数。
[0068] 进一步，根据降雨侵蚀力R，计算年平均降雨侵蚀力Rm的公式为：
[0069]
[0070] 其中，Rm为年平均降雨侵蚀力，单位为MJ mm h‑1ha‑1yr‑1；m为年份，m≥2。
[0071] 步骤3、计算校正系数：将步骤2中得到的年平均降雨侵蚀力数据与步骤1中获取的站点插值降雨侵蚀力数据在气候带数据上进行一元线性回归分析，通过计算得到校正系数，并对校正系数进行可靠性验证。
[0072] 利用R语言将上述计算得到的年平均降雨侵蚀力数据输出为ASCII格式的文件，然后利用ArcGIS软件读取上述ASCII文件以及步骤1中获取的站点插值降雨侵蚀力数据，根据不同类型的子气候带数据对计算得到的年平均降雨侵蚀力数据以及站点插值降雨侵蚀力数据分别进行裁剪；然后在对同一子气候带类型的不同降雨侵蚀力数据进行一元线性回归分析。
[0073] 选取年平均降雨侵蚀力数据作为自变量x，站点插值降雨侵蚀力数据作为因变量y，进行回归分析，其一元线性回归方程为：
[0074] y＝θ0+θ1x (5)
[0075] 其中，θ0为截距项；θ1为校正系数。
[0076] θ0与θ1为待定系数，根据最小二乘法计算求得θ0与θ1，计算过程为：
[0077] 设有q组样本点：(xa，ya)，a＝1，…，q；
[0078] 对每个xa都有一个线性模型预测值 公式如下所示：
[0079]
[0080] 线性模型预测值 与真实值ya之间的差值,称为预测误差 其计算公式为：
[0081]
[0082] 为使得预测值 与真实值ya更接近，故需要使得预测误差 更小，通过下述公式(8)来调整θ0与θ1的取值，使得预测误差 最小：
[0083]
[0084] 根据公式计算确定的θ0与θ1即为所求，从而得到该子气候带类型的校正系数θ1。
[0085] 为验证所求得的校正系数的可靠性，需要对一元线性回归方程模型计算其决定系2 2 2
数R ，R 大小表征为：模型可以解释为多大程度是自变量导致因变量的改变。即例如R 为
2
0.7，则可以解释为模型有70％是自变量导致因变量改变的。一般来说，当R 大于0.5时，可以认为一元线性回归结果较好，即通过可靠性验证。
[0086] 决定系数R2的计算公式为：
[0087]
[0088] 其中，SSR为回归平方和；SST为总平方和；yi为真实值； 为预测值； 为y平均值。
[0089] 根据上述计算方法，依次计算不同子气候带类型上的校正系数。
[0090] 步骤4、结果校正：根据步骤3中得到的校正系数对步骤2中得到的年平均降雨侵蚀力数据进行校正，得到校正后的降雨侵蚀力数据。
[0091] 在计算得到的年平均降雨侵蚀力数据基础上，按照不同子气候带类型将每个点位上的数值乘以对应的校正系数，得到不同子气候带上校正后的降雨侵蚀力数据结果，然后利用ArcGIS软件将不同子气候带降雨侵蚀力校正结果进行拼接，则得到校正后的全球范围内降雨侵蚀力计算结果。
[0092] 实施例一
[0093] 本实施例为上述方法的具体应用实例。
[0094] 本实施例选用的卫星降水数据为时间精度为30min，空间精度为0.1°的2001‑2020年间GPM卫星降水数据；按照上述方法步骤计算得到全球20年的年平均降雨侵蚀力数据；站点插值降雨侵蚀力数据为通过在欧洲土壤数据中心(ESDAC)下载的全球范围的站点插值降雨侵蚀力数据；气候带数据为通过在国家地球系统科学数据中心下载的 ‑Geiger气候带数据。
[0095] 在ArcGIS中根据八个子气候带数据对年平均降雨侵蚀力数据以及站点插值降雨侵蚀力数据分别进行裁剪，分别为：A，热带气候；Bs，半干旱草原气候；Bw，沙漠气候；Cf，海洋性气候；Cs，地中海气候；Cw，亚热带湿润气候；Df，温带大陆性湿润气候；Dw，温带季风气候。
[0096] 将年平均降雨侵蚀力数据与站点插值降雨侵蚀力数据在不同子气候带类型上进行一元线性回归分析，依次计算不同子气候带类型上的校正系数，结果如图2所示。本实施例在R语言里采用lm(y～x‑1)函数，从而舍弃了截距θ0，得到过原点的一元线性回归方程：Y2
＝θ1*x，θ1即为本实施例所计算的校正系数。决定系数R 利用R语言中的ggplot2工具计算得到，结果如表1所示。
[0097] 表1
[0098]
[0099] 根据上述计算，最后得到如图2所示的回归结果、回归方程以及一元线性回归方程2
模型的决定系数R。
[0100] 根据上述计算得到的校正系数，对全球20年的年平均降雨侵蚀力数据进行校正，得到校正后的降雨侵蚀力数据，利用ArcGIS软件出图，得到校正后的全球20年的年平均降雨侵蚀力数据分布图。
[0101] 为验证本发明所述方法的结果优越性，在本实施例得到校正后的全球20年的年平均降雨侵蚀力数据分布图后，紧接着利用ArcGIS软件查看本实施例计算结果(GPM)、站点插值降雨侵蚀力数据结果(GloREDa)以及以往基于卫星降水数据计算所得降雨侵蚀力数据结果(CMORPH)的平均值和标准差，同时利用ArcGIS软件将GPM结果与CMORPH结果分别与GloREDa结果相减，最后再在R语言中将GPM与CMORPH分别与GloREDa作相关性分析，选取的是斯皮尔曼(Spearman)相关性分析，其计算公式如下：
[0102]
[0103] 其中，β为Spearman秩相关系数；s为对象数量；di为对应变量的秩之差，即两个变量分别排序后成对的变量位置(等级)差。
[0104] β为本实施例的秩相关系数结果，其结果如表2所示。
[0105] 表2
[0106]
[0107] 从表2中可以看出本实施例结果与以往基于卫星降水数据计算所得降雨侵蚀力数据结果相比发现，本实施例结果更接近于站点插值降雨侵蚀力数据结果，准确度更高。
[0108] 最后应说明的是，以上所述仅用以说明本发明的技术方案而非限制，尽管参照较佳布置方案对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的精神和范围。
