CN109344865B - 一种多数据源的数据融合方法 - Google Patents

一种多数据源的数据融合方法 Download PDF

Info

Publication number
CN109344865B
CN109344865B CN201810974820.9A CN201810974820A CN109344865B CN 109344865 B CN109344865 B CN 109344865B CN 201810974820 A CN201810974820 A CN 201810974820A CN 109344865 B CN109344865 B CN 109344865B
Authority
CN
China
Prior art keywords
data
field
daily
cfsr
time
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810974820.9A
Other languages
English (en)
Other versions
CN109344865A (zh
Inventor
谢刚
彭岩波
韩子叻
宋杰
邓杰
王国强
王溥泽
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shandong Institute of ecological environment planning
Beijing Normal University
Original Assignee
Shandong Academy For Environmental Planning
Beijing Normal University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shandong Academy For Environmental Planning, Beijing Normal University filed Critical Shandong Academy For Environmental Planning
Priority to CN201810974820.9A priority Critical patent/CN109344865B/zh
Publication of CN109344865A publication Critical patent/CN109344865A/zh
Application granted granted Critical
Publication of CN109344865B publication Critical patent/CN109344865B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/25—Fusion techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/24—Classification techniques

Abstract

本发明实施例提出了一种多数据源的数据融合方法，包括：基础数据处理步骤，用于对台站气象数据、CFSR气象数据、Princeton大气驱动数据、遥感数据进行处理以获取处理后的基础数据；数据集建立步骤，用于通过根据处理后的基础数据建立近地面气温场、近地面相对湿度场、近地面风场、近地面气压场、降水场、辐射场。

Description

一种多数据源的数据融合方法
技术领域
本发明涉及数据处理技术领域，尤其是涉及一种多数据源的数据融合方法，其典型的可用于气象数据的融合。
背景技术
随着社会的发展，越来越多领域都开始使用数据分析和数据处理技术。在使用海量数据的很多领域，都需要对多个数据源得到的数据进行融合，然后再进行后续的分析和处理。气象数据是一种极其典型的海量数据集，气象数据要综合多种数据源的海量数据，因此必须先对数据进行融合以使这些数据可以统一使用。举例来说，我国所采用的气象数据通常包括：中国大陆地区台站观测资料、CFSR再分析资料、Princeton大气驱动资料、GEWEXSRB下行短波辐射资料。
在生成数据集之前需要对海量数据进行融合，但是现有技术中缺少可靠、高效的手段对这些海量数据进行融合，导致最终在生成数据集时存在很多问题。
发明内容
针对当前的数据融合技术不够完善的问题，本发明实施例提出了一种为了实现上述目的，本发明实施例提供了一种多数据源的数据融合方法，包括：
基础数据处理步骤，用于对台站气象数据、CFSR气象数据、Princeton大气驱动数据、遥感数据进行处理以获取处理后的基础数据；
数据集建立步骤，用于通过根据处理后的基础数据建立近地面气温场、近地面相对湿度场、近地面风场、近地面气压场、降水场、辐射场。
进一步的，所述基础数据处理步骤具体包括：
台站气象数据处理子步骤，用于获取来自于中国气象局气象信息中心常规气象要素740站的观测数据，其中所述观测数据中包括以下的至少一种观测变量：近地面1.5米气温、气压、相对湿度、近地面10米风速、累积降水和日照时数(用于产生辐射场)；并将中国大陆地区台站观测资料分成四个分区的台站观测资料；
CFSR气象数据再分析子步骤：用于获取NCEP CFSR再分析资料；其中CFSR再分析资料中包括：近地面气温、相对湿度、风速和气压变量的3小时分辨率的数据，以为大气、海洋、陆地和海冰模式提供初始场；
Princeton大气驱动数据再分析子步骤：基于NCEP-1再分析资料，再利用月气候资料对NCEP-1再分析资料进行订正，然后进行空间降尺度；其中Princeton大气驱动数据中至少包括：近地面气温、相对湿度、风速、气压和短波辐射数据；
遥感数据处理子步骤：用于对CMORPH降水资料进行降维处理，并对GWEXE SRB下行短波辐射资料中的地表下行短波辐射资料进行处理。
进一步的，所述数据集建立步骤中通过以下方法建立近地面气温场：
步骤11、建立趋势面，在每个有观测的时刻，建立如下的趋势面模型，
1958：1978：t(x，y)＝f(x，y)+β·z(x，y)+ε(x，y)
1979-2010：t(x，y)＝f(x，y)+β1·z(x，y)+β2·tcfsr(x，y)+ε(x，y) (2-11)
式中t代表近地面2米气温，z为高程，tcfsr是CFSR再分析数据在点的线性插值；其中(x，y)是经纬度，{ψr，r＝1，2，...，p}是给定的协变量函数(例如高程)，βr是回归系数，f是2阶薄板平滑样条函数，ε为误差，其中误差ε是独立同分布的；
其中，某一时刻的驱动变量u分解为
Figure GDA0003469179360000021
其中(x，y)是经纬度，{ψr，r＝1，2，...，p}是给定的协变量函数，βr是回归系数，f是2阶薄板平滑样条函数；ε为误差，是独立同分布的；公式中的f(x，y)和
Figure GDA0003469179360000022
合称为变量的趋势面；只要样条函数的系数和回归系数确定了，趋势面相对于经纬度的函数关系就确定了；则在任一点的插值只需将该点的经纬度代入趋势面的函数即可获得；在本步骤中驱动变量u为近地面2米气温t；其中2阶薄板平滑样条函数f是针对于在经度{x1，x2，...xn}和纬度{y1，y2，...yn}上的n维观测向量定义的；其表现为
Figure GDA0003469179360000031
和限制条件
Figure GDA0003469179360000032
其中x1,x2,...xn和y1,y2,...yn是样条系数；d1、d2、d3为一阶系数，ci为二阶系数；i＝1,2...,n；
由公式(2-1)、公式(2-2)可知，2阶薄板平滑样条函数f的自由度与观测的个数相同；因此在计算样条系数和回归系数时，不但要使趋势面靠近观测，也要使样条尽量平滑以控制趋势面在无观测地区的误差；在实际操作中，这些系数被取为使得以下目标函数的最小的值，
Figure GDA0003469179360000033
其中第一项代表误差的规模，第二项
Figure GDA0003469179360000034
J2(f)代表2阶薄板平滑样条函数f的光滑性；λ是光滑参数，用于平衡误差规模和光滑性；zi为高程；j＝0,1,2；
估计样条系数和回归系数的关键技术是对光滑参数λ的估计；因为λ确定后，目标函数(2-2)只是样条系数和回归系数的2次函数，使它在约束条件(2-3)下极小的系数有显示表达；
基于最小交叉验证原则对光滑参数λ进行验证：
步骤12、对趋势面进行订正，具体包括：
用简单克里金方法对趋势面的残差场插值，得到插值结果为
Figure GDA0003469179360000041
其中,||x,y；xi,yi||表示点(x,y)和(xi,yi)之间的欧氏距离，c是一个单调非负函数的协方差函数。然后将残差场的插值叠加到趋势面后就得到最后的插值场
Figure GDA0003469179360000042
其中协方差函数c的估计采用以下方法：将观测集合按两点间的距离分为若干点对的集合，在每一个集合内计算点对的协方差，以构造一个随距离变化的序列；用一维薄板平滑样条拟合这个序列就可到协方差函数c；
步骤13、将得到的近地面气温场与CFSR气象数据、Princeton大气驱动数据、遥感数据进行对比以确定数据的准确性。
处理子步骤：用于对CMORPH降水资料进行降维处理，并对GWEXE SRB下行短波辐射资料中的地表下行短波辐射资料进行处理。
进一步的，所述数据集建立步骤中通过以下方法建立近地面相对湿度场：
步骤21、对于相对湿度q，在每个有观测的时刻，建立如下的趋势面模型，
1958-1978:q(x，y)＝f(x，y)+ε(x，y)
1979-2010:q(x，y)＝f(x，y)+β1qcfsr(x，y)+ε(x，y) (2-12)
其中qcfsr是CFSR再分析数据在点的线性插值，其中2阶薄板平滑样条函数f的计算步骤与步骤11相同；
步骤22、对趋势面进行订正，订正方式与与步骤12相同；
步骤23、对日尺度插值结果和6小时尺度的插值结果进行时间降尺度，具体包括：
在每一日的监测数据中确定有监测数据的时刻，并在每两个相邻的监测数据之间进行插值，插值的为相邻时刻的驱动场的平均值；
利用插值后的监测数据对辅助资料进行订正以达到时间降尺度的目的；其中该辅助资料为以下的一种：CFSR再分析资料、Princeton大气驱动资料、GEWEX SRB下行短波辐射资料；
其中利用插值后的监测数据对辅助资料进行订正具体包括：对于风速的气象变量，获取其辅助资料，将3小时的辅助资料除以它的日平均再乘以插值的日资料；
步骤24、将得到的地面相对湿度场与CFSR气象数据、Princeton大气驱动数据、遥感数据进行对比以确定数据的准确性。
进一步的，所述数据集建立步骤中通过以下方法建立近地面风场：
步骤31、对于近地面风场w，在每个有观测的时刻，采用如下模型：
w(x，y)＝f(x，y)+ε(x，y) (2-13)
步骤32、对日尺度插值结果和6小时尺度的插值结果进行时间降尺度，具体包括：
在每一日的监测数据中确定有监测数据的时刻，并在每两个相邻的监测数据之间进行插值，插值的为相邻时刻的驱动场的平均值；
利用插值后的监测数据对辅助资料进行订正以达到时间降尺度的目的；其中该辅助资料为以下的一种：CFSR再分析资料、Princeton大气驱动资料、GEWEX SRB下行短波辐射资料；
其中利用插值后的监测数据对辅助资料进行订正具体包括：对于风速的气象变量，获取其辅助资料，将3小时的辅助资料减去它的日平均再加上插值的日资料。
进一步的，所述数据集建立步骤中通过以下方法建立近地面气压场：
步骤41、对于近地面风场w，在每个有观测的时刻，采用如下模型：
p(x，y)＝f(x，y)+s(z(x，y))+ε(x，y) (2-14)
其中s是一维平滑样条函数；
步骤42、对日尺度插值结果和6小时尺度的插值结果进行时间降尺度，具体包括：
在每一日的监测数据中确定有监测数据的时刻，并在每两个相邻的监测数据之间进行插值，插值为相邻时刻的驱动场的平均值；
利用插值后的监测数据对辅助资料进行订正以达到时间降尺度的目的；其中该辅助资料为以下的一种：CFSR再分析资料、Princeton大气驱动资料、GEWEX SRB下行短波辐射资料；
其中利用插值后的监测数据对辅助资料进行订正具体包括：对于气压的气象变量，获取其辅助资料，将3小时的辅助资料减去它的日平均再加上插值的日资料。
进一步的，所述数据集建立步骤中通过以下方法建立降水场：
步骤51、在整个1958-2010年用日降水资料建立趋势面，在以下时间段的趋势面模型分别为
1958-1978：r(x，y)＝f(x，y)+ε(x，y)
1979-1997：r(x，y)＝f(x，y)+α·rcfsr(x，y)+ε(x，y)
1998-2010：r(x)＝f(x，y)+α·rcmorph(x，y)+ε(x，y) (2-15)
其中r代表日观测降水量，rcdsr(x.y)和rcmorph(x，y)分别为日平均的CFSR再分析资料和CMORPH降水产品在(x，y)处的线性插值；
步骤52、对日尺度插值结果和6小时尺度的插值结果进行时间降尺度，具体包括：
在每一日的监测数据中确定有监测数据的时刻，并在每两个相邻的监测数据之间进行插值，插值的为相邻时刻的驱动场的平均值；
利用插值后的监测数据对辅助资料进行订正以达到时间降尺度的目的；其中该辅助资料为以下的一种：CFSR再分析资料、Princeton大气驱动资料、GEWEX SRB下行短波辐射资料；
其中利用插值后的监测数据对辅助资料进行订正具体包括：对于风速的气象变量，获取其辅助资料，将3小时的辅助资料除以它的日平均再乘以插值的日资料。
进一步的，所述数据集建立步骤中通过以下方法建立辐射场：利用台站观测资料中的日照时数、日最高气温、日最低气温、日平均相对湿度，建立下行短波辐射。日照时数是指太阳在一地实际照射地面的时数，即地面观测地点受到太阳直接辐射辐照度等于和大于120W/m2的累计时间；
具体包括：
步骤61、建立日照时数在两个时间段的回归模型：
1958-1982和2008-2010：s(x，y)＝f(x，y)+β1td(x，y)+β2F(q(x，y))+ε(x，y)
1983-2007：s(x，y)＝f(x，y)+β1td(x，y)+β2F(q(x，y))+β3G(x，y)+ε(x，y) (2-16)
其中s为日照时数，td为温度在北京时间的日变化，q为在北京时间的日平均相对湿度，；
Figure GDA0003469179360000061
G(x，y)是短波辐射遥感产品GEWEX SRB在处的插值；利用5公里格点的经纬度信息、温度日变化、日均相对湿度以及GEWEX SRB产品生成5公里格点场的日尺度日照时数；
步骤62、利用太阳辐射模型计算5公里3小时的下行短波辐射；具体包括：
根据前面已经建立的近地面气温场，相对湿度场和气压场，计算日平均气温、日平均相对湿度和日平均气压，并结合日照时数场，计算日尺度的下行短波辐射；
根据太阳高度角的变化，将下行短波辐射的日值分配到3小时的时间尺度上，最后，生成时空分辨率为5公里3小时的下行短波辐射场，并将北京时间转换成世界时间；
其中下行长波辐射场是由短波辐射、近地面气温、相对湿度和气压场直接计算生成。
进一步的，所述方法还包括：
分类数据的数字化处理步骤，用于对土壤数据进行分类处理；具体包括：
土壤平均粒径计算步骤，根据以下公式计算土壤平均粒径：
Figure GDA0003469179360000071
其中Si是第i种土壤粒径的序数，ri是这种粒径的土壤所占的比例；
地表植物的平均根深计算步骤，利用地表植物的序数和这种地表植物所占的比例，计算地表植物的平均根深；
数据聚类分析步骤，用于将数字化的土壤和地表植物的数字化参数与其他流域数据/参数整合在一起，进行多源数据的聚类分析，以获取流域的降水特点与土壤平均粒径、植被类型、海拔高度之间的关系；
根据上中下游的不同的地理条件选择不同的多源降水数据融合方案以及水文模型的参数化方案；通过不同的降水数据和参数化方案的模拟效果进行比较，以获取与子流域属性相关的最优的降水数据融合以及水文模型的参数化方案。
本发明的技术方案具有以下优势：
上述方案提出了一种多数据源的数据融合方法，其融合了中国大陆地区台站观测资料、CFSR再分析资料、Princeton大气驱动资料和GEWEX SRB下行短波辐射资料，构造中国大陆地区1958-2010年逐3小时5公里×5公里分辨率的大气驱动数据集。
附图说明
通过下面结合附图对本发明的一个优选实施例进行的描述，本发明的技术方案及其技术效果将变得更加清楚，且更加易于理解。其中：
图1为本发明实施例的流程图。
具体实施方式
以下将结合所附的附图对本发明的一个优选实施例进行描述。
本发明实施例通过融合中国大陆地区台站观测资料、CFSR再分析资料、Princeton大气驱动资料和GEWEX SRB下行短波辐射资料，构造中国大陆地区1958-2010年逐3小时5公里×5公里分辨率的大气驱动数据集(称为BNU中国地区陆面模型大气驱动场)，包括近地面气温、相对湿度、风速、气压、降水、下行短波辐射和下行长波辐射这7个要素，并对BNU中国地区陆面模型大气驱动场的结果进行评估。本发明实施例示例性的采用交叉验证的评估方法对2003年每月3号的6小时数据进行验证。
基础数据
1.1台站观测资料及其预处理
本发明实施例用于完成1958-2010年中国大陆区域驱动数据的研制应用的观测来自于中国气象局(China Meteorological Administration,CMA)气象信息中心常规气象要素740站的观测数据(其中2007-2010年的产品是在NMIC计算完成)。涉及的观测变量有：近地面1.5米气温、气压、相对湿度、近地面10米风速、累积降水和日照时数(用于产生辐射场)。由国家气象信息中心(National Meteorological Information Center of China,NMIC)提供的观测数据是经过质量控制的并且误差控制在2％以内(任芝花等,2007；中国气象局,2003)。为了计算更准确也更方便，本发明实施例基于地形和站点观测密度的因素考虑，并将中国大陆地区台站观测资料分成四个分区的台站观测资料。
这套观测数据的不同变量在不同时段的观测频次也不同：其中降水和日照时数这两个变量在整个1958-2010年驱动数据的制备过程中是日观测数据。对于近地面气温、相对湿度、风速和气压的观测数据，1958-1989年是日尺度观测；1990-1997年是6小时频率的观测；1998-2010年的观测资料是3小时频率。
此外，中国地区的近地面气温、相对湿度、风速和气压的台站观测数据在1998-2006年期间虽然是3小时观测，但是各时刻之间的观测数量很不均匀(3，9，15，21时的观测数据量比0，6，12，18时的数据量少)。本发明实施例利用薄板平滑样条模型的理论(见1.2.1)来将这个时间段3，9，15，21时刻的观测补全。
1.2再分析资料
1.2.1CFSR再分析资料
本发明实施例应用NCEP CFSR再分析资料中的近地面气温、相对湿度、风速和气压变量的3小时分辨率的数据(0、6、12、18时是用再分析数据，3、9、15、21是用的预报数据)。NCEP CFSR资料是NCEP最新的一套全球耦合再分析数据(Saha et al.,2010；2013)，发展该数据集的主要目标就是为大气、海洋、陆地和海冰模式提供初始场。
CFSR再分析资料中的大气变量的空间分辨率为0.3125°×0.3125°(约38Km)，37个气压层，时间覆盖范围是1979-2011，同化的时间间隔是6小时并且有逐小时的预报产品。
CFSR资料在整个时间段内都对卫星观测的辐射率进行了同化，这是NCEP首次将卫星辐射率直接同化进它的全球再分析产品中；CFSR资料采用的模式分辨率是T382，水平分辨率约为38公里，垂直分层是64层等面，最高0.266hPa，相比于NCEP之前的模式有较大的提高；CFSR在产生6小时初猜场时还耦合了海洋模式，并加入了交互的海冰模式，在辐射的参数化方案中考虑了CO2、气溶胶和其它痕量气体的变化；CFSR的陆面模块采用Noah陆面模式。
本发明实施例用到了1979-2010年CFSR数据集中的近地面气温、气压、相对湿度、风速和降水变量。
1.2.2Princeton大气驱动资料
Princeton陆面模式驱动数据是普林斯顿大学发展的一套专门用于驱动陆面水文和能量收支模型的长时间序列的覆盖全球的近地表气象要素数据集(Sheffiled et al.,2004，2006)。时间分辨率是3小时，空间分辨率是1°×1°。时间覆盖范围是1948-2006年。该数据集包括了7个变量：近地面气温、气压、比湿、近地面风速、向下短波辐射和向下长波辐射和降水率。Princeton数据集制作过程中需要的数据集有：NCEP/NCAR(NCEP-1)再分析资料、CRU(Climatic Research Unit monthly climate variables)月气候资料、GPCP(Global capitationion Climatology Project daily precipitation)日尺度的降水资料、TRMM(Tropical Rainfall Measuring Mission)3小时降水资料和NASA Langley的月尺度地表辐射收支数据。
Princeton数据是基于NCEP-1再分析资料，利用CRU的月尺度数据进行订正，再做空间降尺度得到。本发明实施例用到了Princeton数据集中的近地面气温、相对湿度、风速、气压和短波辐射数据。
1.3遥感资料
1.3.1CMORPH降水资料
本发明实施例在建立1998年以后的中国区域降水格点场时，用到了美国NOAA CPC(National Oceanic and Atmospheric Administration Climate Prediction Center)开发的CMORPH(CPC Morphing Technique)高时空分辨率的卫星反演降水产品。这种降水估测技术(Joyce et al.,2004)全部采用了低轨道卫星的微波观测估计，而低轨道卫星微波观测降水估计特点是用来自静止卫星红外数据来进行空间传递信息。目前为止CMORPH降水资料中，融合了DMSP 13,14&15(SSM/I),NOAA-15,16,17&18(AMSU-B)，AMSR-E，TMI aboardNASA's Aqua和TRMM等卫星上的被动微波降水估计。这些降水估计利用Ferraro(1997)反演SSM/I,Ferraro et al.(2000)反演AMSU-B and Kummerow et al.(2001)反演TMI数据用到的方法。值得注意的是这种技术不是一种降水估计方法，而是现有的微波降水估计的方法的融合后的平均。因此，这种方法非常灵活，来自任意一颗卫星的微波降水估计都可以包括进来。
原始CMORPH资料的空间分辨率是8km，时间分辨率是30分钟，本发明实施例所采用的CMORPH数据的空间分辨率是0.25°×0.25°，时间分辨率是3小时。
1.3.2GWEXE SRB下行短波辐射资料
美国航天局的GEWEX SRB(Global Energy and Water Exchanges SurfaceRadiation Budget)项目是GEWEX对辐射研究的一个重要组成部分，该项目旨在产生能够预测短暂的气候变化和长期的气候趋势的精确的地表、大气层顶、大气短波长波辐射通量。本发明实施例使用该项目提供的地表下行短波辐射资料。该数据通过GEWEX-SRB辐射方法得到，本发明实施例使用的短波数据集的版本是3.0，资料的时间长度是1983-2007年，空间分辨率是1°×1°。
方法原理
2.1驱动变量趋势面的建立
在本发明实施例中，某一时刻的驱动变量u(例如温度)可以假定分解为
Figure GDA0003469179360000101
其中(x,y)是经纬度，{ψr,r＝1,2,...,p}是给定的协变量函数(例如高程)，βr是回归系数，f是2阶薄板平滑样条函数(以下将详细介绍)，ε为误差，是独立同分布的。式右边前两项合称为变量的趋势面。只要样条函数的系数和回归系数确定了，趋势面相对于经纬度的函数关系就确定了。在任一点的插值只需将该点的经纬度代入趋势面的函数即可获得。
一个2阶薄板平滑样条是针对于在经度{x1,x2,...xn}和纬度{y1,y2,...yn}上的n维观测向量定义uo＝[uo(x1，y1)，uo(x2，y2)，...uo(xn，yn)]的，它具有如下形式
Figure GDA0003469179360000111
和限制条件
Figure GDA0003469179360000112
其中{x1,x2,...xn}和{x1,x2,...xn}是样条系数。
由上两式可知，样条f的自由度与观测的个数相同。所以不可用通常的最小2乘法估计样条的系数。因为这做会使估计的样条通过每一个观测点，而对没有观测的地方误差将很大。在估计样条系数和回归系数时，不但要使趋势面靠近观测，也要使样条尽量平滑以控制趋势面在无观测地区的误差。在实际操作中，这些系数被取为使得以下目标函数的最小的值，
Figure GDA0003469179360000113
其中第一项代表误差的规模，第二项
Figure GDA0003469179360000114
J2(f)代表样条的光滑性；λ是光滑参数，用于平衡误差规模和光滑性。其中δ是偏微分符号，d为微分符号。
估计样条系数和回归系数的关键技术是对光滑参数λ的估计。因为λ确定后，目标函数(2-2)只是样条系数和回归系数的2次函数，使它在约束条件(2-3)下极小的系数有显示表达。在本发明实施例中，对λ的估计是基于最小交叉验证原则。
在实际中观测的空间分布往往被设计得比较均匀。相对于去掉某个台站的观测集合来说，可认为被去掉的台站位于观测稀疏的地区。因为趋势面函数的确定用到了交叉验证原则，所以插值在观测稀疏的地区误差最小。另一方面，趋势面与观测的残差往往显著大于观测的仪器误差。这是因为观测集合相对于变量u的空间变化往往比较稀疏，使得有一部分本应属于趋势面的变化归到误差中去了。
2.2趋势面的订正
前一小节提到，用交叉验证原则估计的趋势面与观测的残差往往过大。这在早期被用户们视为薄板平滑样条插值的缺点。而其它插值方法，例如Cressman插值和Barnes插值等，就不存在这个问题(Zheng and Basher 1995)。在本节中，本发明实施例将用残差订正方法克服这一缺点。
首先用简单克里金方法对趋势面的残差场插值，得到在点处的插值结果为
Figure GDA0003469179360000121
其中,||x,y；xi,yi||表示点(x,y)和(xi,yi)之间的欧氏距离，c是一个单调非负函数的协方差函数。然后将残差场的插值叠加到趋势面后就得到最后的插值场
Figure GDA0003469179360000122
在本发明实施例中协方差函数c的估计是在“和的协方差只与两点之间的距离||x,y,xi,yi||有关”的假定下进行的。具体方法如下：将观测集合按两点间的距离分为若干点对的集合，在每一个集合内计算点对的协方差，以构造一个随距离变化的序列；用一维薄板平滑样条拟合这个序列就可得到协方差函数c。
2.3时间降尺度
在某些情况下，观测资料只是一日4次(即0、6、12、18世界时)。这样只能在这4个时刻进行插值。这时其它4个时刻(3、9、15、21世界时)的驱动场将估计为他们的相邻时刻的驱动场的平均。
如果只存在日平均观测资料，本发明实施例只能得到插值的日资料。这时本发明实施例用插值的日资料去订正某个3小时的辅助资料，以达到时间降尺度的目的。这个辅助资料视情况可以是CFSR再分析资料、Princeton大气驱动资料、GEWEX SRB下行短波辐射资料等等。对气温、气压、风速和辐射等变量，将3小时的辅助资料减去它的日平均再加上插值的日资料；对降水和湿度等变量，将3小时的辅助资料除以它的日平均再乘以插值的日资料。这样订正后的辅助资料的日平均与插值的日资料相等。
2.4空间拼接
如果一个区域中的数据较多，用薄板平滑样条拟合该区域的数据计算量会较大。本发明实施例将中国地区分别划分为多个区域。空间插值先在每一个区域上进行，然后再拼接各区域上的插值结果。
2.5评价指标
评价两个插值方案的基本思路是：将观测视为真值的参考值，评价两种插值方案哪一个的插值结果离观测更近。一个评价指标是均方根误差(root-mean-square error，RMSE)
Figure GDA0003469179360000131
这里N是研究区域的所有时刻的台站数量总和，
Figure GDA0003469179360000132
是在(xi,yi)位置上的估计值，它可以是估计的趋势面，也可以是趋势面订正后的插值结果，也可以是再分析资料在台站上的插值结果。另一个评价指标是交叉验证(Cross Validation，CV)，
Figure GDA0003469179360000133
这里
Figure GDA0003469179360000134
的意义与定义(8)中
Figure GDA0003469179360000135
相同，只是下标-i表示点(xi,yi)处的估计值没有用到观测u0(xi,yi)。
RMSE和CV都能代表拟合的优度。RMSE在估计点(xi,yi)的值的时候用到了该点的观测值u0(xi,yi)，但是计算CV时却没有用到。因此从某种意义上来讲，RMSE代表了在台站观测比较稠密的区域的拟合误差，然而CV可以看作是用来评估插值方法在独立数据上的验证。
将粗格点上的再分析资料插值到观测点上会造成很大插值误差。以上用台站数据直接验证再分析资料和本发明实施例插值资料的优劣对再分析资料不太公道。本发明实施例采用以下指标评估再分析资料和本发明实施例插值资料在再分析资料的粗格点上的优劣，
Figure GDA0003469179360000141
这里m是包含观测的一个再分析资料的格，M是区域中这种格的总数，
Figure GDA0003469179360000142
是格m里观测的均值，ur(m)是格m上的再分析数据，
Figure GDA0003469179360000143
是利用本发明实施例的插值方法在格m中所有5公里分辨率的子格点的平均，但格m中观测数据在估计
Figure GDA0003469179360000144
时被剔除了。
中国大陆地区大气驱动数据集的建立及验证
3.1近地面气温场
3.1.1格点场的建立
利用前面的方法建立趋势面，在每个有观测的时刻，建立如下的趋势面模型，
1958：1978：t(x，y)＝f(x，y)+β·z(x，y)+ε(x，y)
1979-2010：t(x，y)＝f(x，y)+β1·z(x，y)+β2·tcfsr(x，y)+ε(x，y) (2-11)
式中t代表近地面2米气温，z为高程，tcfsr是CFSR再分析数据在点的线性插值，其它符号与公式(2-1)相同。同时用1.2.2小节方法对趋势面模型残差进行订正。
其中在1990年之前，利用日观测数据建立模型，1990-1997年，用6小时观测资料建立模型，1998-2010年是用3小时观测资料建立模型。最后对日尺度和6小时时间尺度的插值结果利用1.2.3节描述的方法进行时间降尺度插值。其中1958-1978年借助于3小时的Princeton资料作为辅助资料，1979-1989年采用CFSR资料作为辅助资料。
3.2.2方法评估与精度验证
(1)与其他方法的比较
本部分也采用2003年每月3号的观测数据进行验证。
一种比较流行的空间插值方法是：用原始的再分析资料或遥感反演资料做趋势面，再把趋势面残差场的空间插值结果叠加在趋势面上，以形成最后的格点场(Xie andXiong，2011)。不可否认，这种插值方法是有效的。例如，用Princeton的气温数据集作为趋势面、再用简单克里金的方法做残差订正，温度的4个区平均CV是3.4℃，而Princeton数据的平均RMSE是4.29℃。然而，BNU数据的平均CV值是1.97℃，明显小于3.4℃。两套方案采用了同样的残差订正方法，BNU产品的优势在于选用了精度更高的趋势面。这说趋势面的质量直接影响最后产品的质量。
最近，中国科学院青藏高原所发展了一套中国区域的大气驱动数据集(ITPCAS,He,2010；Chen et al.,2011)。其中的趋势面选用Princeton大气驱动资料，并用薄板平滑样条插值对残差场进行订正。ITPCAS气温数据的全国平均CV值是2.5℃，仍高于而BNU的1.97℃。其中一个重要原因也是因为BNU采用的质量更好的趋势面。
Hutchinson et al.(2009)认为三维的薄板平滑样条(以经度、纬度和高程作为样条函数的自变量)用于气温场的模拟可能比用高程做协变量的二维样条更好。这里本发明实施例来验证这个结论。用CFSR数据作为协变量的三维薄板平滑样条模型建立的趋势面的CV值在1，2，3，4区分别是：3.01℃,2.24℃,1.92℃和1.40℃。而BNU产品对这四个区域的CV值分别是2.80℃，2.04℃，1.82℃,1.41℃。显然Hutchinson et al.(2009)的结论对前三个区并不正确。虽然从理论上三维样条可以模拟出温度递减率的空间变化率(Hutchinson,2009)，但可能由于前三个区的观测比较稀疏，不能够支持更复杂的三维样条模型的参数估计，反而比较简单的二维样条模型效果更好。4区的观测更密集些，三维样条的CV与二维样条的几乎一样，这也进一步验证了本发明实施例的结论。
Barnes空间插值(Barnes 1964；1978)也是一种常用的方法。为了消除地形对气温的影响，先把气温应用0.65℃/(100m)的温度递减率降到海平面上，然后应用Barnes方法进行插值，再把最后的插值的气温结果升到地形高度上。表2-1列出了Barnes方法的CV值，它在每一个区都大于BNU产品的CV值。可见本发明实施例提出的方法优于传统的Barnes空间插值。
表2-1 Barnes方法和本发明实施例提出的方法的CV比较(℃)
Figure GDA0003469179360000151
(2)与其他格点数据集的比较
Princeton驱动数据集是现在应用于陆面模式最为流行的大气驱动数据，而CFSR再分析资料是最近几年发展的高时空分辨率的再分析资料，所以本发明实施例将BNU驱动资料和这两套驱动资料分别在观测台站和格点上进行比较。
首先是在观测台站上的比较将CFSR和Princeton气温数据插值到观测台站上的RMSE结果。结果表明，BNU气温数据的CV要比CFSR以及Princeton资料的RMSE要小很多，这说明了BNU的气温数据在观测台站尺度上比相应的CFSR和Princeton驱动数据集的精度要高。而且，本发明实施例在气温趋势面模型的基础上进行了残差订正，所以最后的BNU格点驱动数据集在附近有台站观测的格点上的误差要比趋势面的CV还要小。
其次在格点上的比较。验证BNU相对湿度产品在0.3125°×0.3125°(CFSR尺度上)和1°×1°(Princeton尺度上)上的精度，根据(2-19)得出的BNU气温数据的方差比Princeton数据和CFSR数据的方差降低量的估计，见表2-2。从结果来看，这些值都是正的并且都比BNU温度场趋势面的CV要大，说明BNU驱动数据在0.3125°×0.3125°的空间分辨率上要比CFSR气温场精确，并且在1°×1°的空间分辨率上也比Princeton气温场精确。
表2-2 BNU近地表气温数据与Princeton，CFSR相比降低的方差估计的平方根(℃)
Figure GDA0003469179360000161
3.2近地面相对湿度场
3.2.1格点场的建立
对于相对湿度q，在不同时间段采用不同的趋势面模型，在每个有观测的时刻，应用如下趋势面模型。
1958-1978:q(x，y)＝f(x，y)+ε(x，y)
1979-2010:q(x，y)＝f(x，y)+β1qcfsr(x，y)+ε(x，y) (2-12)
符号含义同前，应用1.2.2节的方法对趋势面模型残差进行订正。
其中在1990年之前，利用日观测数据建立模型，1990-1997年，用6小时观测资料建立模型，1998-2010年是用3小时观测资料建立模型。最后对于1990年之前和1990-1997两个时间段的插值结果利用1.2.3节描述的方法进行时间降尺度插值。借助的辅助资料与前面气温场相同。
3.2.2方法评估与精度验证
(1)与其他方法的比较
本小节的研究思路与气温场中“与其他方法的比较”大致相同。
第一个参与对比的方案是：Princeton的相对湿度数据集作为趋势面，然后应用简单克里金的方法对Princeton趋势面的残差进行简单克里金订正。Princeton相对湿度数据的RMSE为21.19％，经过简单克里金订正之后的CV为15.9％，而应用本发明实施例提出的方法算出的CV为10.54％。这说明对Princeton相对湿度数据进行订正确实是有效的。但是由于Princeton数据作为趋势面，其精度远远不如BNU数据趋势面的精度，所以基于Princeton相对湿度数据订正后的结果会比BNU数据差一些。
其次是与ITPCAS驱动数据制作方案做对比，即以Princeton数据集作为背景场，基于台站观测数据并利用不含协变量的薄板平滑样条模型对背景场的残差场进行订正。应用该方法计算得出的区域平均的CV值为12.28％，而本发明实施例提出的方法的平均CV结果为10.54％。产生这个差异的一个主要原因也是趋势面的质量。
最后将本发明实施例提出的方案与Barnes方法进行对比，验证结果见表2-3，可以看出，本发明实施例建立的相对湿度场方案要比Barnes方案精确的多。
表2-3本发明实施例建立相对湿度场的方案与Barnes方案的精度对比(％)
Figure GDA0003469179360000171
(2)与其他格点数据集的比较
将BNU相对湿度驱动资料分别和CFSR、Princeton大气驱动数据集中的相对湿度数据进行站点尺度上和格点尺度上的对比。
首先是在观测台站上的比较，将CFSR和Princeton数据的相对湿度数据插值到中国大陆地区的观测台站上的RMSE。结果表明，BNU相对湿度数据的CV要比CFSR以及Princeton资料的RMSE要小很多，这说明了BNU的相对湿度数据在观测台站尺度上比相应的CFSR和Princeton数据的精度要高。而且，本发明实施例在相对湿度趋势面模型的基础上进行了残差订正，所以最后的BNU格点驱动数据集在附近有台站观测的格点上的误差要比趋势面的CV还要小。
其次在格点上的比较。验证BNU相对湿度产品在0.3125°×0.3125°(CFSR尺度上)和1°×1°(Princeton尺度上)上的精度，BNU相对湿度数据的方差比Princeton数据和CFSR数据的方差降低的估计，见表2-4。从结果来看，这些值都是正的并且都比BNU相对湿度场趋势面的CV要大，说明BNU驱动数据在0.3125°×0.3125°的空间分辨率上要比CFSR相对湿度场精确，并且在1°×1°的空间分辨率上也比Princeton相对湿度场精确。
表2-4 BNU近地表气温数据与CFSR，Princeton相比降低的方差估计的平方根(％)
Figure GDA0003469179360000172
Figure GDA0003469179360000181
3.3近地面风场
3.3.1格点场的建立
在1958-2010年期间，在每个有观测的时刻，采用如下模型：
w(x，y)＝f(x，y)+ε(x，y) (2-13)
对趋势面残差不进行订正，并且应用1.2.3节的方法对日尺度插值结果和6小时尺度的插值结果进行时间降尺度(和气温场完全相同)。
3.3.2方法评估与精度验证
(1)与其他方法的比较
第一种参与对比的方案：应用Princeton的风速数据集作为趋势面，并对Princeton风速趋势面的残差进行简单克里金订正，Princeton数据在台站上的全国平均RMSE的结果为3.19m/s，对Princeton数据的残差进行克里金订正后的结果为1.3m/s，而本发明实施例建立的风场的CV结果是0.51m/s。从中可以看出，虽然基于Princeton数据进行残差订正是有效果的，但是，利用Princeton资料作为趋势面远远不如本发明实施例提出的薄板平滑样条模型作为趋势面的精度高，而且也说明趋势面的质量能直接影响到最后的风速产品的质量。
其次是与ITPCAS驱动数据制作方案做对比：即以Princeton数据集作为背景场，基于台站观测数据并利用不含协变量的薄板平滑样条模型对背景场的残差场进行订正。应用该方法计算得出的区域平均的CV值分别为0.76m/s，而本发明实施例提出的方法的平均CV结果为0.51m/s。产生这种差异的主要原因也是趋势面的质量。
最后将本发明实施例提出的方案与Barnes方法进行对比，结果如表2-5，从结果中可以看出，本发明实施例建立的风速场方案要比Barnes方案精确的多。
表2-5 Barnes方法和BNU中国地区风场的CV比较(m/s)
Figure GDA0003469179360000182
(2)与其他格点数据集的比较
将BNU风场数据与Princeton、CFSR的风场数据分别进行观测台站尺度上和格点尺度上的比较。
首先在观测台站上的比较。将CFSR和Princeton风速数据插值到台站上的结果。结果表明，BNU风速场的CV要比CFSR以及Princeton资料的RMSE要小很多，这说明了BNU风速场的精度在观测台站尺度上比相应的CFSR和Princeton驱动数据集的精度要高。
其次在格点上的比较。验证BNU风速数据集在CFSR产品的分辨率(0.3125°×0.3125°)上和Princeton风速数据的分辨率(1°×1°)上的精度。表2-6为BNU风速数据的方差比相应的CFSR以及Princeton数据的方差降低的平方根。从表中可以看出，这些值都是正的且都比BNU数据趋势面的CV要大，说明BNU风速驱动数据集在CFSR风速数据的空间分辨率上比CFSR数据精确的多，而且在Princeton风速数据的空间分辨率上也比Princeton数据精确的多。
表2-6 BNU近地表风速数据与CFSR数据以及Princeton数据相比降低的方差估计的平方根(m/s)
Figure GDA0003469179360000191
3.4近地面气压场
3.4.1格点场的建立
在1958-2010年期间，在每个有观测的时刻，采用如下模型：
p'(x,y)＝f(x,y)+s'(z(x,y))+ε(x,y) (2-14)
其中s'是一维平滑样条函数，其他的符号的含义见近地面气温场的建立。应用3.3节的方法对日尺度插值结果和6小时尺度的插值结果进行时间降尺度(和气温场完全相同)。
3.4.2方法评估与精度验证
(1)与其他方法的比较
第一种参与比较的方案：应用Princeton的气压数据集作为趋势面，并对Princeton气压趋势面的残差进行订正。Princeton气压数据在台站上的区域平均RMSE为4.04hPa，利用简单克里金方法对Princeton进行订正后的CV结果为3.9hPa，而本发明实施例提出的方法的平均CV结果是1.42hPa。因此可以看出，对Princeton气压数据进行残差订正确实是有效的，但用这种方法做出气压场的CV值要高于本发明实施例采用的方法。原因在于利用Princeton资料作为趋势面，其精度不如本发明实施例提出的薄板平滑样条模型作为趋势面的精度高。
其次是与ITPCAS驱动数据制作方案做对比：以Princeton数据集作为背景场，基于台站观测数据并利用不含协变量的薄板平滑样条模型对背景场的残差场进行订正，应用该方法计算得出的区域平均的CV值为3.06hPa，本发明实施例提出的方法的平均CV结果是1.42hPa，趋势面的质量是产生这种差异的主要原因。
最后将本发明实施例提出的方案与Barnes方法进行对比并进行验证，结果如表2-7，从结果可以看出，本发明实施例建立的气压场方案要比Barnes方案精确的多。
表2-7 Barnes方法和BNU气压场建立方法的CV比较(hPa)
Figure GDA0003469179360000201
(2)与其他格点数据集的比较
将BNU近地面气压数据集与CFSR、Princeton数据集中的气压数据进行站点尺度上和格点尺度上的比较。
首先在观测台站上的比较。将CFSR和Princeton气压数据插值到观测站点上的RMSE结果如表8所示。气压数据的插值过程高程的订正采用Slonosky(2001)提出的方法，由表2-8可知，BNU气压数据的趋势面的CV要比CFSR以及Princeton资料的RMSE要小很多，这说明了BNU气压数据的精度在观测台站尺度上比相应的CFSR和Princeton驱动数据集的精度要高。
表2-8 Princeton、CFSR、和BNU气压数据的精度比较(hPa)
Figure GDA0003469179360000202
其次在格点上的比较。验证BNU气压数据集在CFSR产品的分辨率(0.3125°×0.3125°)上和Princeton气压数据的分辨率(1°×1°)上的精度。表2-9为BNU气压数据的方差比相应的CFSR数据以及Princeton数据的方差降低估计的平方根。从结果来看，BNU驱动数据在CFSR产品的空间分辨率上的精度比CFSR驱动数据的精度要高，而且在Princeton数据的空间分辨率上的精度比Princeton的精度高。
表2-9 BNU气压数据与CFSR、Princeton数据相比降低的方差估计的平方根(hPa)
Figure GDA0003469179360000211
参考两个时间点的气压台站观测数据、BNU数据产品、CFSR再分析资料气压数据和Princeton气压数据的中国大陆空间分布图。总的来说，BNU、CFSR、Princeton的气压格点数据都和台站观测数据比较吻合，这是因为气压变量在很大程度上受高程所决定，所以这三套资料的气压空间分布和观测基本一致。
3.5降水场
3.5.1格点场的建立
在整个1958-2010年用日降水资料建立趋势面，在以下时间段的趋势面模型分别为
1958-1978：r(x，y)＝f(x，y)+ε(x，y)
1979-1997：r(x，y)＝f(x，y)+α·rcfsr(x，y)+ε(x，y)
1998-2010：r(x)＝f(x，y)+α·rcmorph(x，y)+ε(x，y) (2-15)
其中r代表日观测降水量，r(x.y)和rcmorph(x，y)分别为日平均的CFSR再分析资料和CMORPH降水产品在(x，y)处的线性插值，α为计算系数；其它符号与公式(1)相同。类似于气压场的建立，不对降水趋势面的残差场进行订正。
利用1.2.3节描述的方法，对日降水资料向3小时时间分辨率进行时间降尺度插值。其中1958-1978年借助于3小时的Princeton资料作为辅助资料，1979-1997年采用CFSR资料作为辅助资料，1998-2010是采用CMORPH资料作为辅助资料。
3.5.2方法评估与精度验证
结果显示，BNU降水数据集的精度要高于CMORPH资料，说明了卫星观测的精确度不如基于观测数据建立的降水驱动数据集，同时也表明了基于观测台站建立降水驱动数据集的重要意义。
本发明实施例挑选了2003年7月份降水量比较大的四天(8日，23日，26日，27日)。对比了气象台站观测资料、BNU产品、CFSR产品和CMORPH产品。
3.6辐射场
3.6.1格点场的建立
本发明实施例应用台站观测资料日照时数，日最高气温，日最低气温和日平均相对湿度建立下行短波辐射。日照时数是指太阳在一地实际照射地面的时数，即地面观测地点受到太阳直接辐射辐照度等于和大于120W/m2的累计时间。这里需要注意的是，日照时数的观测是统一采用北京时间。
第一步是利用1.2.1节方法建立日照时数在两个时间段的回归模型，如下：
1958-1982和2008-2010：s(x，y)＝f(x，y)+β1td(x，y)+β2F(q(x，y))+ε(x，y)
1983-2007：s(x，y)＝f(x，y)+β1td(x，y)+β2F(q(x，y))+β3G(x，y)+ε(x，y) (2-16)
其中s为日照时数，td为温度在北京时间的日变化(即日最高气温减去日最低气温)，q为在北京时间的日平均相对湿度，
Figure GDA0003469179360000221
G(x，y)是短波辐射遥感产品GEWEX SRB在处的插值。(20)和(21)的拟合方案同温度场(11)和(12)的估计方法，并利用5公里格点的经纬度信息、温度日变化、日均相对湿度以及GEWEX SRB产品生成5公里格点场的日尺度日照时数。
第二步是利用中国科学院青藏高原研究所阳坤教授等所提出的太阳辐射模型计算5公里3小时的下行短波辐射(Yang et al.，2001,2006)。计算步骤为：首先，根据前面已经建立的近地面气温场，相对湿度场和气压场，计算日平均气温、日平均相对湿度和日平均气压(北京时间)，并结合日照时数场，计算日尺度的下行短波辐射。其次，根据太阳高度角的变化，将下行短波辐射的日值分配到3小时的时间尺度上，最后，生成时空分辨率为5公里3小时的下行短波辐射场，并将北京时间转换成世界时间。
下行长波辐射场是采用Crawford(1999)中提出的算法，由短波辐射、近地面气温、相对湿度和气压场直接计算生成。
3.6.2方法评估与精度验证
本发明实施例将使用2003年全国90个下行短波辐射观测站的日观测数据作为“真值”来验证产品的质量。为了验证GEWEX SRB数据对于建立日照时数模型的作用，本发明实施例利用回归模型(20)和(21)分别计算了2003年的日照时数，分别用这两种日照时数的结果计算了下行短波辐射，两者的比较结果见表2-10。从结果来看，用GEWEX SRB数据作为协变量在各个区域都能提高对短波辐射的精度，尤其是在西部地区解释方差提高了31％左右，在整个中国大陆区域平均解释方差提高了20％左右。
同时，为了同Princeton和GEWEX SRB的遥感反演产品的短波辐射精度做对比，表10中也给出了Princeton资料和GEWEX SRB遥感资料的RMSE。从结果可以看出，GEWEX SRB的短波辐射产品精度远好于Princeton再分析资料。表2-10还给出了用和不用GEWEXSRB协变量的短波辐射的CV值。除了东南部(4区)，在中国其它地区用GEWEXSRB作协变量的短波辐射的CV值小于GEWEX SRB的CV值，这就提供了一种更好地应用GEWEX SRB改进短波辐射的新途径。不用GEWEX SRB协变量的短波辐射的CV值(34.44W/m2)略大于GEWEX SRB的RMSE(31.41W/m2)，但明显地小于Princeton短波辐射的RMSE(57.65W/m2)。由此可见在没有GEWEX遥感产品的年份，本发明实施例作出的短波辐射场的精度大抵相当于GEWEX产品的质量，而在有GEWEX数据的年份，可以得到优于GEWEX产品的短波辐射数据。
表2-10应用2003年每月3号的数据对下行短波辐射的精度验证(W/m2)
Figure GDA0003469179360000231
由于下行长波辐射无观测值，因此无法在站点上验证。
1.分类数据的数字化
4.1分类数据的数字化方法
在土地利用和土壤类型中，有一类数据很常见，就是分类数据。例如：土壤类型分粘土、沙土等，植被类型分森林、灌木、农作物、草地等，这些类别的数据往往没有得到很好的利用，现有文献中，人们往往把不同的类型分别来进行分析，这样做的不足就是不同类型里面暗含的“次序”信息没有被很好地利用到，比如，森林的根深比草地深；沙土的粒径比粘土大。若这一种次序的信息能和百分比信息一起整合到本发明实施例的数据里，就能对流域的地表作出更精准的描述。现在本发明实施例定义两种新的数据，一是土壤平均粒径，一是地表植物的平均根深。这两种通过定义秩序来将植被类型数据，土壤分类数据进行尺度统一，例如土壤平均粒径的概念：
Figure GDA0003469179360000232
其中Si是第i种土壤粒径的序数，Ri是这种粒径的土壤所占的比例。地表植被平均根深数据的定义也同公式(23)。这样对于某一具体的流域而言，有两个连续的变量来分别刻画土壤和植被的特征，这样当流域样本量足够多时，本发明实施例可以从数据中提取一些有用的信息，
4.2数字化的分类数据的应用举例
将分类数据数字化之后，本发明实施例将所的到的连续变量与其他流域数据甚至是参数整合在一起，进行多源数据的聚类分析：
聚类分析的结果显示，流域的降水特点和土壤平均粒径以及植被类型，海拔高度有尤为密切的关系，可见本发明实施例要根据上中下游的不同的地理条件选择不同的多源降水数据融合方案以及水文模型的参数化方案。通过几种不同的降水数据和参数化方案的模拟效果进行比较，本发明实施例得出了和子流域属性相关的最优的降水数据融合以及水文模型的参数化方案。
除此之外，同化技术的应用也是一种融合多源数据，充分利用有利信息的有效手段。目前在陆面、水文模型中，集合卡尔曼滤波的应用十分广泛。集合卡尔曼滤波的计算代价较小，计算灵活，适合在地形地貌十分复杂，信息来源比较多样的流域如雅砻江流域进行利用。通过历史数据的训练，本发明实施例可以对于不同的子流域选取适合其特点的最优同化方案。
2.服务于流域水文模型模拟的土壤参数库的构建与评估
以中科院南京土壤所全国土种数据为基础，探讨了流域水文模型建模所需的土壤物理参数与水力学特征参数的估算方法。利用模拟精确度较高的土壤粒径分布模型与土壤传递函数进行土壤粒径的转换与水文参数的估算，进而构建了服务于流域水文模型模拟的土壤参数库。并在普适似然不确定性估计方法(GLUE)框架内，以典型流域水文过程模拟效果评估了土壤参数数据库的有效性。
评估是通过对比两种估算方法下参数在模型模拟中的不确定性实现的。建立两种估算方法下的分布式水文模型：一是土壤参数分三种土壤全部通过率定获取参数值，即19个参数参与率定，设为CAL_19；二是土壤参数通过土壤传递函数计算得到，即除土壤参数外的10个参数参与率定，设为CAL_10。两种估算方法下参数的后验分布与平均分布的先验分布差别均不明显。参与评估的各土壤参数中，大部分参数间相关性较低，个别参数间相关性在0.01显著性水平下显著相关。通过分析两种估算方法下的模型模拟效果，在CAL_10条件下，模型模拟的精度较高，不确定性指标值较低；95％置信区间不确定带较窄；大于临界值的参数组模拟得到的NSE值整体较高，即在CAL_10条件下NSE出现最高概率的值比在CAL_19条件下的值高。通过以上几个角度的分析，使用本土壤参数库参数值可以提高模型模拟精度并降低模型的不确定性，该数据库的建立对雅砻江流域内的水文过程模拟具有较强应用价值。
对于所属技术领域的技术人员而言，随着技术的发展，本发明构思可以不同方式实现。本发明的实施方式并不仅限于以上描述的实施例，而且可在权利要求的范围内进行变化。

Claims (2)

1.一种多数据源的数据融合方法，其特征在于，包括：
基础数据处理步骤，用于对台站气象数据、CFSR气象数据、Princeton大气驱动数据、遥感数据进行处理以获取处理后的基础数据；
数据集建立步骤，用于通过根据处理后的基础数据建立近地面气温场、近地面相对湿度场、近地面风场、近地面气压场、降水场、辐射场；
所述基础数据处理步骤具体包括：
台站气象数据处理子步骤，用于获取来自于中国气象局气象信息中心常规气象要素740站的观测数据，其中所述观测数据中包括以下的至少一种观测变量：近地面1.5米气温、气压、相对湿度、近地面10米风速、累积降水和日照时数；
CFSR气象数据再分析子步骤：用于获取CFSR再分析资料；其中CFSR再分析资料中包括：近地面气温、相对湿度、风速和气压变量的3小时分辨率的数据，以为大气、海洋、陆地和海冰模式提供初始场；
Princeton大气驱动数据再分析子步骤：基于NCEP-1再分析资料，再利用月气候资料对NCEP-1再分析资料进行订正，然后进行空间降尺度；其中Princeton大气驱动数据中至少包括：近地面气温、相对湿度、风速、气压和短波辐射数据；
遥感数据处理子步骤：用于对CMORPH降水资料进行降维处理，并对GWEXE SRB下行短波辐射资料中的地表下行短波辐射资料进行处理；
所述数据集建立步骤中通过以下方法建立近地面气温场：
步骤11、建立趋势面，在每个有观测的时刻，建立如下的趋势面模型，
1958-1978：t(x,y)＝f(x,y)+β·z(x,y)+ε(x,y)
1979-2010：t(x,y)＝f(x,y)+β1·z(x,y)+β2·tcfsr(x,y)+ε(x,y) (2-11)
式中t代表近地面2米气温，z为高程，tcfsr是CFSR再分析数据在点的线性插值；
其中，某一时刻的驱动变量u分解为
Figure FDA0003469179350000011
其中(x,y)是经纬度，{ψr,r＝1,2,...,p}是给定的协变量函数，βr是回归系数，f是2阶薄板平滑样条函数；ε为误差，是独立同分布的；公式中的
Figure FDA0003469179350000021
合称为变量的趋势面；只要样条函数的系数和回归系数确定了，趋势面相对于经纬度的函数关系就确定了；则在任一点的插值只需将该点的经纬度代入趋势面的函数即可获得；在本步骤中驱动变量u为近地面2米气温t；其中2阶薄板平滑样条函数f是针对于在经度{x1,x2,...xn}和纬度{y1,y2,...yn}上的n维观测向量定义的；其表现为
Figure FDA0003469179350000022
和限制条件
Figure FDA0003469179350000023
d1、d2、d3为一阶系数，ci为二阶系数；
由公式(2-1)、公式(2-2)可知，2阶薄板平滑样条函数f的自由度与观测的个数相同；因此在计算样条系数和回归系数时，不但要使趋势面靠近观测，也要使样条尽量平滑以控制趋势面在无观测地区的误差；在实际操作中，这些系数被取为使得以下目标函数的最小的值，
Figure FDA0003469179350000024
其中第一项代表误差的规模，第二项
Figure FDA0003469179350000025
J2(f)代表2阶薄板平滑样条函数f的光滑性；λ是光滑参数，用于平衡误差规模和光滑性；zi为高程；δ是偏微分符号；d为微分符号；
估计样条系数和回归系数的关键技术是对光滑参数λ的估计；因为λ确定后，目标函数(2-2)只是样条系数和回归系数的2次函数，使它在约束条件(2-3)下极小的系数有显示表达；
基于最小交叉验证原则对光滑参数λ进行验证：
步骤12、对趋势面进行订正，具体包括：
用简单克里金方法对趋势面的残差场插值，得到插值结果为
Figure FDA0003469179350000031
其中，||x，y；xi，yi||表示点(x，y)和(xi，yi)之间的欧氏距离，c是一个单调非负函数的协方差函数；然后将残差场的插值叠加到趋势面后就得到最后的插值场
Figure FDA0003469179350000032
其中协方差函数c的估计采用以下方法：将观测集合按两点间的距离分为若干点对的集合，在每一个集合内计算点对的协方差，以构造一个随距离变化的序列；用一维薄板平滑样条拟合这个序列就可得到协方差函数c；
步骤13、将得到的近地面气温场与CFSR气象数据、Princeton大气驱动数据、遥感数据进行对比以确定数据的准确性；
处理子步骤：用于对CMORPH降水资料进行降维处理，并对GWEXE SRB下行短波辐射资料中的地表下行短波辐射资料进行处理；
所述数据集建立步骤中通过以下方法建立近地面相对湿度场：
步骤21、对于相对湿度q，在每个有观测的时刻，建立如下的趋势面模型，
1958-1978：q(x，y)＝f(x，y)+ε(x，y)
1979-2010：q(x，y)＝f(x，y)+β1qcfsr(x，y)+ε(x，y)
(2-12)
其中qcfsr是CFSR再分析数据在点的线性插值，其中2阶薄板平滑样条函数f的计算步骤与步骤11相同；
步骤22、对趋势面进行订正，订正方式与步骤12相同；
步骤23、对日尺度插值结果和6小时尺度的插值结果进行时间降尺度，具体包括：
在每一日的监测数据中确定有监测数据的时刻，并在每两个相邻的监测数据之间进行插值，插值的为相邻时刻的驱动场的平均值；
利用插值后的监测数据对辅助资料进行订正以达到时间降尺度的目的；其中该辅助资料为以下的一种：CFSR再分析资料、Princeton大气驱动资料、GEWEX SRB下行短波辐射资料；
其中利用插值后的监测数据对辅助资料进行订正具体包括：对于风速的气象变量，获取其辅助资料，将3小时的辅助资料除以它的日平均再乘以插值的日资料；
步骤24、将得到的地面相对湿度场与CFSR气象数据、Princeton大气驱动数据、遥感数据进行对比以确定数据的准确性；
所述数据集建立步骤中通过以下方法建立近地面风场：
步骤31、对于近地面风场w，在每个有观测的时刻，采用如下模型：
w(x，y)＝f(x，y)+ε(x，y) (2-13)
步骤32、对日尺度插值结果和6小时尺度的插值结果进行时间降尺度，具体包括：
在每一日的监测数据中确定有监测数据的时刻，并在每两个相邻的监测数据之间进行插值，插值的为相邻时刻的驱动场的平均值；
利用插值后的监测数据对辅助资料进行订正以达到时间降尺度的目的；其中该辅助资料为以下的一种：CFSR再分析资料、Princeton大气驱动资料、GEWEX SRB下行短波辐射资料；
其中利用插值后的监测数据对辅助资料进行订正具体包括：对于风速的气象变量，获取其辅助资料，将3小时的辅助资料减去它的日平均再加上插值的日资料；
所述数据集建立步骤中通过以下方法建立近地面气压场：
步骤41、对于近地面气压场p′，在每个有观测的时刻，采用如下模型：
p′(x，y)＝f(x，y)+s′(z(x，y))+ε(x，y) (2-14)
其中s′是一维平滑样条函数；
步骤42、对日尺度插值结果和6小时尺度的插值结果进行时间降尺度，具体包括：
在每一日的监测数据中确定有监测数据的时刻，并在每两个相邻的监测数据之间进行插值，插值为相邻时刻的驱动场的平均值；
利用插值后的监测数据对辅助资料进行订正以达到时间降尺度的目的；其中该辅助资料为以下的一种：CFSR再分析资料、Princeton大气驱动资料、GEWEX SRB下行短波辐射资料；
其中利用插值后的监测数据对辅助资料进行订正具体包括：对于气压的气象变量，获取其辅助资料，将3小时的辅助资料减去它的日平均再加上插值的日资料；
所述数据集建立步骤中通过以下方法建立降水场：
步骤51、在整个1958-2010年用日降水资料建立趋势面，在以下时间段的趋势面模型分别为
1958-1978：r′(x，y)＝f(x，y)+ε(x，y)
1979-1997：r′(x，y)＝f(x，y)+α·rcfsr(x，y)+ε(x，y)
1998-2010：r′(x，y)＝f(x，y)+α·rcmorph(x，y)+ε(x，y) (2-15)
其中r’代表日观测降水量，rcfsr(x.y)和rcmorph(x，y)分别为日平均的CFSR再分析资料和CMORPH降水产品在(x，y)处的线性插值；α为计算系数；
步骤52、对日尺度插值结果和6小时尺度的插值结果进行时间降尺度，具体包括：
在每一日的监测数据中确定有监测数据的时刻，并在每两个相邻的监测数据之间进行插值，插值的为相邻时刻的驱动场的平均值；利用插值后的监测数据对辅助资料进行订正以达到时间降尺度的目的；其中该辅助资料为以下的一种：CFSR再分析资料、Princeton大气驱动资料、GEWEX SRB下行短波辐射资料；其中利用插值后的监测数据对辅助资料进行订正具体包括：对于风速的气象变量，获取其辅助资料，将3小时的辅助资料除以它的日平均再乘以插值的日资料；
所述数据集建立步骤中通过以下方法建立辐射场：利用台站观测资料中的日照时数、日最高气温、日最低气温、日平均相对湿度，建立下行短波辐射；日照时数是指太阳在一地实际照射地面的时数，即地面观测地点受到太阳直接辐射辐照度等于和大于120W/m2的累计时间；
具体包括：
步骤61、建立日照时数在两个时间段的回归模型：
1958-1982和2008-2010：s(x，y)＝f(x，y)+β1td(x，y)+β2F(q(x，y))+ε(x，y)
1983-2007：s(x，y)＝f(x，y)+β1td(x，y)+β2F(q(x，y))+β3G(x，y)+ε(x，y) (2-16)
其中s为日照时数，td为温度在北京时间的日变化，q为在北京时间的日平均相对湿度；
Figure FDA0003469179350000061
G(x,y)是短波辐射遥感产品GEWEX SRB在处的插值；利用5公里格点的经纬度信息、温度日变化、日均相对湿度以及GEWEX SRB产品生成5公里格点场的日尺度日照时数；
步骤62、利用太阳辐射模型计算5公里3小时的下行短波辐射；具体包括：
根据前面已经建立的近地面气温场，相对湿度场和气压场，计算日平均气温、日平均相对湿度和日平均气压，并结合日照时数场，计算日尺度的下行短波辐射；
根据太阳高度角的变化，将下行短波辐射的日值分配到3小时的时间尺度上，最后，生成时空分辨率为5公里3小时的下行短波辐射场，并将北京时间转换成世界时间；
其中下行长波辐射场是由短波辐射、近地面气温、相对湿度和气压场直接计算生成。
2.根据权利要求1所述的多数据源的数据融合方法，其特征在于，所述方法还包括：
分类数据的数字化处理步骤，用于对土壤数据进行分类处理；具体包括：
土壤平均粒径计算步骤，根据以下公式计算土壤平均粒径：
Figure FDA0003469179350000062
其中Si’是第i’种土壤粒径的序数，Ri’是这种粒径的土壤所占的比例；
地表植物的平均根深计算步骤，利用地表植物的序数和这种地表植物所占的比例，计算地表植物的平均根深；
数据聚类分析步骤，用于将数字化的土壤和地表植物的数字化参数与其他流域数据/参数整合在一起，进行多源数据的聚类分析，以获取流域的降水特点与土壤平均粒径、植被类型、海拔高度之间的关系；
根据上中下游的不同的地理条件选择不同的多源降水数据融合方案以及水文模型的参数化方案；通过不同的降水数据和参数化方案的模拟效果进行比较，以获取与子流域属性相关的最优的降水数据融合以及水文模型的参数化方案。
CN201810974820.9A 2018-08-24 2018-08-24 一种多数据源的数据融合方法 Active CN109344865B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810974820.9A CN109344865B (zh) 2018-08-24 2018-08-24 一种多数据源的数据融合方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810974820.9A CN109344865B (zh) 2018-08-24 2018-08-24 一种多数据源的数据融合方法

Publications (2)

Publication Number Publication Date
CN109344865A CN109344865A (zh) 2019-02-15
CN109344865B true CN109344865B (zh) 2022-03-04

Family

ID=65291830

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810974820.9A Active CN109344865B (zh) 2018-08-24 2018-08-24 一种多数据源的数据融合方法

Country Status (1)

Country Link
CN (1) CN109344865B (zh)

Families Citing this family (11)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109886354B (zh) * 2019-03-06 2020-12-22 国家卫星海洋应用中心 风场融合方法及装置
CN111948736A (zh) * 2019-05-14 2020-11-17 中国电力科学研究院有限公司 一种基于大数据平台的高维天气预报数据降维方法
CN110609923A (zh) * 2019-07-31 2019-12-24 象辑知源（武汉）科技有限公司 一种分布式的多算法融合的气象数据插值方法
CN111126466B (zh) * 2019-12-16 2020-09-22 西安科技大学 一种多源pwv数据融合方法
CN111369483B (zh) * 2020-03-05 2020-11-13 北京师范大学 一种融合多源遥感数据生成高时空分辨率遥感数据的方法
CN112069449B (zh) * 2020-09-04 2021-07-16 中科三清科技有限公司 基于初值集合的天气预报方法及装置
CN112579725B (zh) * 2021-01-04 2022-11-15 四创科技有限公司 一种河湖多源异构数据聚合关联方法及系统
CN112946368A (zh) * 2021-02-02 2021-06-11 长春工程学院 一种多源参数融合分析的杆塔接地电阻同步检测方法
CN114330146B (zh) * 2022-03-02 2022-06-28 北京英视睿达科技股份有限公司 一种卫星气体数据补全方法和系统
CN115292554B (zh) * 2022-09-27 2023-01-17 中国民用航空局空中交通管理局航空气象中心 一种航空气象四维数据集的构建方法及系统
CN116722544B (zh) * 2023-08-02 2023-10-20 北京弘象科技有限公司 分布式光伏短期预测方法、装置、电子设备及存储介质

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103426026A (zh) * 2013-09-10 2013-12-04 信阳师范学院 一种混合神经网络预测及识别景区气象要素的方法
CN103488871A (zh) * 2013-08-27 2014-01-01 国家电网公司 一种无径流资料地区的洪水预报方法
WO2016076394A1 (ja) * 2014-11-14 2016-05-19 国立研究開発法人海洋研究開発機構 画像処理装置、画像処理方法、画像処理プログラム

Family Cites Families (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105868529A (zh) * 2016-03-18 2016-08-17 北京师范大学 一种基于遥感的近地表日均大气温度反演方法
CN108416031B (zh) * 2018-03-12 2021-07-27 南京恩瑞特实业有限公司 气象多源探测资料融合分析系统

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103488871A (zh) * 2013-08-27 2014-01-01 国家电网公司 一种无径流资料地区的洪水预报方法
CN103426026A (zh) * 2013-09-10 2013-12-04 信阳师范学院 一种混合神经网络预测及识别景区气象要素的方法
WO2016076394A1 (ja) * 2014-11-14 2016-05-19 国立研究開発法人海洋研究開発機構 画像処理装置、画像処理方法、画像処理プログラム

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
LI Tao et al..Mapping Near-surface Air Temperature, Pressure, Relative Humidity and Wind Speed over Mainland China with High Spatiotemporal Resolution.《ADVANCES IN ATMOSPHERIC SCIENCES,》.2014, *
Mapping Near-surface Air Temperature, Pressure, Relative Humidity and Wind Speed over Mainland China with High Spatiotemporal Resolution;LI Tao et al.;《ADVANCES IN ATMOSPHERIC SCIENCES,》;20140930;第1-8页 *
MAPPING RAINFALL FIELDS AND THEIR ENSO VARIATION IN DATA-SPARSE TROPICAL SOUTH-WEST PACIFIC OCEAN REGION;REID E. BASHER et al.;《INTERNATIONAL JOURNAL OF CLIMATOLOGY》;19981231;第239-242页 *
Thin-Plate Smoothing Spline Modeling of Spatial CVmate Data and Its Application to Mapping South Pacifac Rainfalls;xiaogu zheng etal;《MONTHLY WEATHER REVIEW》;19951031;第3087-3091页 *
一种综合反映土地利用变化对流域土壤侵蚀及泥沙量反馈影响的建模方法及应用;方青青等;《北京师范大学学报（ 自然科学版）》;20141231;第483-490页 *

Also Published As

Publication number Publication date
CN109344865A (zh) 2019-02-15

Similar Documents

Publication Publication Date Title
CN109344865B (zh) 一种多数据源的数据融合方法
Wang et al. Evaluation of the GPM IMERG satellite-based precipitation products and the hydrological utility
Charbit et al. Assessment of the Greenland ice sheet–atmosphere feedbacks for the next century with a regional atmospheric model coupled to an ice sheet model
CN102539336B (zh) 基于环境一号卫星的可吸入颗粒物估算方法及系统
CN105674901B (zh) 基于统计学分段的大气气溶胶反演方法
CN112785024B (zh) 一种基于流域水文模型的径流计算和预测方法
CN105912836A (zh) 一种纯遥感数据驱动的流域水循环模拟方法
CN114265836A (zh) 一种基于云区温湿度廓线反演的卫星微波温湿度计全天候同化方法及装置
CN110909447B (zh) 一种高精度电离层区域短期预报方法
CN108983324A (zh) 一种基于卡尔曼滤波的气温预报方法及系统
Zheng et al. Improved forecast skill through the assimilation of dropsonde observations from the Atmospheric River Reconnaissance program
ZHANG et al. Developing a process-based and remote sensing driven crop yield model for maize (PRYM–Maize) and its validation over the Northeast China Plain
CN115420690A (zh) 近地表痕量气体浓度反演模型及反演方法
Yin et al. Evaluation of ORCHIDEE-MICT-simulated soil moisture over China and impacts of different atmospheric forcing data
Hong et al. 15 Global Precipitation Estimation and Applications
Wright Rainfall information for global flood modeling
CN116050163A (zh) 一种基于气象站的生态系统碳水通量计算方法及系统
CN112836449B (zh) 一种用于率定水文模型的方法
Guo et al. Mapping the spatial distribution and time evolution of snow water equivalent with passive microwave measurements
Wu et al. Triple collocation-based estimation of spatially correlated observation error covariance in remote sensing soil moisture data assimilation
El Badaoui et al. The prediction of moisture through the use of neural networks MLP type
Ekici et al. Total Global Solar Radiation Estimation with Relative Humidity and Air Temperature Extremes in Ireland and Holland
De Blasi Scale dependence of hydrological effects from different climatic conditions on glacierized catchments
Hao et al. Evaluation of snow processes over the Western United States in E3SM land model
Li et al. Spatiotemporal estimation of model error to improve soil moisture analysis in ensemble Kalman filter data assimilation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CP01 Change in the name or title of a patent holder
CP01 Change in the name or title of a patent holder

Address after: No. 3377, Jingshi Road, Licheng District, Jinan City, Shandong Province 250010

Patentee after: Shandong Institute of ecological environment planning

Patentee after: Beijing Normal University

Address before: No. 3377, Jingshi Road, Licheng District, Jinan City, Shandong Province 250010

Patentee before: SHANDONG ACADEMY FOR ENVIRONMENTAL PLANNING

Patentee before: Beijing Normal University