An Extended Kriging Method to Interpolate Near-Surface Soil Moisture Data Measured by Wireless Sensor Networks

Zhang, Jialin; Li, Xiuhong; Yang, Rongjin; Liu, Qiang; Zhao, Long; Dou, Baocheng

doi:10.3390/s17061390

开放式访问第条

无线传感器网络近地表土壤水分数据的扩展克里格插值方法

¹

北京师范大学全球变化与地球系统科学学院全球变化研究中心，北京100875

²

北京师范大学与中国科学院遥感与数字地球研究所联合主办的遥感科学国家重点实验室，北京师范大学，北京100101

^三

中国环境科学研究院，北京100012

^*

信件应寄给的作者。

传感器 2017,17(6), 1390;https://doi.org/10.3390/s17061390

收到的提交文件：2017年2月28日/修订日期：2017年6月9日/接受日期：2017年6月10日/发布日期：2017年6月15日

（本文属于特刊传感器和农业土地系统的智能传感)

下载

浏览地物

版本注释

摘要

:

在插值无线传感器网络（WSN）网格测量的近地表土壤水分的实践中，传统的带辅助变量的克里格方法，如Co-Kriging和带外部漂移的克里格（KED），由于土壤水分的异质性及其与辅助变量的相关性很低，因此不能取得令人满意的结果。本研究开发了一种扩展克里格方法，用于借助遥感图像进行插值。其基本思想是通过引入谱变量来扩展传统的克里格方法，并在空间和谱组合空间上进行操作。该算法已应用于HiWATER活动中WSN测量的土壤水分数据，以生成2012年6月10日至7月15日的每日地图。为了进行比较，采用了三种传统的克里格方法：普通克里格法（OK），它只使用WSN数据，协同克里格法和KED，两者都将遥感数据作为协变量进行综合。目视检查表明，与OK、Co-Kriging和KED相比，Extended Kriging的结果显示了更多的空间细节。在四种插值结果中，扩展克里格法的均方根误差（RMSE）最小。这表明该方法在结合遥感信息和地面测量数据进行土壤水分插值方面具有优势。

关键词：

无线传感器网络;克里格插值;土壤湿度;光谱变量

1.简介

无线传感器网络（WSN）是一种新的地面数据采集技术，目前需求量很大。它已应用于水文、土壤环境、大气环境、森林气象、火灾等多个领域[1]. WSN是一种跨学科技术，集成了传感器、自动控制、通信和数据分析[2]. 与传统的地面数据观测方法相比，由WSN组成的多个站点有助于测量参数的区域分布。此外，它的体积小，价格相对较低[1]允许在研究区域安装更多节点以满足研究需求。这有助于我们获得更广泛、更密集的特定地表参数观测数据。与卫星遥感技术相比，无线传感器网络是一个更加灵活的平台，可以支持更多种类的传感器[2,三]. 例如，他们可以测量每个观测节点的叶面积指数（LAI）、土壤湿度、土壤温度、电导率和大气红外辐射温度[4,5]. 与遥感不同，遥感数据是在立交桥期间获得的，无线传感器网络的监测是连续和实时的。这补充了遥感测量的面积数据，也为卫星和航空遥感的验证提供了参考数据[6].

然而，无线传感器网络仅测量点数据，而在环境研究和管理中，为了做出有效的决策或合理的解释，越来越需要特定区域的空间连续数据[7]. 因此，为了确定陆面参数的空间分布，必须进行数据插值。一些传统的插值方法，如逆距离加权法（IDW）、多项式法和克里金法，已被广泛用于将点数据转换为空间分布。IDW方法和多项式方法相对快速且易于计算，但不能探索节点之间关系的空间变化，并且估计精度有限[8]. 相比之下，克里金方法使用变异函数分析估计空间变异结构，并考虑空间自相关[9,10]. 克里金方法是一组随机插值方法，包括简单克里金方法、普通克里金（OK）方法、通用克里金方法，协同克里金方法和回归克里金方法以及残差克里金方法[11]. 这些方法已用于各种研究。例如，Maruyama等人[12]使用简单克里格法估计峰值地面速度。Wu等人[13]利用带有纬度、经度和海拔的输入变量的残差克里金法估计美国的月平均温度。Liang等人[14]使用钴克里格法估算每日NO_三-在每日流量的帮助下，农业河流中的N荷载。Cantet[15]使用带外部漂移的克里格（KED）方法绘制了位于小安的列斯群岛的一个小岛的月平均降雨量。通过地质统计插值方法的发展，通过对目标参数的统计或与其他参数的相关性，探索了有关目标参数分布规律或模式的更多信息。

土壤湿度是农业、生态和水文研究中的一个重要因素[16,17,18,19]. 农田土壤水分动态监测是精准农业的重要组成部分。尽管最近的插值算法取得了进展，但近地表土壤水分的空间分布估计仍然不令人满意，因为土壤水分在空间和时间尺度上都是高度不均匀的，即使在短距离内也是如此[20]. 许多研究人员已应用遥感技术反演土壤水分。例如，Fan等人[21]改进了基于紧凑型机载光谱成像仪（CASI）/热机载光谱成像机（TASI）数据的Ts/VI空间反演土壤水分的性能。El Hajj等人[22]使用神经网络（NNs）从灌溉草地区域的X波段SAR数据估计地表土壤水分。Ponnurangam等人[23]使用紧凑的极化分解和表面成分反演来估计裸露和植被覆盖农业土壤的土壤水分。然而，由于光学遥感对土壤水分不直接敏感，而热红外和被动微波遥感通常具有较低的空间分辨率，因此土壤水分在定量遥感反演中存在很大困难。虽然雷达具有提供高空间分辨率的能力，但几十米左右的雷达对地表粗糙度、地形特征和植被更为敏感，这意味着土壤水分反演过程非常困难。地面观测数据，如WSN测量数据，一直是遥感反演土壤水分的必要补充数据。相反，可以将遥感信息视为辅助数据，以帮助插值地面测量的土壤湿度。在这方面，Kang等人[24]借助多源遥感信息，利用回归克里金方法对异质农田上WSN测得的土壤水分进行了放大。Liao等人[25]分析了各种不确定性来源，如土壤特性和地形指数，同时借助协克里金法在两个典型山坡上估算近地表土壤水分含量。然而，由于土壤水分受到许多因素的影响，例如海拔、植被、温度和灌溉，因此很难找出一个显著线性相关的因素来辅助插值。一些研究人员还尝试将数据同化方法用于土壤水分估算。Gao等人[26]基于WSN数据和ASTER（Terra）地表温度测量的辅助信息，使用贝叶斯最大熵（BME）方法估计土壤水分的空间格局。

本文通过在插值算法中引入高分辨率遥感图像光谱变量，对克里格方法进行了扩展。将该算法与OK、Co-kriging和KED算法进行了分析和比较，然后将其应用于WSN在HiWATER活动中获取的土壤水分数据，生成30 m空间分辨率和日时间分辨率的土壤水分时间序列图。

2.材料和方法

2.1. 研究区域

研究区域，如所示图1位于中国西北部黑河流域中游张掖绿洲。HRB是第二大内河流域，以大面积高寒干旱景观和小部分绿洲农业用地为特征。潜在蒸发量每年从1200毫米到1800毫米不等，人工绿洲的年平均降雨量约为177毫米。灌溉是作物的主要水源。

HRB长期以来一直是综合流域研究和水文试验的试验台[27]. HiWATER是一项正在进行的流域生态水文综合实验，始于2012年[28]. 为了提高综合观测能力，作为HiWATER基础实验的一部分，安装了生态水文WSN。2012年5月至2012年9月，在中游主要绿洲5.5km×5.5km重点试验区安装了50个WATERNET节点。在本研究中，我们选择了一个约4.5 km×5.0 km的区域作为实验区，该区域由48个WATERNET节点组成(图1). 所有节点都安装在玉米地里。

2.2. 数据资源

本研究使用了WSN数据和遥感数据。我们使用了2012年6月10日至2012年7月15日期间的数据，当时研究区域的植被覆盖从稀疏变为稠密。由于该地区的主要作物是玉米，其大部分根系分布在0~15cm区域[29]由于对灌溉的依赖性很大，因此需要对灌溉管理的研究进行近地表水分监测。在研究区域部署了48个WATERNET节点，测量两层（4cm、10cm）的土壤水分和土壤温度，每5分钟收集一次数据[30]. 在本研究中，我们使用WSN测量的10 cm深度的土壤水分进行插值。土壤湿度测量基于使用Hydro Probe II（HP-II）传感器的频域反射法[26]. Jin等人提供了WATERNET设计的详细信息和其他信息[31]. 由于无线数据传输过程中，一些WSN节点采集的数据受到传感器噪声和其他异常情况的影响，因此需要进行平滑和降噪处理。此外，当土壤水分接近饱和时，土壤水分传感器无法正常工作，有时会给出异常值[32]. 因此，这部分数据也需要删除。我们首先通过将零值、负值和异常高值（土壤含水量>50%）指定为无效数据NaN来排除异常数据。然后，我们对全天的数据进行平均，并将平均结果用作每个节点的最终土壤水分值。

遥感数据作为辅助信息，用于估算研究区近地表土壤水分。中国HJ微型卫星（HJ-CCD）上的CCD摄像机在四个光谱带（蓝色、绿色、红色、NIR）上对研究区域进行遥感观测，空间分辨率为30 m，重访频率约为2天。对HJ-CCD图像进行了几何和辐射预处理。以1TP Landsat5 TM影像为参考，利用自动匹配的控制点进行了大地校正，控制点的残差小于0.5像素。利用HJ-CCD图像元数据中的系数进行校准，并利用HiWATER战役期间部署在研究区域的CIMEL CE318太阳光度计获取的每日气溶胶和水汽参数，通过6s大气模型进行大气校正。

本研究使用了两个遥感变量NDVI（归一化植被指数）和反照率。它们来自预处理的HJ-CCD图像，如下所示：

N个 天 V（V） 我 = \frac{ρ_{4} - ρ_{三}}{ρ_{4} + ρ_{三}}

(1)

一 我 b条 e（电子） o个 = - 0.5878 ρ_{1} + 0.6543 ρ_{2} + 0.4183 ρ_{三} + 0.38914 ρ_{4}

(2)

哪里

ρ_{1}

到

ρ_{4}

分别表示四个光谱带中的表面反射率，方程（2）的系数由Gao等人获得[33].

由于云层的影响，在这段时间内，只能使用以下日期的五幅清晰的图像（如所示图2)：2012年6月15日、6月19日、9月29日、7月8日和7月13日。由于插值时间为2012年6月10日至2012年7月15日，我们对NDVI和反照率进行了每日线性插值，以进行克里格计算。

气象数据是在Daman Superstation采集的[34,35]它靠近研究区的中心。图3显示了研究区域的降水量和日平均气温。所用灌溉数据由兰州寒冷干旱地区科学数据中心提供(http://westdc.westgis.ac.cn). 尽管记录不完整（存在一些遗漏），但可用零件仍足以满足本研究的需求。从研究区选取6个典型灌区，将其土壤水分估算结果与降水量和灌溉数据进行比较。这些灌溉场的名称和位置如所示图1：金城六号（JC-6）、石桥一号（SQ-1）、商头闸一号（STZ-1）、五星四号（WX-4）、小满二号（XM-2）和中华六号（ZH-6）。

在本研究中，我们引入了土壤湿度图作为计算半变差函数和一些插值分析的模拟数据。该土壤湿度图来源于CASI/TASI数据的机载传感器[21]空间分辨率为3m。

2.3. 方法

本文提出的新的空间插值方法是基于普通克里格算法的。该方法将传统的X和Y空间坐标扩展到空间和光谱组合坐标，并在插值算法中利用遥感衍生光谱变量NDVI和反照率作为补充信息。在本节中，首先给出了传统克里格算法的概述，然后解释了扩展克里格算法所采用的技术。在拟合土壤湿度的半变异函数时，由于WSN节点的数量不足以收集稳健的统计数据，我们使用了中提到的遥感衍生土壤湿度图第2.2节计算方差函数。

2.3.1. 传统克里格方法

（a）基本公式

克里格法是一种源自区域化变量理论的插值方法，它继承了地质统计学的概念[11]. 它被用于在未采样位置提供线性无偏预测，并取决于变量的空间变化在半变异函数中的表达[36,37]. 该方法量化并减少了估计的不确定性，最大限度地减少了自我估计的预测误差[26]. 克里金的核心是一个最优线性无偏估计量，可以表示为[38]:

{Z轴}^{*} ({v（v）}_{0}) = \sum_{我 = 1}^{n个} λ_{我} Z轴 ({v（v）}_{我})

(3)

哪里

{Z轴}^{*}

是位置处变量的估计值

{v（v）}_{0}

,n个是用于插值的最近相邻采样数据点的数量，

λ_{我}

是分配给每个观测值的克里格权重

{Z轴 (v（v）}_{我})

.

最佳估计要求误差方差最小：

σ_{κ}^{2} = 变量 [Z轴 ({v（v）}_{0}) - {Z轴}^{*} ({v（v）}_{0})] = E类 {{[Z轴 ({v（v）}_{0}) - \sum_{我 = 1}^{n个} λ_{我} Z轴 ({v（v）}_{我})]}^{2}} = 最小值

(4)

为了确保无偏估计，以下约束必须满足以下等式：

\sum_{我 = 1}^{n个} λ_{我} = 1

(5)

为了解决这个约束优化问题，采用了拉格朗日乘子法（LMM）。以方程（4）为目标函数，以方程（5）为约束条件，LMM最小化以下成本函数：

（f） (λ_{1} {, λ}_{2}, \dots {, λ}_{n个}, μ) = \frac{1}{2} E类 {{[Z轴 ({v（v）}_{0}) - \sum_{我 = 1}^{n个} λ_{我} Z轴 ({v（v）}_{我})]}^{2}} + μ (1 - \sum_{我 = 1}^{n个} λ_{我})

(6)

哪里μ是拉格朗日乘数。在成本函数的最小点（f）对于它的每个变量都是零。因此，优化问题分解为求解以下方程组之一：

{\begin{matrix} \frac{\partial （f）}{\partial λ_{我}} = 0, 我 = 1, 2, \dots, n个, \\ \frac{\partial （f）}{\partial μ} = 0, \end{matrix}

(7)

区分成本函数，我们有：

{\begin{matrix} \frac{\partial （f）}{\partial λ_{我}} = λ_{我} E类 [{Z轴}^{2} ({v（v）}_{我})] + \underset{j个 \neq 我}{\sum_{j个 = 1}^{n个}} λ_{j个} E类 [Z轴 ({v（v）}_{我}) Z轴 ({v（v）}_{j个})] - E类 [Z轴 ({v（v）}_{0}) Z轴 ({v（v）}_{我})] - μ = 0, \\ \frac{\partial （f）}{\partial μ} = 1 - \sum_{我 = 1}^{n个} λ_{我} = 0, \end{matrix}

(8)

如果我们知道

E类 [{Z轴}^{2} ({v（v）}_{我})]

,

E类 [Z轴 ({v（v）}_{0}) Z轴 ({v（v）}_{我})]

、和

E类 [Z轴 ({v（v）}_{我}) Z轴 ({v（v）}_{j个})]

，然后可以求解方程。这些值由中的半变异函数估计第2.3.1节（b） ●●●●。

误差最小方差(

σ_{κ}^{2}

)如式（4）所示，可作为估算中的质量指标[39]. 它可以从算法本身评估不确定性的内在估计。

（b）估计半方差和半方差图

半方差和半方差图是地质统计学中的重要概念，包含空间相关性信息。变量在某些位置的半方差由半方差函数估计，半方差函数是两个位置之间距离的函数。通常，对观测数据应用去趋势预处理。经过预处理后，假设变量的空间分布是平稳的，这意味着半方差不随位置变化。基于这个假设，可以从随机变量在空间上作为分离距离的函数具有良好相关性的数据中估计半方差。半方差(γ)第页，共页Z轴在两个数据点之间定义为：

γ ({x个}_{我}, {x个}_{0}) = γ (小时) = \frac{1}{2} 变量 [Z轴 ({x个}_{我}) - Z轴 ({x个}_{0})]

(9)

哪里小时是点之间的距离

{x个}_{我}

和

{x个}_{0}

、和

γ (小时)

是半变异函数[40].

半变异函数通常根据样本点的统计数据进行估算，如下所示：

\hat{γ} (小时) = \frac{1}{2 n个} \sum_{我 = 1}^{n个} {(Z轴 ({x个}_{我}) - Z轴 ({x个}_{我} + 小时))}^{2}

(10)

哪里n个是由距离分隔的采样点对数小时[41].

由于WSN节点的数量不足以收集稳健的统计数据，因此本文使用航空高光谱遥感反演的土壤湿度图来推导半变异函数。在计算中，我们从土壤湿度图中随机抽取了9000个点，并排除了居民区的这些点。在计算半方差之前，应删除空间趋势并检查平稳性[42]. 考虑到趋势，我们对空间趋势进行了线性拟合，并将其从数据中删除，最后将残差用于计算实验变异函数。空间距离被划分为20个箱，并导出每个箱中点对的平均半方差，以表示箱的半方差值（如所示图4).

本文比较了三种常用的半变异函数模型：球面（Sph）模型、指数（Exp）模型和高斯（Gau）模型，以找到适合半变异函数拟合的最佳模型。模型方程如方程（11）-（13）所示：

γ {(小时)}_{S公司 第页 小时} = {\begin{matrix} {C类}_{0}, 小时 = 0 \\ {C类}_{0} {+ C类}_{1} [1.5 * (\frac{小时}{一}) - 0.5 * {(\frac{小时}{一})}^{三}], 0 \leq 小时 \leq 一 \\ {C类}_{0} {+ C类}_{1}, 小时 > 一 \end{matrix}

(11)

γ {(小时)}_{E类 x个 第页} = {\begin{matrix} {C类}_{0}, 小时 = 0 \\ {C类}_{0} {+ C类}_{1} [{1 - e（电子）}^{- \frac{小时}{一}}], 小时 > 0 \end{matrix}

(12)

γ {(小时)}_{G公司 一 u个} = {\begin{matrix} {C类}_{0}, 小时 = 0 \\ {C类}_{0} {+ C类}_{1} [{1 - e（电子）}^{- {(\frac{小时}{一})}^{2}}], 小时 > 0 \end{matrix}

(13)

哪里

γ (小时)

是半方差；

{C类}_{0}

代表块金，即在零距离处观察到的最小变化或“噪音”；

{C类}_{1}

是结构方差，

{C类}_{0} {+ C类}_{1}

表示门槛方差；和一是地质统计学中表示相关长度的范围。

当空间距离超过一定程度时，由于采样数据量不足，其平均半方差值不稳定。在进行半变异函数拟合之前，应删除这些无效数据。因此，我们仅在以下情况下使用数据小时装配过程中≤4000 m。正如米纳斯尼和麦克布拉特尼建议的那样[43]由于金块可以估计为最短可能间隔距离处的半方差，我们将金块固定为最短间隔30m处的平均半方差，其值为2.0747。这三种型号的装配线如所示图4与RMSE（均方根误差）一起评估模型的拟合性能。从中可以看出图4指数模型比其他两个模型更好地估计了半方差曲线，其RMSE值为0.4791，也是最小的。因此，我们最终选择了指数模型（带参数

{C类}_{0} = 2.0747

,

{C类}_{1} = 9.7697

,

一 = 2400 米

)作为半变差函数模型。在此，我们假设该半变异函数可用于内插2012年6月10日至2012年7月15日期间WSN测量的土壤水分。

（c）带辅助变量的传统克里格方法

现有的几种方法在克里金理论框架下将辅助信息融入插值中，例如协克里金方法和KED方法。所有这些算法在引入辅助变量时都是在空间空间上运行的。

Co-kriging方法通过易观测变量的相关性来估计目标变量[44]而隐含了目标变量和辅助变量之间的相关性。特别是，当辅助变量的数量较低，并且在研究区域的所有网格中都不可用时，可以使用协克里格来改进插值预测。KED方法起源于石油和天然气勘探[45]. 它的漂移不是在通用克里金方法中使用坐标的单项式，而是通过辅助变量在外部定义的[46,47]. 有必要强调，以外部漂移形式纳入的辅助变量应与感兴趣的变量高度线性相关[45].

在本研究中，为了检验所提出的扩展克里金方法插值结果的合理性，我们在接下来的分析中使用了协克里金方法和KED方法作为对比。

2.3.2. 光谱变量的选择

在本研究中，我们选择NDVI和反照率作为辅助信息来辅助插值。它基于以下三个主要考虑因素：（1）NDVI和反照率可以直接或间接反映表层土壤水分的变化[48,49]; （2） NDVI和反照率是几乎所有高分辨率遥感数据源都很容易获得的光谱指数；（3） NDVI和反照率代表了可见光和近红外光谱范围内的大部分遥感信息。虽然选择NDVI和反照率作为光谱变量超出了实际考虑范围，但仍有必要分析土壤水分与NDVI/反照率之间的相关性，如下所示图5.

在图5a、 b、从五个可用的HJ卫星图像中收集NDVI和反照率数据，以与WSN观测到的相应日期的土壤水分数据进行比较。从对比结果可以看出，NDVI和反照率与土壤水分有一定的相关性，但相关系数不显著（绝对值小于0.5）。这一结果的解释可能是NDVI可以反映植被状况，当土壤水分充足时，植被通常生长良好[50]. 然而，在灌溉土地上，所有的农田都得到了足够的灌溉；土壤水分与NDVI的相关性减弱。在植被稀疏的土地上，干土的反照率通常高于湿土[48]. 因此，土壤水分和反照率之间存在相关性，但这种相关性可能会受到土壤类型和植被存在的干扰。实际上，人们通常认为不可能从可见光和近红外遥感数据中估计土壤水分。

2.3.3. 扩展克里格方法以融合遥感信息

为了反映土壤水分空间分布模式的更多细节，我们提出了一种新的算法，将NDVI和反照率等遥感变量纳入基本克里格方法。传统的插值空间是由x个和年协调。新算法将插值空间扩展到组合空间和光谱空间，其中NDVI和反照率被视为坐标，就像x和y一样。组合空间中的距离由空间距离和光谱距离表征，如下所示：

小时 = \sqrt{{Δ x个}^{2} {+ Δ 年}^{2}}

(14)

秒 = \sqrt{{(\frac{Δ NDVI公司}{σ_{NDVI公司}})}^{2} + {(\frac{Δ 反照率}{σ_{反照率}})}^{2}}

(15)

哪里小时是空间距离，

Δ x个

和

Δ 年

是两个采样点之间的坐标差；秒表示光谱距离，

Δ NDVI公司

和

Δ 反照率

是两个采样点之间NDVI和反照率值的差异，以及

σ_{NDVI公司}

和

σ_{反照率}

是两个标准化因子（在本研究中，我们简单地将其值设置为0.1、0.1）。

相应地，将半方差计算和半变异函数模型扩展到空间和光谱组合空间。我们还利用土壤湿度图推导了半方差统计，并在计算之前进行了去趋势处理。空间距离和光谱距离分别被划分为20个箱子。样本数量不足的bin值在半方差图拟合之前被去除。我们仅在以下情况下使用数据小时≤4000 m和秒装配过程中≤2。半方差，作为小时和秒，如所示图6我们估计金块的最短延迟为30米，最短光谱延迟为0.005，这是最大光谱距离的千分之一。金块

({C类}_{0})

值为1.8875。

在拟合半变异函数时，我们还比较了球面（Sph）模型、指数（Exp）模型和高斯（Gau）模型，并计算了它们的RMSE值以选择合适的模型。由于指数模型的RMSE值为0.5612，小于其他模型（RMSE_速度=0.6312，RMSE_高（Gau）=0.9054），我们还使用指数模型进行半方差拟合。扩展的指数模型用以下等式表示：

γ_{1} (小时) {= C类}_{1} [{1 - e（电子）}^{- \frac{小时}{一_{1}}}], 小时 > 0

(16)

γ_{2} (秒) {= C类}_{2} [{1 - e（电子）}^{- \frac{秒}{一_{2}}}], 秒 > 0

(17)

γ (小时, 秒) {= γ}_{1} (小时) {+ γ}_{2} (秒) {+ C类}_{0}

(18)

哪里

γ_{1}

和

γ_{2}

是关于小时和秒,γ是总体半变异函数，以及

一_{1}

和

一_{2}

是空间和光谱变量的滞后距离。

使用上述等式（16）–（18）中的半变异函数模型，可以获得拟合半方差图，如所示图7. The

一_{1}

的空间距离预设为2400 m

一_{2}

光谱距离预设为2.5。然后，拟合参数

{C类}_{1}

、和

{C类}_{2}

分别为7.3909和2.9696。

3.结果和讨论

将提出的扩展克里格方法应用于HiWATER活动中的水网测量，以生成2012年6月10日至7月15日4.5 km×5 km绿洲地区的每日土壤水分图。为了评估扩展克里金方法的性能，我们首先目视检查插值结果，并将其与OK方法、Co-Kriging方法和KED方法获得的结果进行比较。然后，利用离开交叉验证方法中的RMSE进行定量不确定度分析。最后，我们插值了2012年6月10日至2012年7月15日期间每天的土壤水分分布，并检查了降雨量和灌溉数据，以确定插值土壤水分的时间变化是否与水分输入一致。

3.1. 空间趋势分析和交叉验证

插值后的土壤水分空间图可以直接反映插值方法获取参数变化模式的能力。我们使用了从第2.3节、2012年7月10日10 cm深处WATERNET的土壤水分数据（45个节点）和2012年7月份10日的遥感数据（NDVI和反照率），通过线性插值获得，以估算试验区的土壤水分分布。用MATLAB语言实现了扩展克里格插值和OK插值；借助NDVI和反照率进行扩展克里格插值，并仅使用WSN数据进行OK插值。利用R语言编写的GSTAT软件包，借助NDVI和反照率，实现了协克里金插值和KED插值。这四种方法的插值结果如所示图8a–d，分别是扩展克里格法、OK法、联合克里格法和KED法的结果。所有结果图都被屏蔽，以排除居民区。显然，与OK和Co-Kriging方法相比，扩展Kriging法的结果显示了更多关于田间土壤水分分布的细节。OK和Cokriging方法的插值结果只能显示平滑且不断变化的趋势面(图8b、 c），而扩展克里金方法的插值结果(图8a）地表有几条粗糙的纹理，田间土壤湿度有细微变化。虽然KED方法的结果(图8d）显示出一些异质性特征，但它仍然不如我们的扩展克里格方法的结果那么明显。

采用留一交叉验证方法来估计插值精度。在留样交叉验证方法中，每个样本都作为测试样本；假设样本数为N，则每次有N−1个训练样本，最后是N个测试结果。这些N个测试结果的RMSE用于评估插值精度。交叉验证结果如所示表1这表明，OK的RMSE值高于扩展克里格法，低于协克里格法和KED法。

3.2. 时间趋势分析及与降雨量/灌溉数据的相关性

除了同时测量来自多个站点的数据外，WSN还能够获得实时和连续的观测结果。这使我们能够在空间和时间域上估计研究区域的连续土壤水分分布。我们插值了2012年6月10日至2012年7月15日期间每天的土壤水分分布。因此，我们得到了36个插值图。由于一些WSN节点出现故障，在此期间可以使用的有效WSN节点的平均数量约为40个，从36个到45个不等。图中显示了六幅间隔七天的地图图9可以看出，土壤含水量的空间分布随时间发生了显著变化，表明有必要对土壤含水量进行动态监测。在插值地图中可以清楚地确定干旱地区及其边界，这对于农业管理是非常有用的信息。

为了验证插值土壤水分图时间序列的真实性，我们选择了研究区域内的六个灌溉农田，并将扩展克里格法和其他三种方法的估计结果与2012年6月10日至2012年7月15日时间序列上的降水和灌溉数据进行了比较。图9显示了六个灌溉场的土壤水分空间分布。比较结果如所示图10如图所示，降雨量记录显示，6月17日、6月26日、6月月27日、7月3日、7月月6日、7月份8日和7月15日有雨。此外，不同农田的灌溉日期不同。很明显，土壤水分值随着降雨量或灌溉量的增加而增加，然后随着农田的干旱而逐渐降低。

遗漏交叉验证也用于准确性评估。图11是所有可用WSN节点的平均RMSE的时间序列，由leave-on-out方法导出。从中可以看出图11OK的RMSE值高于扩展Kriging，而低于Co-Kriging和KED。这意味着，在土壤水分与这两个遥感变量之间相关性较弱的情况下，我们的扩展克里金方法可以利用这两个间接相关变量来改进土壤水分分布的估计，而协克里金方法和KED方法则表现不佳。这表明，当目标变量和辅助变量之间的相关性相对较弱时，我们的扩展克里格方法更适用。

4.结论

随着无线传感器网络等地面对地观测技术的快速发展，我们现在能够实时、连续、多采样点地监测环境参数。然而，仍然需要插值来将点测量扩展到区域内相应参数的空间分布。由于卫星遥感是获取区域地球观测数据的有效方法，因此需要将遥感和地面观测网络的信息结合起来。

本研究提出的扩展克里格方法将遥感图像光谱信息引入到传统的插值方法中。NDVI和反照率是算法中使用的光谱变量。这些光谱变量与空间变量（即x和y）的处理方式相同。因此，插值基本上是相同的克里格算法，但操作在空间维度和光谱维度的组合空间上。将半变异函数模型推广到组合空间。本文利用遥感土壤水分图拟合半变异函数模型。然而，只要数据集足够大，可以导出关于半方差的稳健统计数据，该土壤湿度图就可以由其他样本来源代替。

将该算法应用于绿洲农业区土壤水分传感器网络（WATERNET）获取的土壤水分数据集，该地区是HiWATER运动的重点试验区。由于WATERNET在48个分散点上提供了连续的近地表土壤湿度测量，插值结果为2012年6月10日至2012年7月15日的每日土壤湿度图，覆盖面积约4.5 km×5.0 km。视觉检查表明，所提出的扩展克里格算法的插值结果比OK和Co克里格算法呈现出更多的空间细节，并且比KED算法呈现出更多的对比度。多个灌区长时间序列的田间平均土壤水分与降水量和灌溉量相关，这些水分输入可以很好地解释土壤水分的时间变化。用leave-on-out方法进行定量不确定性分析表明，当目标变量与辅助变量之间的相关性相对较弱时，扩展克里格算法比协克里格算法和KED算法更适用。目前，建议将NDVI和反照率作为光谱变量，以帮助插值，因为它们可以很容易地从大多数高分辨率卫星图像中得出。然而，我们在讨论中证明，可以将更多相关的光谱变量（如LST）纳入该算法以提高其性能。然而，如何选择信息谱变量仍是该算法的一个公开课题。

尽管给出了协克里金法和KED法的一些结果，但出于以下两个考虑，我们不希望在精度方面将所提出的算法与其他复杂算法进行比较：首先，扩展克里金算法比协克里金、KED和数据同化方法简单得多。因此，它可能适用于其他复杂算法的前提条件不满足的情况。其次，作为一种新的算法，扩展克里格法仍需改进。例如，可以细化光谱变量的归一化因子，而组合空间中的半变异函数模型过于简单。然而，我们倾向于向读者展示最简单的算法形式，并将这些改进留给未来的研究人员。

在这篇短文中，我们无法更深入地讨论其他方面。其中之一是WSN测量土壤水分的质量或准确性。正如我们所知，在表面下方4厘米和10厘米处准确安装Hydro Probe II（HP-II）传感器在技术上是困难的；WSN传感器或多或少受到其内部电压和温度的影响；传感器可以在土壤湿度极低或极高时达到饱和状态[26,32]. 同样，遥感数据可能会受到不准确校准和大气校正的影响。除了数据的缺陷外，还应考虑WSN节点的足迹与30 m分辨率遥感像素之间的尺度失配。本研究中的半变异函数不是从野外测量中得出的，而是来自高光谱遥感图像反演的土壤水分图。这两种半变异函数之间的尺度差异很小，可能会给WSN测量值的插值带来额外的不确定性。然而，目前我们没有可靠的数据来支持这些方面的分析。

另一个相关的有趣主题是WSN节点的期望密度，以及如何优化节点的位置。幸运的是，其他研究人员也提出了类似的问题[51]虽然它们的研究领域和目标参数不同。在本文中，我们简单地使用了所有有效的WATERNET数据。

该算法是经典克里金方法的发展。虽然本文建议插值土壤水分数据，但它可能适用于其他环境参数。随着新观测技术的广泛应用，越来越多的高质量测量数据将在小范围内的多个地点提供给公众，扩展克里格法的潜力可以进一步探索。

鸣谢

本研究得到了国家自然科学基金项目（41331171号和41476161号）、中央高校基本科研基金项目（105599GK号）基于WSN的作物生长远程动态监测报警和WSN应用扩展项目（OFSLRSS201626号）的支持国家遥感科学实验室开放基金项目。地面数据由兰州寒冷干旱地区科学数据中心提供。我们还要感谢中国科学院遥感与数字地球研究所（RADI）的研究人员进行了宝贵的讨论，并提供了数据处理支持。

作者贡献

J.Z.和Q.L.构思、设计并执行了研究的主要部分；J.Z.进行了实验并分析了数据；X.L.、B.D.和L.Z.参与了实验，并为材料和分析工具做出了贡献；J.Z.写了这篇论文；Q.L.和R.Y.修订了论文。

利益冲突

作者声明没有利益冲突。

工具书类

奥斯曼医学博士。；Shazali，K.无线传感器网络应用：环境监测系统的研究。Procedia工程。 2012,41, 1204–1210. [谷歌学者] [交叉参考]
Gong，P.无线传感器网络是一种用于环境监测的新型地面遥感技术。J.遥感。 2007,11, 545–551. [谷歌学者]
李，J。；Gao，H.传感器网络研究综述。J.计算。资源开发。 2008,45, 1–15. [谷歌学者]
窦，B。；温，J。；李，X。；刘，Q。；彭杰。；肖，Q。；张，Z。；Tang，Y。；吴，X。；Lin，X.典型地表参数的无线传感器网络及其在粗分辨率遥感像素中的初步应用。Int.J.经销商净敏感度。 2016,2016. [谷歌学者] [交叉参考]
李，X。；刘，Q。；Yang，R。；张，H。；张杰。；蔡，E.叶面积指数传感器的设计与实现。传感器 2015,15, 6250–6269. [谷歌学者] [交叉参考] [公共医学]
李，X。；刘，Q。；Yang，R。；温，J。；张杰。；蔡，E。；Zhang，H.怀来县地面传感网络和卫星遥感的结合。IEEE传感器J。 2016,16, 3819–3826. [谷歌学者] [交叉参考]
李，J。；Heap，A.D.环境科学中空间插值方法比较研究综述：性能和影响因素。经济。通知。 2011,6, 228–241. [谷歌学者] [交叉参考]
Lu，G.Y。；Wong，D.W.一种自适应逆距离加权空间插值技术。计算。地质科学。 2008,34, 1044–1055. [谷歌学者] [交叉参考]
Krige，D.G.《威特沃特斯兰德一些矿山估价和相关问题的统计方法》。1951年，南非约翰内斯堡威特沃特斯兰德大学硕士论文。[谷歌学者]
Matheron，G.地质统计学原理。经济。地质。 1963,58, 1246–1266. [谷歌学者] [交叉参考]
奥利弗，文学硕士。；Webster，R.Kriging：地理信息系统的插值方法。国际地质杂志。信息系统。 1990,4, 313–332. [谷歌学者] [交叉参考]
Y.丸山。；山崎，F。；瑞穗，K。；Tsuchiya，Y。；Yogai，H.基于日本最近地震后损伤数据集的高速公路路堤脆性曲线。土壤动力学。接地。工程师。 2010,30, 1158–1167. [谷歌学者] [交叉参考]
Wu，T。；Li，Y.使用残差克立格法对美国温度进行空间插值。应用。地理。 2013,44, 112–120. [谷歌学者] [交叉参考]
梁，X。；席林，K。；张永康。；Jones，C.Co-Kriging农业河流中硝态氮负荷的估算。水资源管理。 2016,30, 1771–1784. [谷歌学者] [交叉参考]
Cantet，P.使用克里金和外部漂移绘制小岛的月平均降雨量图。西奥。应用。气候。 2017,127, 31–44. [谷歌学者] [交叉参考]
滴落，W。；Bradbury，K。一个简单的每日土壤-水平衡模型，用于估算温带湿润地区地下水补给的时空分布。氢。J。 2007,15, 433–444. [谷歌学者] [交叉参考]
托斯，E。；Farkas，C。；纳吉，V。；A.Hagyo。；Stekauerova，V.土壤-土壤系统中土壤水分状况的空间变化评估。谷物研究委员会。 2008,36, 307–310. [谷歌学者]
王，M。；Shi，S。；林，F。；郝，Z。；江，P。；Dai，G.土壤水分和氮素对东北地区水曲柳幼苗生长和光合反应的影响。公共科学图书馆 2012,7，e30754。[谷歌学者] [交叉参考] [公共医学]
Yu，X。；查·T。；彭，Z。；吴，B。；王，X。；陈，G。；李，C。；曹，J。；贾，G。；Li，X.中国50年生侧柏人工林土壤呼吸对土壤温度和水分的响应。公共科学图书馆 2011,6，e28397。[谷歌学者] [交叉参考] [公共医学]
哥梅兹广场，A。；阿尔瓦雷兹·罗格尔，J。；Albaladejo，J。；Castillo，V.，半干旱环境中不同尺度土壤水分的空间格局和时间稳定性。水文学。程序。 2000,14, 1261–1277. [谷歌学者] [交叉参考]
范，L。；肖，Q。；温，J。；刘，Q。；Tang，Y。；你，D。；Wang，H。；龚，Z。；Li，X.评估用于估算近地表土壤水分的机载CASI/TASI Ts-VI空间方法。远程传感器。 2015,7, 3114–3137. [谷歌学者] [交叉参考]
El Hajj，M。；北巴格达迪。；Zribi，M。；贝劳，G。；切维隆，B。；Courault博士。；Charron，F.使用X波段SAR数据反演灌溉草地上的土壤水分。遥感环境。 2016,176, 202–218. [谷歌学者] [交叉参考]
G.G.蓬努拉加姆。；Jagdhuber，T。；Hajnsek，I。；Rao，Y.S.使用RISAT-1的混合极化SAR数据估算土壤湿度。IEEE传输。地质科学。远程传感器。 2016,54, 2033–2049. [谷歌学者] [交叉参考]
Kang，J。；金·R。；Li，X.基于回归kriging的无线传感器网络土壤水分测量和异质农田多资源遥感信息的放大。IEEE地质科学。遥感快报。 2015,12, 92–96. [谷歌学者] [交叉参考]
Liao，K。；赖，X。；刘，Y。；Zhu，Q.两个典型土地利用山坡近地表土壤水分估算的不确定性分析。J.土壤沉积物 2016, 1–13. [谷歌学者] [交叉参考]
高，S。；朱，Z。；刘，S。；金·R。；Yang，G。；Tan，L.利用遥感辅助数据，基于贝叶斯最大熵方法估算土壤水分的空间分布。国际期刊申请。地球观测地理信息。 2014,32, 54–66. [谷歌学者] [交叉参考]
Cheng，G。黑河流域水生态经济系统的综合管理; 科学出版社：北京，中国，2009。[谷歌学者]
李，X。；Cheng，G。；刘，S。；肖，Q。；马，M。；金·R。；Che，T。；刘，Q。；Wang，W。；齐，Y.黑河流域联合遥测实验研究（HiWATER）：科学目标和实验设计。牛市。Am.Meteorol公司。Soc公司。 2013,94, 1145–1160. [谷歌学者] [交叉参考]
Mengel，D。；Barber，S.田间条件下玉米根系的发育和分布。阿格隆。J。 1974,66, 341–344. [谷歌学者] [交叉参考]
金·R。；李，X。；Yan，B。；李，X。；罗，W。；马，M。；郭杰。；Kang，J。；朱，Z。；Zhao，S.，用于捕获中国黑河流域中游地区表面异质性的嵌套生态水文无线传感器网络。IEEE地质科学。遥感快报。 2014,11, 2015–2019. [谷歌学者] [交叉参考]
金·R。；李，X。；Yan，B。；罗，W。；李，X。；郭杰。；马，M。；Kang，J。；Zhang，Y.黑河流域生态水文无线传感器网络简介。地球科学。 2012,27, 993–1005. [谷歌学者] [交叉参考]
张杰。；刘，Q。；李，X。；牛，H。；蔡，E。；Chang，C.无线传感器网络的校准和数据验证。SPIE公司 2015. [谷歌学者] [交叉参考]
高，B。；贾，L。；Wang，T.通过结合HJ-1A/B反射率观测和MODIS BRDF产品，推导高分辨率下的地表反照率。远程传感器。 2014,6, 8966–8985. [谷歌学者] [交叉参考]
刘，S。；徐，Z。；宋，L。；赵（Q.Zhao）。；Ge，Y。；徐，T。；马云（Ma，Y.）。；朱，Z。；贾，Z。；Zhang，F.在非均匀陆地表面上从多站点到卫星像素尺度的蒸发蒸腾测量。农业。对于。美托洛尔。 2016,10. [谷歌学者] [交叉参考]
徐，Z。；刘，S。；李，X。；Shi，S。；Wang，J。；朱，Z。；徐，T。；Wang，W。；Ma，M.HiWATER-MUSOEXE期间使用的表面能通量测量系统的相互比较。《地球物理学杂志》。Res.大气。 2013,118. [谷歌学者] [交叉参考]
伯吉斯，T。；Webster，R.土壤特性的最佳插值和无规则制图。土壤科学杂志。 1980,31, 333–341. [谷歌学者] [交叉参考]
克瑞格语的起源。数学。地质。 1990,22, 239–252. [谷歌学者] [交叉参考]
Journel，A.G。；C.J.Huijbregts。采矿地质统计学; 学术出版社：美国纽约州纽约市，1978年。[谷歌学者]
Yamamoto，J.K.普通克里格估计可靠性的另一种度量方法。数学。地质。 2000,32, 489–509. [谷歌学者] [交叉参考]
韦伯斯特，R。；文学硕士Oliver。环境科学家地理统计学; 约翰·威利父子公司：美国纽约州纽约市，2007年。[谷歌学者]
巴鲁·P。；麦克唐奈，R。；劳埃德，C。地理信息系统原理; 牛津大学出版社：英国牛津，2014年。[谷歌学者]
齐默尔曼，B。；Zehe，E。；哈特曼，N。；Elsenber，H.分析空间数据：使用土壤水力数据对假设、新方法和不确定性进行评估。水资源。物件。 2008,44. [谷歌学者] [交叉参考]
米纳斯尼，B。；McBratney，A.B.Matérn函数是土壤变异函数的一般模型。Geoderma公司 2005,128, 192–207. [谷歌学者] [交叉参考]
Yates，S。；Warrick，A.使用协克里金法估算土壤含水量。土壤科学。美国社会期刊。 1987,51, 23–30. [谷歌学者] [交叉参考]
Hudson，G。；Wackernagel，H.使用克里金和外部漂移绘制温度图：苏格兰的理论和实例。国际气候杂志。 1994,14, 77–91. [谷歌学者] [交叉参考]
Chiles，J.-P。；Delfiner，P。地理统计学：建模空间不确定性; John Wiley&Sons:Hoboken，新泽西州，美国，2009年；第497卷。[谷歌学者]
Wackernagel，H。多元地质统计学：应用简介; 施普林格科技与商业媒体：德国柏林，2013年。[谷歌学者]
艾德索，S。；杰克逊，R。；里贾纳托，R。；Kimball，B。；Nakayama，F.裸地反照率对土壤含水量的依赖性。J.应用。美托洛尔。 1975,14, 109–113. [谷歌学者] [交叉参考]
彼得斯，A.J。；Walter-Shea，E.A。；纪磊。；维纳，A。；海耶斯，M。；Svoboda，M.D.使用基于NDVI的标准化植被指数进行干旱监测。照片。工程远程传感器。 2002,68, 71–75. [谷歌学者]
Benedetti，R。；Rossini，P.关于NDVI概况作为农业统计工具的使用：艾米利亚-罗马涅小麦产量估计和预测的案例研究。遥感环境。 1993,45, 311–326. [谷歌学者] [交叉参考]
吴，X。；温，J。；肖，Q。；刘，Q。；彭杰。；窦，B。；李，X。；你，D。；Tang，Y。；Liu，Q.非均匀无雪地面上的粗尺度原位反照率观测和验证策略：中国北方MODIS反照率产品初步验证案例。遥感环境。 2016,184, 25–39. [谷歌学者] [交叉参考]

图1。黑河流域中游研究区：图像来自HJ卫星，由波段3（RED）、波段4（NIR）和波段2（GREEN）组合而成；编号的黄色多边形表示研究区域内选定的灌溉田：A为JC-6，B为SQ-1，C为STZ-1，D为WX-4，E为XM-2，F为ZH-6。

图2。2012年6月10日至2012年7月15日的五份可用NDVI地图：(一)6月15日；(b条)6月19日；(c（c）)6月29日；(d日)7月8日；和(e（电子）)7月13日。

图3。2012年6月10日至2012年7月15日，研究区域的降雨量和日平均温度。

图4。采样土壤水分数据的半方差和土壤水分半方差图的拟合曲线：半方差由土壤水分图中的9000个随机采样点计算得出；半变异函数由球面模型、指数模型和高斯模型拟合，数据为小时≤4000m，范围设定为2400m。

图5。辅助变量（NDVI和反照率）与土壤水分的相关性分析：(一)NDVI数据来自五个可用的HJ卫星图像，并与WSN观测到的相应日期的土壤湿度数据进行了比较；(b条)反照率数据来自五个可用的HJ卫星图像，并与WSN观测到的相应日期的土壤湿度数据进行了比较。

图5。辅助变量（NDVI和反照率）与土壤水分的相关性分析：(一)NDVI数据来源于五个可用的HJ卫星图像，并与WSN观测到的相应日期的土壤湿度数据进行了比较；(b条)反照率数据来源于五个可用的HJ卫星图像，并与WSN观测到的相应日期的土壤湿度数据进行了比较。

图6。采样土壤水分数据相对于空间和光谱距离的半方差：X（X）-轴是空间距离Y（Y）-轴是光谱距离，每个网格中的颜色代表土壤水分的平均半方差值；只有有效的半方差数据用于半方差图拟合(小时≤4000 m和秒≤ 2).

图7。使用指数模型在空间和光谱维对土壤水分半变异函数的拟合结果：X（X）-轴是空间距离Y（Y）-轴是光谱距离，每个网格中的颜色表示土壤水分的平均半方差值。

图8。使用WSN测量数据，于2012年7月10日在10 cm深处插值土壤水分（SM）：(一–d日)分别是扩展克里格法、OK法、协同克里格法和KED法的结果；地图中的圆圈表示WATERNET节点。

图9。采用扩展克里格法绘制研究区域内插土壤水分（SM）图：子图对应于6月10日(一)6月17日(b条)6月24日(c（c）)，7月1日(d日)，7月8日(e（电子）)和7月15日(（f）)2012年；图中显示了六个灌溉田：A是JC-6，B是SQ-1，C是STZ-1，D是WX-4，E是XM-2，F是ZH-6。

图9。采用扩展克里格法绘制研究区域内插土壤水分（SM）图：子图对应于6月10日(一)6月17日(b条)6月24日(c（c）)，7月1日(d日)7月8日(e（电子）)和7月15日(（f）)2012年；图中显示了六个灌溉田：A是JC-6，B是SQ-1，C是STZ-1，D是WX-4，E是XM-2，F是ZH-6。

图10。土壤水分（SM）变化趋势与所选灌溉农田的降雨量和灌溉数据的比较：(一)JC-6；(b条)SQ-1；(c（c）)STZ-1；(d日)WX-4；(e（电子）)XM-2；和(（f）)ZH-6；灰色条代表降雨量；红色箭头代表灌溉。

图10。土壤水分（SM）变化趋势与所选灌溉农田的降雨量和灌溉数据的比较：(一)JC-6；(b条)SQ-1；(c（c）)STZ-1；(d日)WX-4；(e（电子）)XM-2；和(（f）)ZH-6；灰色条代表降雨量；红色箭头表示灌溉。

图11。通过扩展克里格法、OK法、Co-Kriging法和KED法计算的RMSE的比较。

表1。2012年7月10日，对扩展克里格法、OK法、协克里格法和KED法的插值结果进行了比较。

表1。2012年7月10日扩展克里格法、OK法、联合克里格法和KED法插值结果的比较。

方法	$σ_{k个}$ (%)	RMSE（%）
扩展克里格	2.7684	2.4185
克里金	2.9234	2.6218
协同克立格	2.9438	3.0110
KED公司	2.9454	3.1686

分享和引用

MDPI和ACS样式

张杰。；李，X。；Yang，R。；刘，Q。；赵，L。；B·窦。无线传感器网络测量近地表土壤水分数据的扩展克里格插值方法。传感器 2017,17, 1390.https://doi.org/10.3390/s17061390

AMA风格

张杰，李旭，杨锐，刘强，赵力，窦波。无线传感器网络测量近地表土壤水分数据的扩展克里格插值方法。传感器. 2017; 17(6):1390.https://doi.org/10.3390/s17061390

芝加哥/图拉宾风格

张嘉林、李秀红、杨荣进、刘强、赵龙、窦宝成。2017.“一种扩展的克里金方法，用于插值无线传感器网络测量的近地表土壤水分数据”传感器17，编号6:1390。https://doi.org/10.3390/s17061390

请注意，从2016年第一期开始，该杂志使用文章编号而不是页码。请参阅更多详细信息在这里.

文章菜单