一种基于二维卷积神经网络的土壤水分产品降尺度方法

1.本发明属于遥感技术领域，具体涉及一种基于二维卷积神经网络的土壤水分产品降尺度方法。


背景技术：

2.土壤水分(soil moisture,sm)是水文学、气象学、农业应用中的重要参数，其时空分布对于研究气候变化、农业干旱、水资源管理等都有着重要的作用。传统的通过站点来监测sm的手段难以获得大面积的sm数据，其应用十分有限。而研究者们发现微波遥感数据与sm有着良好的关系，因此提出了利用微波遥感数据来对sm进行反演。
3.到目前为止，利用被动微波遥感数据反演sm的研究已经取得了很大的突破，相应的sm产品也陆续被发布，例如smap(soil moisture active passive)，smos(soil moisture and ocean salinity)全球sm产品。但是，这些产品的空间分辨率大多在数十公里级，难以满足更为精细化的应用需求，因此需要对这些产品进行降尺度，以提升其空间分辨率。
4.关于sm降尺度方法的研究已十分众多，总结起来，可以分为：(1)基于经验模型的降尺度方法；(2)基于物理模型的降尺度方法；(3)基于空间插值的降尺度方法；(4)基于机器学习的降尺度方法。方法(1)-(3)虽然也可以通过辅助数据引入新的信息，但一般来说它们用到的辅助数据较少，并且依赖的是辅助数据影像中单个像元与降尺度之后的目标影像中单个像元的一一对应关系，也即“像元-像元”的思想。方法(4)虽然可以利用机器学习的方法，加入大量辅助数据进行训练，但目前所使用的机器学习模型依然只是基于“像元-像元”的思想，而忽略了同一景影像中像元与像元之间的内在联系。但显然，一景影像中的每个像元并不是相互独立的。例如，由于坡度与坡向的影响，导致地表所接受的太阳能量有所不同，从而影响地表的蒸散；同时，地表的高度差还会影响土壤微粒的分布和土壤水的渗漏，不同的土壤微粒对sm不同的吸附作用，又会进一步导致sm的空间变异。因此，传统的降尺度方法在本质上丢失了一些有用的信息。
5.二维卷积神经网络(convolutional neural networks,cnn)最初用于图像识别与分类，其可以设置一定大小(例如3*3，5*5，7*7等)的感受野(也即卷积核)，来对整景图像进行学习，如此一来不仅仅可以学习到每个像元包含的信息，而且可以考虑到每个感受野中像元与像元之间的相互作用，这恰恰可以弥补传统的sm降尺度方法的不足。但是，二维cnn却很难应用于遥感反演领域，最主要的原因是，二维cnn网络需要使用大量空间连续的数据来进行训练，而在遥感中这样的数据是难以获取的。
6.因此，本发明旨在针对传统sm降尺度方法无法考虑一景影像中像元与像元之间的信息的不足之处，以及充分考虑了将二维cnn应用于遥感反演中的问题的解决办法，提出了一种基于二维cnn技术的新的sm降尺度方法，并以smap sm产品为例进行说明。
7.检索国内外专利和文献显示，目前尚无同样的降尺度方法发明或发表。


技术实现要素：

8.本发明针对上述传统sm降尺度方法无法考虑一副影像中像元与像元之间信息的不足之处，以及充分考虑了将二维cnn应用于遥感反演中的问题的解决办法，提出了一种基于二维cnn的sm降尺度方法。
9.本发明的技术方案为：
10.一种基于二维卷积神经网络的土壤水分产品降尺度方法，其步骤包括：
11.1)获取目标区域的被动微波土壤水分产品影像以及相应的辅助数据影像；并建立辅助数据影像与被动微波土壤水分产品影像之间的对应关系；
12.2)根据所述对应关系，构建训练数据；然后利用所述训练数据对cnn网络进行训练；
13.3)利用训练后的cnn网络对该目标区域的待处理被动微波土壤水分产品影像进行降尺度处理，得到目标降尺度产品。
14.进一步的，建立辅助数据影像与被动微波土壤水分产品影像之间的对应关系的方法为：首先将辅助数据影像重采样到被动微波土壤水分产品影像的空间分辨率，使得辅助数据影像的像元与被动微波土壤水分产品影像的像元一一对应；设置一个大小为m
×
m的窗口，用以裁剪辅助数据影像，其中m为奇数；然后取辅助数据影像每一窗口作为“图”，将该“图”的中心像元所对应的被动微波土壤水分产品影像作为“像元”，以此建立起“图-像元”的关系作为所述对应关系。
15.进一步的，空间变异性越大的地区对应的被动微波土壤水分产品影像所采用的窗口越大；对于辅助数据影像边缘的像元，在该辅助数据影像边缘周围引入新的像元，并将新像元的值都设为该辅助数据影像中所有像元的平均值；然后对辅助数据影像边缘的像元进行裁剪。
16.进一步的，对引入新像元后的辅助数据影像设置一权重层，该权重层中该辅助数据影像的原始像元权重设为h，引入的新像元的权重为h；其中h大于h；然后将该权重层与引入新像元后的辅助数据影像组合在一起，作为特征层。
17.进一步的，构建所述训练数据的方法为：根据每一“图-像元”关系确定一训练数据，其中所述训练数据中的“像元”为标签，“图”对应的特征层作为特征数据。
18.进一步的，构建所述训练数据的方法为：
19.11)将根据每一“图-像元”关系确定的训练数据中的“像元”对应的被动微波土壤水分产品影像下采样，得到低分辨率的被动微波土壤水分产品影像；
20.12)将低分辨率的被动微波土壤水分产品影像中的每个像元进行复制扩充，再次得到与被动微波土壤水分产品影像相同分辨率的影像；
21.13)将“像元”对应的被动微波土壤水分产品影像与步骤12)得到的影像进行相减，得到δsm影像，并以此作为训练数据的标签；
22.14)将“图-像元”关系确定的训练数据中的“图”对应的辅助数据影像分别重采样，得到被动微波土壤水分产品影像的分辨率下的辅助数据影像t和所述低分辨率下的辅助数据影像t；
23.15)将辅助数据影像t进行复制扩充，使其分辨率与辅助数据影像t分辨率相同；
24.16)将步骤15)处理后的辅助数据影像t中的像元设置权重生成对应的特征层，作
为特征数据。
25.进一步的，利用第i张被动微波土壤水分产品影像前后相邻的共n张被动微波土壤水分产品影像，对该第i张被动微波土壤水分产品影像进行空值点消除。
26.进一步的，依据被动微波土壤水分产品影像随时间变化的快慢确定n值；对于第i张被动微波土壤水分产品影像的每一空值点，取该n张被动微波土壤水分产品影像中该空值点位置对应的像元平均值，作为第i张被动微波土壤水分产品影像中该空值点位置的像元值。
27.进一步的，利用所述训练数据对cnn网络进行训练过程中，空间变异性越大的地区对应的被动微波土壤水分产品影像作为输入时，设置cnn卷积核越大。
28.进一步的，所述辅助数据包括srtm提供的dem数据、modis提供的反射率数据以及地表温度产品数据；得到目标降尺度产品的方法为：利用训练后的cnn网络对该目标区域的待处理被动微波土壤水分产品影像进行预测，然后根据smd＝smo+sm
cnn
得到目标降尺度产品smd；其中，smo表示待处理被动微波土壤水分产品影像，sm
cnn
表示cnn网络预测得到的差值影像δsm。
29.本发明基于二维cnn的sm降尺度方法步骤如下：
30.s1.下载研究区域的被动微波sm产品影像以及相应的辅助数据影像。
31.s2.将所有数据进行6天合成，以消除影像中的空值点。
32.s3.利用6天合成处理后得到的数据产品，建立辅助数据影像与sm影像之间的对应关系。
33.s4.利用s3中建立的对应关系，构建对应的训练数据。
34.s5.搭建cnn网络，并将s4得到的训练数据带入网络中进行训练。
35.s6.利用训练好的网络进行sm的降尺度工作，得到目标降尺度产品。
36.s7.对降尺度产品进行精度验证。
37.进一步地，步骤s1中，辅助数据包括srtm(shuttle radar topography mission)提供的数字高程模型(digital elevation model,dem)数据、中分辨率成像光谱仪(moderate-resolution imaging spectroradiometer,modis)提供的反射率数据以及地表温度产品数据。并且，需要对所有的数据进行预处理，以得到标准平面投影下的geotiff格式遥感影像。
38.进一步地，步骤s2中，对不同的数据要依据其随时间变化的快慢采用不同的数据合成方法，例如对sm这种随时间变化较快的数据采用平均值合成，也即，对于某一个特定的像元，取6天中这一像元的sm的平均值，作为最终6天的合成值。而对植被反射率这种随时间变化较慢的数据要采用ndvi(normalized difference vegetation index)最大值合成，也即，对于某一个特定的像元，取6天中ndvi值达到最大的那一天，将其反射率的值作为最终6天的合成值。
39.进一步地，步骤s3中，对于辅助数据影像与sm影像之间对应关系的建立，为了区别于传统方法中所用的“像元-像元”的思想，本发明采用了一种“图-像元”的方法，并且为了提升精度，在辅助数据中引入了权重层。
40.首先，“图-像元”关系的建立方法如图1所示：
41.首先需要将6天合成后的辅助数据重采样到sm产品的空间分辨率，使得辅助影像
的像元可以与sm影像的像元一一对应。再设计一个大小为m
×
m的窗口(图1中的窗口大小为3
×
3)，用以裁剪辅助产品影像，得到多个窗口；取辅助数据影像每一窗口作为“图”，将该“图”的中心像元所对应的被动微波土壤水分产品影像作为“像元”，以此建立起“图-像元”的关系作为对应关系。其中，m必须为奇数，其大小可以根据经验进行确定，例如对于地表均一的平原地区，可以令m为3或5，对于空间变异性较大的山地等，可以令m为7或9；也可以通过设置多组m，最后通过与验证数据的对比来确定最优值。并且，取该辅助“图”的中心像元所对应的sm“像元”，与该辅助“图”相联系，以此建立起“图-像元”的关系。如图1中的左侧图所示，一个窗口由像元1，2，3，5，6，7，9，10，11这9个像元组成，作为“图”，像元6为“图”的中心像元，像元6对应的sm作为“像元”，以此建立起“图-像元”的关系作为对应关系，如此便可以确定一个3
×
3大小的辅助“图”与sm“像元”的对应关系。
42.对于影像边缘的像元，例如像元1，由于其左侧与上侧已经没有其他的像元存在，因此上述的裁剪方法无法进行。此时则需要将辅助产品影像进行扩充，也即在辅助产品影像周围引入新的像元，如图1中的右侧图所示。这些像元的值都设为辅助产品影像中所有像元的平均值。为了使得最外的像元也可以裁剪出对应的m
×
m窗口，扩充的宽度l应满足：
[0043][0044]
其次，权重层的引入方法如图2所示。
[0045]
根据上文所述，“图-像元”的方法会在辅助影像周围引入实际中并不存在的扩充像元，为了削弱这些扩充像元的影响，本发明同时引入了一个权重层。如图2所示，权重层的大小与扩充之后的辅助影像大小相同，并且将原始像元的权重设为0.95，扩充像元的权重设为0.05。接着，将权重层与辅助影像组合在一起，一起组成特征层。
[0046]
为了方便，在下面的叙述中会将sm产品称为标签，将权重层与辅助影像组合而成的特征层，称为特征。
[0047]
进一步地，步骤s4中，训练数据的构建是建立在“图-像元”关系的基础上的。如图2右下角所示，当确定了特征与标签影像之后，就可以利用步骤s3中所建立的“图-像元”关系，对每一个sm“像元”都可以得到与之相对应的一个辅助特征“图”。将“图”作为cnn的输入，“像元”作为cnn的输出，因此每一组“图-像元”就构成了cnn的一组训练数据。
[0048]
进一步地，步骤s5中，cnn的结构如图3所示，主要参数的设置如表1所示。
[0049]
表1为cnn的主要参数
[0050][0051]
*r表示特征影像的分辨率(km)
[0052]
并且，根据不同的地形情况，可以更改cnn卷积核的大小。例如对于地表较为均一的平原地区，可以将卷积核设置成较小的3
×
3，而对于空间变异性较大的山地等，可以将卷积核设置成5
×
5或7
×
7，以获取更多的同一副影像中像元与像元之间的信息。同时，为了提高降尺度结果的精度，本发明对sm的残差(δsm)进行训练与预测，在此基础上，训练数据的生成方法如图4所示：
[0053]
为了得到δsm，本发明采用了两组不同分辨率sm影像，如图3所示，假设较高分辨率影像的大小为4*4，较低分辨率的影像大小为2*2。对于某一特定的被动微波传感器来说，通常只能提供一种粗分辨率的sm影像，因此首先将原始的4*4大小的sm影像下采样到2*2大小，再将2*2大小影像中的每个像元进行复制扩充，再次得到4*4大小的影像，最后将原始的4*4大小的sm影像与降采样以及复制扩充之后的影像进行相减，得到δsm影像，并以此作为cnn网络的标签。同时，将原始高分辨率的辅助数据影像分别重采样到sm对应的两个分辨率下，得到大小分别为4*4与2*2的辅助数据影像。同样地，将2*2大小的辅助数据影像进行复制扩充变为4*4大小，同时可以根据步骤s3的方法得到4*4大小的权重层，将这二者与4*4大小的影像进行组合，作为cnn网络的特征。最后，利用步骤s4中所述的方法，构建训练数据。
[0054]
进一步地，步骤s6中，模型预测得到的也是降尺度之后的sm影像与原始sm影像的差值，为了得到降尺度之后的sm影像，还需要进行一步计算，也即：
[0055]
smd＝smo+sm
cnn
[0056]
式中，smd表示目标降尺度sm影像，smo表示原始分辨率的sm影像，sm
cnn
表示cnn预测得到的δsm影像。
[0057]
进一步地，步骤s7中，对降尺度结果的验证，需要采用统计中的一些指标，本发明使用的是相关系数r，均方根误差rmse以及无偏均方根误差ubrmse。
[0058]
与现有技术相比，本发明的积极效果为：
[0059]
cnn有效地考虑了辅助影像中像元与像元之间的信息，从而显著地提高了降尺度后土壤水分产品的精度。利用地面站点验证结果表明，在将原始的36km smap sm产品降至3km时，本发明的平均r,rmse,ubrmse分别为0.6728,0.0834,0.0497，而传统的利用支持向量回归方法降尺度的r,rmse,ubrmse分别为0.6523,0.0856,0.0525。
附图说明
[0060]
图1是“图-像元”关系的建立方法图；
[0061]
图2是权重层的设计图。
[0062]
图3是cnn的基本结构图。
[0063]
图4是标签影像与对应特征影像的确定方法图。
具体实施方式
[0064]
为使本发明的上述目的、特征和优点能够更加明显易懂，下面以美国的俄克拉荷马州为研究区，以smap卫星提供的2018年1月1日到2018年12月30日的36km sm产品为例，结合附图对本发明作进一步详细说明。此次示例的目标是通过辅助数据对smap的36km产品进行降尺度，得到9km的sm降尺度影像。
[0065]
s1.需要下载的smap sm产品有两种：smap l3 radiometer global daily 36km ease-grid soil moisture v006(l3_sm_p)和smap enhanced l3 radiometer global daily 9km ease-grid soil moisture v003(l3_sm_p_e)，l3_sm_p_e用来对降尺度之后的影像进行验证。考虑到地表反射率和温度对sm含量的响应，以及地形对sm空间变异的影响，本发明利用srtm(shuttle radar topography mission)提供dem(digital elevation model)数据、modis(moderate-resolution imaging spectroradiometer)提供的反射率数据以及地表温度产品数据作为辅助数据。其中，dem数据的分辨率为30m；modis地表反射率产品共有7个波段，分辨率为500m，其产品编号为mod09ga.006；地表日间温度产品的分辨率为1km，其产品编号为mod11a1.006。下载的l3_sm_p产品的空间范围为规则的9*9网格，覆盖的地表为均一的草地，其余产品的空间范围与其相对应。
[0066]
s2.将sm数据、dem数据以及地表日间温度产品按照平均值进行6天合成，将地表反射率产品以ndvi值为指标，按照ndvi最大值进行6天合成。并且，若合成后的地表温度与地表反射率产品中存在空值点，则用空值点四周最近邻的4个像元的平均值代替
[0067]
s3.“图-像元”关系的建立如图2所示，具体的步骤已在发明内容中进行了描述。
[0068]
s4.标签影像与对应特征影像的确定方法如图4所示，具体的步骤也已在发明内容中进行了描述。首先需要确定的是此次示例中的特征与标签。将36km的l3_sm_p产品升尺度至72km，并将72km的影像进行复制扩充，使其与36km的产品大小相同，再将两景影像相减得到δsm，以此作为cnn的标签；接着分别将dem数据、modis中7个波段的地表反射率数据以及modis地表日间温度数据重采样到36km和72km，同样对72km的影像进行复制扩充，并引入对应的权重层，将这3者进行组合，作为cnn的特征。确定了特征与标签影像之后，通过步骤3中叙述“图-像元”关系的建立方法来生成训练数据集，具体的步骤已在发明内容中进行了描述。
[0069]
s5.本次示例所搭建的cnn如发明内容中所述。由于此次所选的研究区较为均一，因此采用了3
×
3大小的卷积核。搭建好cnn之后，将s4中得到的训练数据集输入，进行网络的训练。
[0070]
s6.在已训练好的cnn的基础上，为了得到9km的sm产品，同样需要先确定预测过程中的特征。首先将dem数据、modis中7个波段的地表反射率数据以及modis地表日间温度数据重采样至36km和9km，并将36km的影像复制扩充至与9km影像相同的大小；接着将两个分
辨率的影像与权重层组合在一起，构成预测过程中的特征；再利用s3中的“图-像元”关系的建立方法，生成预测过程中的特征数据集；最后将这些特征数据集输入到训练好的cnn中，得到预测结果。预测得到的结果是9km sm与36km sm的差值，因此，可以通过下式得到9km的sm产品：
[0071]
sm
9km
＝sm
36km
+sm
cnn
[0072]
式中，sm
36km
为36km的smap sm产品(l3_sm_p)，sm
cnn
为预测过程中cnn得到的结果，sm
9km
为9km的降尺度sm产品。
[0073]
s7.本发明利用了smap提供的9km sm产品(l3_sm_p_e)对降尺度后的产品进行了验证。结果表明，9km的降尺度sm产品与l3_sm_p_e产品对比，r为0.9509，rmse为0.0225，ubrmse为0.0224。
[0074]
最后说明的是，以上实施例仅为本发明的较佳实施例而已，用以说明本发明的技术方案，而非限制本发明，本领域的普通技术人员应当理解，可以对本发明的技术方进行改进或者变换，而不脱离本发明原理的前提下，其均应涵盖在本发明所附权利要求的保护范围当中。