CN113935956A - 一种二向混合建模矿区土壤含水量数据缺失修复方法 - Google Patents

一种二向混合建模矿区土壤含水量数据缺失修复方法 Download PDF

Info

Publication number
CN113935956A
CN113935956A CN202111116905.1A CN202111116905A CN113935956A CN 113935956 A CN113935956 A CN 113935956A CN 202111116905 A CN202111116905 A CN 202111116905A CN 113935956 A CN113935956 A CN 113935956A
Authority
CN
China
Prior art keywords
pixel
missing
water content
soil
model
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.)
Granted
Application number
CN202111116905.1A
Other languages
English (en)
Other versions
CN113935956B (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.)
Beijing Shulun Technology Co ltd
China University of Mining and Technology Beijing CUMTB
China Energy Investment Corp Ltd
National Institute of Clean and Low Carbon Energy
Shenhua Beidian Shengli Energy Co Ltd
Original Assignee
Beijing Shulun Technology Co ltd
China University of Mining and Technology Beijing CUMTB
China Energy Investment Corp Ltd
National Institute of Clean and Low Carbon Energy
Shenhua Beidian Shengli Energy Co Ltd
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 Beijing Shulun Technology Co ltd, China University of Mining and Technology Beijing CUMTB, China Energy Investment Corp Ltd, National Institute of Clean and Low Carbon Energy, Shenhua Beidian Shengli Energy Co Ltd filed Critical Beijing Shulun Technology Co ltd
Priority to CN202111116905.1A priority Critical patent/CN113935956B/zh
Publication of CN113935956A publication Critical patent/CN113935956A/zh
Application granted granted Critical
Publication of CN113935956B publication Critical patent/CN113935956B/zh
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/0002—Inspection of images, e.g. flaw detection
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q10/00—Administration; Management
    • G06Q10/04—Forecasting or optimisation specially adapted for administrative or management purposes, e.g. linear programming or "cutting stock problem"
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/02—Agriculture; Fishing; Forestry; Mining
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T7/00—Image analysis
    • G06T7/60—Analysis of geometric attributes
    • G06T7/62—Analysis of geometric attributes of area, perimeter, diameter or volume
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T2207/00—Indexing scheme for image analysis or image enhancement
    • G06T2207/30—Subject of image; Context of image processing
    • G06T2207/30204—Marker

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Business, Economics & Management (AREA)
  • Mathematical Physics (AREA)
  • Data Mining & Analysis (AREA)
  • Mathematical Analysis (AREA)
  • Pure & Applied Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Computational Mathematics (AREA)
  • Economics (AREA)
  • Human Resources & Organizations (AREA)
  • Strategic Management (AREA)
  • Marketing (AREA)
  • Operations Research (AREA)
  • General Business, Economics & Management (AREA)
  • Algebra (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Tourism & Hospitality (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Quality & Reliability (AREA)
  • Game Theory and Decision Science (AREA)
  • Primary Health Care (AREA)
  • Probability & Statistics with Applications (AREA)
  • Animal Husbandry (AREA)
  • Health & Medical Sciences (AREA)
  • Agronomy & Crop Science (AREA)
  • General Health & Medical Sciences (AREA)
  • Computing Systems (AREA)
  • Evolutionary Biology (AREA)
  • Geophysics And Detection Of Objects (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Geometry (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Development Economics (AREA)
  • Mining & Mineral Resources (AREA)

Abstract

本发明公开了一种二向混合建模矿区土壤含水量数据缺失修复方法，其方法如下：A、将整幅矿区土壤含水量影像像元点分别为已知像元和缺失像元，B、根据缺失像元Bj所对应的驱动因子集信息通过预测模型得到缺失像元的预测土壤含水量
Figure DDA0003275660570000011
C、采用地理加权回归模型GWR构建已知像元的关联模型，D、用已知像元的实际土壤含水量SMAi减去模型拟合土壤含水量
Figure DDA0003275660570000012
并得到误差εi再做普通克里金插值计算得到潜在误差εBj，E、根据缺失像元的预测土壤含水量
Figure DDA0003275660570000013
并结合潜在误差εBj对缺失像元按照如下公式进行缺失误差修正。本发明能够提升单独用正向建模预测缺失像元的精度，尤其是在较大面积的缺失情况下正向建模无法覆盖缺失区域，能够高精度修复较大面积成块缺失数据。

Description

一种二向混合建模矿区土壤含水量数据缺失修复方法
技术领域
本发明涉及生态学领域、煤炭领域、遥感及地理信息领域，尤其涉及一种二向混合建模矿区土壤含水量数据缺失修复方法。
背景技术
遥感传感器在成像过程中，不可避免地会受到云、雪、气溶胶和传感器自身性能等因素的影响，使得遥感数据在空间上呈现不连续，即存在缺失数据的情况。矿区场景研究区域范围小(通常仅为几公里到数十公里)、地表空间异质性强，要求空间分辨率高，数据缺失会严重阻碍空间异质性的分析。所以在反演矿区场景下的土壤含水量时，若不能修补缺失数据将严重影响反演结果的实用价值与成果图件的视觉效果。现有关于数据插补方法主要分为2大类：第一类方法是采用滤波法根据时间序列数据估算缺失数据，此类方法依赖于时间序列遥感数据，对于时间分辨率要求低的产品数据修复具有较高的稳定性和准确性，但对于土壤含水量此类要求时间分辨率高的产品来说不适用。第二类方法是采用空间插值法根据邻近数据估算缺失数据，例如反距离权重法(Inverse Distance Weighting，IDW，参见如下文献：张优,王娟,张杰,等.GIS与地统计学的土壤水分空间插值方法[J].四川师范大学学报(自然科学版),2019, 42(05):703-710)、普通克里金插值法(Ordinary Kriging，OK，参见如下文献：张优,王娟,张杰,等.GIS与地统计学的土壤水分空间插值方法[J].四川师范大学学报(自然科学版),2019,42(05):703-710)、径向基函数插值法 (Radial BasisFunction，RBF，参见如下文献：周银明,吴达胜.基于RBF 和IDW的气象要素插值方法比较[J].计算机时代,2019(10):8-10)和回归克里格法(Return Kriging，RK，参见如下文献：姚雪玲,傅伯杰,吕一河,等.基于GIS和统计模型的黄土丘陵沟壑区土壤水分插值方法[J].水土保持学报, 2013,27(006):93-96.)等，这些方法受到区域样点均匀性的影响，对土壤含水量这种具有空间非平稳性的参数也不适用。地理加权回归克里格 (GeographicallyWeighted Regression Kriging，GWRK，参见如下文献：张国峰,杨立荣,瞿明凯,等.基于地理加权回归克里格的日平均气温插值[J].应用生态学报,2015,26(05):1531-1536.)；虽然考虑到空间非平稳性，对已知像元局部建模，根据空间相关性来预测其周边缺失像元处的土壤含水量，但由于土壤含水量空间异质性强、变化剧烈，受到模型最优带宽的限制，已知像元处的模型仅对其带宽以内的像元适用，对带宽以外的缺失像元不具备适用性。若增加带宽，模型适用范围会增大，但对该带宽范围内的缺失像元来说并不是最优带宽模型，失去最优性。综上所述，GWRK虽然解决了空间非平稳性问题，但对于遥感影像反演中较大面积成块缺失土壤含水量数据无法精确预测，无法同时满足模型的大面积适用性与最优性。
发明内容
针对现有技术存在的不足之处，本发明的目的在于提供一种二向混合建模矿区土壤含水量数据缺失修复方法，在正向构建已知像元关联模型的基础上，利用缺失像元周边的已知像元信息反向构建缺失像元处的预测模型，通过二向混合建模的方法充分考虑缺失像元模型的适用性与最优性，弥补直接应用周边已知像元的GWR模型来预测缺失像元土壤含水量的不足。
本发明的目的通过下述技术方案实现：
一种二向混合建模矿区土壤含水量数据缺失修复方法，其方法如下：
A、得到矿区土壤含水量影像，整幅矿区土壤含水量影像像元点分为已知像元Ai和缺失像元Bj，遍历整幅影像像元点，将有土壤含水量信息的像元点标记为已知像元Ai，保存在List A中，将没有土壤含水量信息的像元点标记为缺失像元Bj，保存在List B中；优选地，本发明中的矿区土壤含水量影像为栅格影像。
B、在矿区土壤含水量影像上建立一个滑动窗口，滑动窗口的带宽长度为h0，采用滑动窗口法并利用滑动窗口对保存在ListB中的所有缺失像元Bj进行逐像元点遍历，在List A中搜索已知像元信息，已知像元信息为已知像元Ai所对应的信息，信息包括驱动因子集信息与土壤含水量信息，采用地理加权回归模型 GWR构建土壤含水量与驱动因子之间的定量关系模型，将该定量关系模型作为缺失像元Bj的预测模型；同时在List B中得到缺失像元Bj所对应的驱动因子集信息，根据缺失像元Bj所对应的驱动因子集信息通过预测模型得到缺失像元Bj的预测土壤含水量
Figure BDA0003275660550000031
C、采用滑动窗口法并利用滑动窗口对保存在List A中的所有已知像元逐像元遍历，并在List A中搜索滑动窗口范围内的所有已知像元信息，已知像元信息包括驱动因子集信息与土壤含水量信息，采用地理加权回归模型GWR构建土壤含水量与驱动因子之间的定量关系模型，将该定量关系模型作为已知像元 Ai的关联模型；
D、潜在误差计算：首先，将List A中的所有已知像元信息中的驱动因子集信息代入关联模型计算得到已知像元的模型拟合土壤含水量
Figure BDA0003275660550000032
然后用已知像元的实际土壤含水量SMAi减去模型拟合土壤含水量
Figure BDA0003275660550000033
并得到误差εi，
Figure BDA0003275660550000034
再对所有已知像元的误差εi做普通克里金插值计算得到缺失像元处的潜在误差εBj；
E、根据缺失像元Bj的预测土壤含水量
Figure BDA0003275660550000035
并结合潜在误差εBj对缺失像元 Bj按照如下公式进行缺失误差修正：缺失像元Bj的土壤含水量
Figure BDA0003275660550000036
优选地，本发明步骤B中预测模型优选采用反向构建模型，预测模型中滑动窗口的带宽长度方法如下：
B1、在List B第一个像元上建立一个滑动窗口，滑动窗口的初始带宽长度为h1，利用AICc信息来确定最优带宽并作为滑动窗口的带宽长度，对在滑动窗口带宽内的所有ListA像元做回归，可得到对于带宽长度为h的AICc，其公式如下：
Figure BDA0003275660550000041
其中，n表示观测点个数，
Figure BDA0003275660550000042
表示随机误差项方差的估计值，S表示帽子矩阵，tr(S)表示帽子矩阵S 的迹，帽子矩阵是观测到的y到拟合值
Figure BDA0003275660550000043
的投影矩阵，其中对于地理加权回归模型GWR，这个帽子矩阵的每一行ri为：
ri＝Xi(XTWiX)-1XTWi
其中Xi为自变量X矩阵的第i行，Wi为核函数矩阵，XT为自变量X的矩阵转置。
优选地，本发明步骤B1之后还包括步骤B2；
B2、改变带宽，重新选择List A中的样本参与回归，按照如下公式确定最优带宽并作为滑动窗口的带宽长度h0：
最优带宽公式：
Figure BDA0003275660550000044
优选地，本发明步骤B或/和步骤C中采用地理加权回归模型GWR构建土壤含水量与驱动因子之间的定量关系模型，其表达式如下：
Figure BDA0003275660550000045
式中，(xi1,xi2,…,xim；yi) 表示第i个像元(ui,vi)处的土壤含水量yi和驱动因子xi1,xi2,…,xim的观测值；βk(ui,vi),(k＝0,1,…,m)是第i个像元的第k个回归系数；(ε1,ε2,…,εn)为独立同分布的误差项，通常假定εi～N(0,σ2)，Cov(εi,εj)＝0(i≠j)。根据本发明的一个优选实施例，本发明优选的回归系数采用加权最小二乘法，对任一像元i的回归系数表达式如下：
Figure BDA0003275660550000046
式中，X为驱动因子抽样矩阵，其第一列取值为1，以估计截距项β0(ui,vi)；y为土壤含水量抽样值列向量；
Figure BDA0003275660550000051
为像元(ui,vi)处的回归分析系数向量；W(ui,vi)为空间核函数矩阵，其对角线元素值为每个数据点到像元(ui,vi)的空间权重值：
Figure BDA0003275660550000052
式中，对角线值wij(j＝1,2,…,n)表示像元点j对像元点i影响的权重值。
优选地，本发明步骤B或/和步骤C中采用地理加权回归模型GWR中核函数为：
Figure BDA0003275660550000053
式中，dij表示像元位置i与像元位置j之间的空间距离，dN表示像元位置i到最近邻第N个点的空间距离，h为带宽值，单位均为 m。
优选地，本发明步骤B或/和步骤C或/和步骤D中驱动因子集信息包括气候气象因子、地理因子与人类活动因子，气候气象因子包括降水、气温、蒸发量、太阳辐射，地理因子包括坡度、坡向、海拔，人类活动因子包括城镇开发、煤炭开采。本发明优选的人类活动因子中城镇开发可以采用如下公式表达：
Figure BDA0003275660550000054
式中，Urbanm为像元m处的城镇开发影响大小，n为像元个数，S为城镇面积，dm-Urban表示像元m 到城镇边界的最短欧氏距离。
优选地，本发明优选的人类活动因子中煤炭开采表达公式如下：
Figure BDA0003275660550000055
式中，Minem为像元m处的煤炭开采影响大小，n为像元个数，C为煤炭开采量，dm-Mine表示像元m到矿坑边界的最短欧氏距离。
本发明较现有技术相比，具有以下优点及有益效果：
(1)本发明在正向构建已知像元关联模型的基础上，利用缺失像元周边的已知像元信息反向构建缺失像元处的预测模型，通过二向混合建模的方法充分考虑缺失像元模型的适用性与最优性，弥补直接应用周边已知像元的GWR模型来预测缺失像元土壤含水量的不足。
(2)本发明能够高精度修复较大面积成块缺失数据，本发明通过正反双向建模，能够提升单独用正向建模预测缺失像元的精度，尤其是在较大面积的缺失情况下正向建模无法覆盖缺失区域，模型丧失精确性。
附图说明
图1为实施例一中GWRK与二向混合建模示意图；
图2为实施例一中的实验样区图；
图3为反向构建缺失像元预测模型的原理示意图；
图4为本发明实施例二的流程示意图；
图5为实施例二中举例的城镇开发影响因子空间化示意图；
图6为实施例二中举例的煤炭开采影响因子空间化结果示意图；
图7为实施例二中正向GWR回归模型展示图；
图8为实施例二中潜在误差展示图；
图9为实施例二中二向混合建模方法修复结果图。
具体实施方式
下面结合实施例对本发明作进一步地详细说明：
实施例一
一种二向混合建模矿区土壤含水量数据缺失修复方法，其方法如下：
A、得到矿区土壤含水量影像，整幅矿区土壤含水量影像像元点分为已知像元Ai和缺失像元Bj，遍历整幅影像像元点，将有土壤含水量信息的像元点标记为已知像元Ai，保存在List A中，将没有土壤含水量信息的像元点标记为缺失像元Bj，保存在List B中；优选地，本实施例的矿区土壤含水量影像为栅格影像。
B、在矿区土壤含水量影像上建立一个滑动窗口，滑动窗口的带宽长度为h0，采用滑动窗口法并利用滑动窗口对保存在ListB中的所有缺失像元Bj进行逐像元点遍历，在List A中搜索已知像元信息，已知像元信息为已知像元Ai所对应的信息，信息包括驱动因子集信息与土壤含水量信息，采用地理加权回归模型 GWR构建土壤含水量与驱动因子之间的定量关系模型，将该定量关系模型作为缺失像元Bj的预测模型(预测模型是基于已知像元Ai进行反向构建的预测模型结构)；同时在List B中得到缺失像元Bj所对应的驱动因子集信息，根据缺失像元Bj所对应的驱动因子集信息通过预测模型得到缺失像元Bj的预测土壤含水量
Figure BDA0003275660550000071
为了筛选出本实施例应用矿区场景下和土壤含水量相关性最高的驱动因子 (驱动因子集信息所包括的驱动因子)，本实施例可对所有驱动因子和土壤含水量做多元逐步回归，剔除驱动因子中不太重要又和其他因子高度相关的因子，降低多重共线性程度。经过本发明研究团队不断筛选与试验，得到最优的驱动因子集信息包括：气候气象因子、地理因子和人类活动因子，具体如下：气候气象因子包括降水、气温、蒸发量、太阳辐射，地理因子包括坡度、坡向、海拔，人类活动因子包括城镇开发、煤炭开采。
优选地，步骤B中预测模型采用反向构建模型(参见图3)，预测模型中滑动窗口的带宽长度方法如下：
B1、在List B第一个像元上建立一个滑动窗口，滑动窗口的初始带宽长度为h1，利用AICc信息(或称AICc信息准则法，Corrected Akaike information criterion，简称AICc)来确定最优带宽并作为滑动窗口的带宽长度，对在滑动窗口带宽内的所有List A像元做回归，可得到对于带宽长度为h的AICc，其公式如下：
Figure BDA0003275660550000081
其中，n表示观测点个数，
Figure BDA0003275660550000082
表示随机误差项方差的估计值，S表示帽子矩阵，tr(S)表示帽子矩阵S 的迹，帽子矩阵是观测到的y到拟合值
Figure BDA0003275660550000083
的投影矩阵，其中对于地理加权回归模型GWR，这个帽子矩阵的每一行ri为：
ri＝Xi(XTWiX)-1XTWi
其中Xi为自变量X矩阵的第i行，Wi为核函数矩阵，XT为自变量X的矩阵转置。
在步骤B1之后还包括步骤B2；
B2、改变带宽，重新选择List A中的样本参与回归(即重新选择样本参与回归)，按照如下公式确定最优带宽并作为滑动窗口的带宽长度h0：
最优带宽公式：
Figure BDA0003275660550000084
C、采用滑动窗口法并利用滑动窗口对保存在List A中的所有已知像元逐像元遍历，并在List A中搜索滑动窗口范围内的所有已知像元信息，已知像元信息包括驱动因子集信息与土壤含水量信息，采用地理加权回归模型GWR构建土壤含水量与驱动因子之间的定量关系模型，将该定量关系模型作为已知像元 Ai的关联模型(关联模型是基于已知像元Ai进行正向构建而得)；
D、潜在误差计算：首先，将List A中的所有已知像元信息中的驱动因子集信息代入关联模型计算得到已知像元的模型拟合土壤含水量
Figure BDA0003275660550000085
然后用已知像元的实际土壤含水量SMAi减去模型拟合土壤含水量
Figure BDA0003275660550000086
并得到误差εi，
Figure BDA0003275660550000087
(已知像元的实际值减去拟合值得到误差εi)；再对所有已知像元的误差εi做普通克里金插值计算得到缺失像元处的潜在误差εBj。
E、根据缺失像元Bj的预测土壤含水量
Figure BDA0003275660550000088
并结合潜在误差εBj对缺失像元 Bj按照如下公式进行缺失误差修正：缺失像元Bj的土壤含水量
Figure BDA0003275660550000089
由此得到缺失像元Bj修正后的矿区土壤含水量数据。
根据本发明的一个实施例，本发明步骤B或/和步骤C中采用地理加权回归模型GWR构建土壤含水量与驱动因子之间的定量关系模型，其表达式如下：
Figure BDA0003275660550000091
式中，(xi1,xi2,…,xim；yi) 表示第i个像元(ui,vi)处的土壤含水量yi和驱动因子xi1,xi2,…,xim的观测值；βk(ui,vi),(k＝0,1,…,m)是第i个像元的第k个回归系数；(ε1,ε2,…,εn)为独立同分布的误差项，通常假定εi～N(0,σ2)，Cov(εi,εj)＝0(i≠j)。回归系数采用加权最小二乘法，对任一像元i的回归系数表达式如下：
Figure BDA0003275660550000092
式中，X为驱动因子抽样矩阵，其第一列取值为1，以估计截距项β0(ui,vi)；y为土壤含水量抽样值列向量；
Figure BDA0003275660550000093
为像元(ui,vi)处的回归分析系数向量；W(ui,vi)为空间核函数矩阵，其对角线元素值为每个数据点到像元(ui,vi)的空间权重值：
Figure BDA0003275660550000094
式中，对角线值wij(j＝1,2,…,n)表示像元点j对像元点i影响的权重值。
根据本发明的一个优选实施例，本发明步骤B或/和步骤C中采用地理加权回归模型GWR中核函数为：
Figure BDA0003275660550000095
式中，dij表示像元位置i与像元位置j之间的空间距离，dN表示像元位置i到最近邻第N个点的空间距离， h为带宽值，单位均为m。
根据本发明的一个优选实施例，本发明步骤B或/和步骤C或/和步骤D中驱动因子集信息包括气候气象因子、地理因子与人类活动因子，气候气象因子包括降水、气温、蒸发量、太阳辐射，地理因子包括坡度、坡向、海拔，人类活动因子包括城镇开发、煤炭开采。其中，人类活动因子中城镇开发可以采用如下表达式：
Figure BDA0003275660550000101
式中，Urbanm为像元m处的城镇开发影响大小，n为像元个数，S为城镇面积，dm-Urban表示像元m 到城镇边界的最短欧氏距离。
其中，人类活动因子中煤炭开采可以采用如下表达式：
Figure BDA0003275660550000102
式中，Minem为像元m处的煤炭开采影响大小，n为像元个数，C为煤炭开采量，dm-Mine表示像元m到矿坑边界的最短欧氏距离。
本实施例提供如下表格和数据来论述本发明在矿区土壤含水量数据缺失修复方法的有效性，采用实验对比了RK、GWRK、本申请二向混合建模方法的优缺点、适用性与修复精度。
Figure BDA0003275660550000103
表1三种方法优缺点对比
对于如图1所示大面积成块缺失区域(参见图1缺失像元边界内的缺失像元)，现要预测缺失区域中的A像元，GWRK方法直接应用离该点最近的已知像元B的关联模型，但该已知像元点B的关联模型最优带宽为h无法覆盖到缺失像元A，如图1a所示(GWRK模型最优性分析)，此类情况将导致缺失像元A的预测模型不具备适用性，预测精度降低。若扩大带宽，如图1b所示(GWRK模型适用性分析)，使模型带宽h1覆盖缺失像元A，则已知像元B处的关联模型用于预缺失像元A是适用的，但失去最优性。因此对于此类大面积缺失区域，GWRK 无法准确预测。
本发明提出一种二向混合建模的矿区土壤含水量数据缺失修复方法，如图1c 所示(本申请二向混合建模方法对应示例)，直接在缺失像元A处建立回归模型，将最优带宽h0内的已知像元纳入回归模型，使缺失像元A的预测模型同时具备适用性与最优性。
为了定量分析二向混合建模方法用于土壤含水量缺失数据修复的精度及与其他方法的对比效果，本实施例选取锡林浩特市胜利矿区附近一小区域为实验样区(参见图2)，以2004年完整的土壤含水量数据作为参考，人工模拟成块缺失数据，对该区域分别应用RK、GWRK、二向混合建模方法进行修复，并与原始土壤含水量进行对比，比较修复精度。
以土壤含水量原始数据为参照，计算三种方法用于修复缺失数据的平均相对误差(Mean Relative Estimation Error，MRE)，公式如下：
Figure BDA0003275660550000111
式中，n为样本点的个数，yi为第i个样本点的真实值，
Figure BDA0003275660550000112
为第i个样本点的预测值。
Figure BDA0003275660550000113
表2三种方法修复误差对比
由表2可以看出，回归克里格(RK)方法修复缺失数据结果的平均相对误差为0.008，地理加权回归克里格方法修复结果平均相对误差为0.0041，误差降低48.75％。二向混合建模方法相对GWRK方法修复结果误差降低9.76％。与RK、 GWRK这两种方法相比，二向混合建模的方法修复缺失数据的精度更高，虽然插补结果与实际数值仍然有一定偏差，但在成片数据缺失且不依赖长时序序列数据的情况下，能有效反映处数据缺失区域的整体状况，精度达到99.63％。
本发明在正向构建已知像元关联模型的基础上，利用缺失像元周边的已知像元信息反向构建缺失像元处的预测模型，通过二向混合建模的方法充分考虑缺失像元模型的适用性与最优性，弥补直接应用周边已知像元的GWR模型来预测缺失像元土壤含水量的不足。本发明能够高精度修复较大面积成块缺失数据，本发明通过正反双向建模，能够提升单独用正向建模预测缺失像元的精度，尤其是在较大面积的缺失情况下正向建模无法覆盖缺失区域，模型丧失精确性。
实施例二
如图4所示，一种二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：其方法如下：
A、缺失像元识别：对缺失像元进行插补前需要识别缺失像元的空间位置。具体如下：得到矿区土壤含水量影像，将整幅矿区土壤含水量影像分为缺失像元和已知像元。对整幅矿区土壤含水量影像像元点进行遍历，将有土壤含水量信息的像元点标记为已知像元Ai，保存在List A中；将值为Nodata的像元点标记为缺失像元Bj，保存在List B中。
B、反向构建缺失像元的预测模型：采用滑动窗口法，首先在整幅矿区土壤含水量影像(栅格影像)上建立一个滑动窗口，滑动窗口的带宽长度为h0，采用AICc信息准则法(Corrected Akaike information criterion，AICc)来确定最优带宽。滑动窗口会对ListB中的缺失像元进行逐像元遍历，并在List A 中搜索滑动窗口范围内的所有已知像元信息，包括驱动因子集信息和土壤含水量信息。采用地理加权回归模型(GWR)构建土壤含水量与各驱动因子之间的定量关系模型作为该缺失像元的预测模型。
B1、为了筛选出该矿区场景下和土壤含水量相关性最高的驱动因子，对所有驱动因子和土壤含水量做多元逐步回归，剔除驱动因子中不太重要又和其他因子高度相关的因子，降低多重共线性程度。经过本发明研究团队不断筛选与试验，得到驱动因子集包括：气候气象因子、地理因子和人类活动因子，具体如下：气候气象因子包括降水、气温、蒸发量、太阳辐射，地理因子包括坡度、坡向、海拔，人类活动因子包括城镇开发、煤炭开采。
步骤B中采用地理加权回归模型建立土壤含水量与驱动因子之间的定量关系模型，表达式如下：
Figure BDA0003275660550000131
式中，(xi1,xi2,…,xim；yi)表示第i 个缺失像元(ui,vi)处的土壤含水量yi和驱动因子xi1,xi2,…,xim的观测值；
βk(ui,vi),(k＝0,1,…,m)是第i个缺失像元的第k个回归系数；(ε1,ε2,…,εn)为独立同分布的误差项，通常假定εi～N(0,σ2)，Cov(εi,εj)＝0(i≠j)。
回归系数估算方法为加权最小二乘估计，对任一缺失像元i的估计系数表达式如下：
Figure BDA0003275660550000132
式中，X为驱动因子抽样矩阵，其第一列取值为1，以截距项β0(ui,vi)；y为土壤含水量抽样值列向量；
Figure BDA0003275660550000133
为缺失像元(ui,vi)处的回归分析系数向量；W(ui,vi)为空间核函数矩阵，其对角线元素值为每个数据点到缺失像元(ui,vi)的空间权重值：
Figure BDA0003275660550000141
式中，对角线值wij(j＝1,2,…,n)表示第j个已知像元到缺失像元i的权重值。本发明优选核函数为：
Figure BDA0003275660550000142
式中，dij表示像元位置i与像元位置j之间的空间距离，dN表示像元位置i到最近邻第N个点的空间距离，h为带宽值，单位均为m。
确定最优带宽h的方法为AICc信息准则法，其数学表达式如下：
Figure BDA0003275660550000143
其中，n表示观测点个数，S表示帽子矩阵， tr(S)表示帽子矩阵S的迹，
Figure BDA0003275660550000144
表示随机误差项方差的估计值，
Figure BDA0003275660550000145
最优带宽h0的选取如下：
Figure BDA0003275660550000146
地理加权回归模型输入模型的长时序气候气象因子、地理因子和人类活动因子如表3所示：
Figure BDA0003275660550000147
表3输入模型数据
步骤B1中人类活动因子(城镇开发、煤炭开采)空间化过程如下：
首先通过土地利用分类数据栅格计算器得到该缺失年份的城镇面积S和城镇、矿区边界。
1、城镇开发因子空间化过程：本发明以GDP来表征城镇开发强度，在某一时期，城镇开发强度为一定值，但由于城镇开发对土壤含水量的影响具有距离衰减性，本发明定义某一像元上城镇开发影响因子计算公式如下：
Figure BDA0003275660550000151
式中，Urbanm为像元m处的城镇开发影响大小，n为像元个数，S为城镇面积，dm-Urban表示像元m到城镇边界的最短欧氏距离(参见图5的城镇开发影响因子空间化示意图，以锡林浩特市胜利矿区2004年为例)。
2、煤炭开采因子空间化过程：本发明以煤炭开采量C来表征煤炭开采强度，在某一时期，煤炭开采强度为一定值C，但由于煤炭开采对土壤含水量的影响具有距离衰减性，本发明定义某一像元上煤炭开采影响因子计算公式如下：
Figure BDA0003275660550000152
式中，Minem为像元m处的煤炭开采影响大小，n为像元个数，C为煤炭开采量，dm-Mine表示像元m到矿坑边界的最短欧氏距离(参见图6的煤炭开采影响因子空间化结果示意图，以锡林浩特市胜利矿区2004年为例)。
C、求缺失像元处的土壤含水量模型拟合值：将所有缺失像元的气候气象因子、地理因子和人类活动因子代入步骤B中得到的预测模型，得到每个缺失像元的预测土壤含水量
Figure BDA0003275660550000153
D、构建已知像元的关联模型(参见图7的正向GWR回归模型展示图，以部分像元点为例)：对保存在List A中的所有已知像元逐像元遍历，并在List A 中搜索滑动窗口范围内的所有已知像元信息，同样采用地理加权回归模型(GWR) 构建土壤含水量与各驱动因子之间的定量关系模型作为该已知像元的关联模型。已知像元的关联模型如图8所示，SM为土壤含水量，Tp、T2m、Eva、DEM、 M、U分别为降水、气温、蒸发量、采矿活动、城镇开发。
E、潜在误差计算：首先将所有已知像元的气候气象因子、地理因子和人类活动因子代入步骤D中得到的关联模型，得到每个已知像元的模型拟合土壤含水量
Figure BDA0003275660550000161
然后用已知像元的实际值减去拟合值得到误差
Figure BDA0003275660550000162
其中SMAi为已知像元实际土壤含水量。最后对所有已知像元的误差εi做普通克里金插值得到缺失像元处的潜在误差εBj(参见图8的潜在误差展示图)。
F、缺失像元插补：将步骤E中计算得到的潜在误差加到步骤C模型拟合得到的缺失像元土壤含水量值上，弥补确定性分量中的潜在误差。最后缺失像元土壤含水量为
Figure BDA0003275660550000163
(参见图9的二向混合建模方法修复结果图)。
以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：其方法如下：
A、得到矿区土壤含水量影像，整幅矿区土壤含水量影像像元点分为已知像元Ai和缺失像元Bj，遍历整幅影像像元点，将有土壤含水量信息的像元点标记为已知像元Ai，保存在List A中，将没有土壤含水量信息的像元点标记为缺失像元Bj，保存在List B中；
B、在矿区土壤含水量影像上建立一个滑动窗口，滑动窗口的带宽长度为h0，采用滑动窗口法并利用滑动窗口对保存在ListB中的所有缺失像元Bj进行逐像元点遍历，在List A中搜索已知像元信息，已知像元信息为已知像元Ai所对应的信息，信息包括驱动因子集信息与土壤含水量信息，采用地理加权回归模型GWR构建土壤含水量与驱动因子之间的定量关系模型，将该定量关系模型作为缺失像元Bj的预测模型；同时在List B中得到缺失像元Bj所对应的驱动因子集信息，根据缺失像元Bj所对应的驱动因子集信息通过预测模型得到缺失像元Bj的预测土壤含水量
Figure FDA0003275660540000016
C、采用滑动窗口法并利用滑动窗口对保存在List A中的所有已知像元逐像元遍历，并在List A中搜索滑动窗口范围内的所有已知像元信息，已知像元信息包括驱动因子集信息与土壤含水量信息，采用地理加权回归模型GWR构建土壤含水量与驱动因子之间的定量关系模型，将该定量关系模型作为已知像元Ai的关联模型；
D、潜在误差计算：首先，将List A中的所有已知像元信息中的驱动因子集信息代入关联模型计算得到已知像元的模型拟合土壤含水量
Figure FDA0003275660540000011
然后用已知像元的实际土壤含水量SMAi减去模型拟合土壤含水量
Figure FDA0003275660540000012
并得到误差εi，
Figure FDA0003275660540000013
再对所有已知像元的误差εi做普通克里金插值计算得到缺失像元处的潜在误差εBj；
E、根据缺失像元Bj的预测土壤含水量
Figure FDA0003275660540000014
并结合潜在误差εBj对缺失像元Bj按照如下公式进行缺失误差修正：缺失像元Bj的土壤含水量
Figure FDA0003275660540000015
2.按照权利要求1所述的二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：步骤B中预测模型采用反向构建模型，预测模型中滑动窗口的带宽长度方法如下：
B1、在List B第一个像元上建立一个滑动窗口，滑动窗口的初始带宽长度为h1，利用AICc信息来确定最优带宽并作为滑动窗口的带宽长度，对在滑动窗口带宽内的所有List A像元做回归，可得到对于带宽长度为h的AICc，其公式如下：
Figure FDA0003275660540000021
其中，n表示观测点个数，
Figure FDA0003275660540000022
表示随机误差项方差的估计值，S表示帽子矩阵，tr(S)表示帽子矩阵S的迹，帽子矩阵是观测到的y到拟合值
Figure FDA0003275660540000023
的投影矩阵，其中对于地理加权回归模型GWR，这个帽子矩阵的每一行ri为：
ri＝Xi(XTWiX)-1XTWi
其中Xi为自变量X矩阵的第i行，Wi为核函数矩阵，XT为自变量X的矩阵转置。
3.按照权利要求2所述的二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：步骤B1之后还包括步骤B2；
B2、改变带宽，重新选择ListA中的样本参与回归，按照如下公式确定最优带宽并作为滑动窗口的带宽长度h0：
最优带宽公式：
Figure FDA0003275660540000024
4.按照权利要求1所述的二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：步骤B或/和步骤C中采用地理加权回归模型GWR构建土壤含水量与驱动因子之间的定量关系模型，其表达式如下：
Figure FDA0003275660540000031
式中，(xi1,xi2,…,xim；yi)表示第i个像元(ui,vi)处的土壤含水量yi和驱动因子xi1,xi2,…,xim的观测值；
βk(ui,vi),(k＝0,1,…,m)是第i个像元的第k个回归系数；(ε1,ε2,…,εn)为独立同分布的误差项，通常假定εi～N(0,σ2)，Cov(εi,εj)＝0(i≠j)。
5.按照权利要求4所述的二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：所述回归系数采用加权最小二乘法，对任一像元i的回归系数表达式如下：
Figure FDA0003275660540000032
式中，X为驱动因子抽样矩阵，其第一列取值为1，以估计截距项β0(ui,vi)；y为土壤含水量抽样值列向量；
Figure FDA0003275660540000033
为像元(ui,vi)处的回归分析系数向量；W(ui,vi)为空间核函数矩阵，其对角线元素值为每个数据点到像元(ui,vi)的空间权重值：
Figure FDA0003275660540000034
式中，对角线值wij(j＝1,2,…,n)表示像元点j对像元点i影响的权重值。
6.按照权利要求4所述的二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：步骤B或/和步骤C中采用地理加权回归模型GWR中核函数为：
Figure FDA0003275660540000035
式中，dij表示像元位置i与像元位置j之间的空间距离，dN表示像元位置i到最近邻第N个点的空间距离，h为带宽值，单位均为m。
7.按照权利要求1所述的二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：步骤B或/和步骤C或/和步骤D中驱动因子集信息包括气候气象因子、地理因子与人类活动因子，气候气象因子包括降水、气温、蒸发量、太阳辐射，地理因子包括坡度、坡向、海拔，人类活动因子包括城镇开发、煤炭开采。
8.按照权利要求1所述的二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：步骤A中的矿区土壤含水量影像为栅格影像。
9.按照权利要求7所述的二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：人类活动因子中城镇开发表达公式如下：
Figure FDA0003275660540000041
式中，Urbanm为像元m处的城镇开发影响大小，n为像元个数，S为城镇面积，dm-Urban表示像元m到城镇边界的最短欧氏距离。
10.按照权利要求7所述的二向混合建模矿区土壤含水量数据缺失修复方法，其特征在于：人类活动因子中煤炭开采表达公式如下：
Figure FDA0003275660540000042
式中，Minem为像元m处的煤炭开采影响大小，n为像元个数，C为煤炭开采量，dm-Mine表示像元m到矿坑边界的最短欧氏距离。
CN202111116905.1A 2021-09-23 2021-09-23 一种二向混合建模矿区土壤含水量数据缺失修复方法 Expired - Fee Related CN113935956B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202111116905.1A CN113935956B (zh) 2021-09-23 2021-09-23 一种二向混合建模矿区土壤含水量数据缺失修复方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202111116905.1A CN113935956B (zh) 2021-09-23 2021-09-23 一种二向混合建模矿区土壤含水量数据缺失修复方法

Publications (2)

Publication Number Publication Date
CN113935956A true CN113935956A (zh) 2022-01-14
CN113935956B CN113935956B (zh) 2022-03-25

Family

ID=79276692

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202111116905.1A Expired - Fee Related CN113935956B (zh) 2021-09-23 2021-09-23 一种二向混合建模矿区土壤含水量数据缺失修复方法

Country Status (1)

Country Link
CN (1) CN113935956B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114723619A (zh) * 2022-02-25 2022-07-08 中国科学院空天信息创新研究院 地表辐射产品修复方法、装置、电子设备及存储介质
CN114817221A (zh) * 2022-05-07 2022-07-29 中国长江三峡集团有限公司 一种双源蒸发数据治理提升方法、系统及存储介质
CN115088596A (zh) * 2022-07-05 2022-09-23 江苏经贸职业技术学院 一种基于物联网的农田灌溉系统及其控制方法

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106504210A (zh) * 2016-10-28 2017-03-15 国网四川省电力公司电力科学研究院 一种modis影像数据缺失修复方法
CN106933776A (zh) * 2017-03-02 2017-07-07 宁波大学 一种modis aod产品缺失数据修复的方法
CN108268735A (zh) * 2018-01-29 2018-07-10 浙江大学 基于多源遥感卫星融合数据的地表土壤水分降尺度方法
CN109902259A (zh) * 2019-02-25 2019-06-18 中国科学院地理科学与资源研究所 一种轻量级的缺失时空数据的重构方法
WO2021000361A1 (zh) * 2019-07-04 2021-01-07 浙江大学 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN106504210A (zh) * 2016-10-28 2017-03-15 国网四川省电力公司电力科学研究院 一种modis影像数据缺失修复方法
CN106933776A (zh) * 2017-03-02 2017-07-07 宁波大学 一种modis aod产品缺失数据修复的方法
CN108268735A (zh) * 2018-01-29 2018-07-10 浙江大学 基于多源遥感卫星融合数据的地表土壤水分降尺度方法
CN109902259A (zh) * 2019-02-25 2019-06-18 中国科学院地理科学与资源研究所 一种轻量级的缺失时空数据的重构方法
WO2021000361A1 (zh) * 2019-07-04 2021-01-07 浙江大学 一种基于经验正交函数分解法的静止海洋水色卫星数据重构方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
梁莹: "地下水位监测值缺失修复混合模型研究", 《中国优秀博硕士学位论文全文数据库（硕士）基础科学辑》 *

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN114723619A (zh) * 2022-02-25 2022-07-08 中国科学院空天信息创新研究院 地表辐射产品修复方法、装置、电子设备及存储介质
CN114817221A (zh) * 2022-05-07 2022-07-29 中国长江三峡集团有限公司 一种双源蒸发数据治理提升方法、系统及存储介质
CN114817221B (zh) * 2022-05-07 2023-06-02 中国长江三峡集团有限公司 一种双源蒸发数据治理提升方法、系统及存储介质
CN115088596A (zh) * 2022-07-05 2022-09-23 江苏经贸职业技术学院 一种基于物联网的农田灌溉系统及其控制方法

Also Published As

Publication number Publication date
CN113935956B (zh) 2022-03-25

Similar Documents

Publication Publication Date Title
CN113935956B (zh) 一种二向混合建模矿区土壤含水量数据缺失修复方法
Marquez et al. Hybrid solar forecasting method uses satellite imaging and ground telemetry as inputs to ANNs
Tveito et al. A GIS-based agro-ecological decision system based on gridded climatology
Köhler et al. Critical weather situations for renewable energies–Part B: Low stratus risk for solar power
Kirstetter et al. Toward a framework for systematic error modeling of spaceborne precipitation radar with NOAA/NSSL ground radar–based National Mosaic QPE
Ahmed et al. Statistical downscaling and bias correction of climate model outputs for climate change impact assessment in the US northeast
Kolassa et al. Soil moisture retrieval from AMSR-E and ASCAT microwave observation synergy. Part 1: Satellite data analysis
Yin et al. Comparison of physical and data-driven models to forecast groundwater level changes with the inclusion of GRACE–A case study over the state of Victoria, Australia
CN109297605B (zh) 一种基于中红外与热红外数据的地表温度反演方法
Hazenberg et al. Radar rainfall estimation of stratiform winter precipitation in the Belgian Ardennes
Pulkkinen et al. Nowcasting of convective rainfall using volumetric radar observations
Furtado et al. Cloud microphysical factors affecting simulations of deep convection during the presummer rainy season in Southern China
Pellarin et al. Using spaceborne surface soil moisture to constrain satellite precipitation estimates over West Africa
Gruber et al. The potential of 2D Kalman filtering for soil moisture data assimilation
Jin et al. Dust emission inversion using himawari‐8 AODs over east Asia: an extreme dust event in may 2017
Vargas Godoy et al. The global water cycle budget: A chronological review
Taylor et al. Satellite retrieval of aerosol microphysical and optical parameters using neural networks: a new methodology applied to the Sahara desert dust peak
Bai et al. Global synthesis of two-decade of research on improving PM2. 5 estimation models: From remote sensing and data science perspectives
Coopman et al. Analysis of the thermodynamic phase transition of tracked convective clouds based on geostationary satellite observations
Alerskans et al. Exploring machine learning techniques to retrieve sea surface temperatures from passive microwave measurements
Xi et al. Comparison of marine boundary layer cloud properties from CERES‐MODIS edition 4 and DOE ARM AMF measurements at the Azores
Visweshwaran et al. Sensitivity‐Based Soil Moisture Assimilation for Improved Streamflow Forecast Using a Novel Forward Sensitivity Method (FSM) Approach
Luo et al. DASSO: a data assimilation system for the Southern Ocean that utilizes both sea-ice concentration and thickness observations
Lee et al. Estimation of maximum daily fresh snow accumulation using an artificial neural network model
Colaninno et al. Towards an operational model for estimating day and night instantaneous near-surface air temperature for urban heat island studies: outline and assessment

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
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20220325

CF01 Termination of patent right due to non-payment of annual fee