水库水储量反演方法和装置

1.本技术涉及水文水资源技术领域，特别涉及一种水库水储量反演方法和装置。


背景技术：

2.水库水储量变化的遥感反演方法主要包括基于水域面积和水库水位计算水储量的方法，在获得完整的水库水域面积或水位后，借助水库的水位-水储量关系或水域面积-水储量关系来计算水库的水储量。
3.然而提取水库水域面积主要依靠光学影像，光学影像提取的水域面积精度较高，但极易受到云雾污染导致无法获得有效观测，时间分辨率大为降低；提取水库的水位则主要依靠雷达或激光测高卫星，雷达测高卫星反演的水位精度较低，在地形复杂地区的小水库可能不适用，激光测高卫星反演的水位精度较高，但是重访期较长、时间分辨率较低。
4.因此，目前计算水库水储量的方法缺乏普适性、实施条件苛刻、计算结果精度低。


技术实现要素：

5.基于此，本技术提出一种水库水储量反演方法、装置、计算机设备和计算机可读存储介质，以提高水库水储量反演的精度。
6.第一方面，本技术提供了一种水库水储量反演方法，所述方法包括：获取目标水库中目标局部水域的合成孔径雷达sar影像序列；根据所述sar影像序列，确定所述目标局部水域的水域面积序列；获取所述目标水库中水位与局部水域面积对应的第一关系；根据所述第一关系将所述水域面积序列转化为目标水位序列；根据水位-水储量关系曲线和所述目标水位序列，得到所述目标水库的水储量序列。
7.在其中一个实施例中，所述根据所述sar影像序列，确定所述目标局部水域的水域面积序列，包括：通过分类算法对所述sar影像序列进行分类，根据分类结果确定所述sar影像序列中的水域像元，所述水域像元为类别为水体的像元；根据所述sar影像序列中的各所述水域像元，确定所述目标局部水域的水域面积序列。
8.在其中一个实施例中，所述获取所述目标水库中水位与局部水域面积对应的第一关系，包括：根据激光测高卫星和/或雷达测高卫星获取所述目标水库的初始水位序列；根据所述初始水位序列对应的时间信息，从所述水域面积序列中获取所述初始水位序列对应的初始水域面积序列；通过多项式回归对所述目标水库的所述初始水位序列和所述初始水域面积序列进行处理，得到水位与局部水域面积的第一关系。
9.在其中一个实施例中，所述水库水储量反演方法包括：获取所述目标局部水域的多个样本影像对，所述样本影像对包括样本光学影像和样本sar影像；根据所述样本光学影像，确定训练区边界，所述训练区边界为所述目标局部水域中水体和陆地的边界；根据所述样本sar影像得到样本特征，所述样本特征包括垂直-垂直的后向散射系数、垂直-水平的后向散射系数、滑动平均后的垂直-垂直的后向散射系数、滑动平均后的垂直-水平的后向散射系数，高程值和坡度值；根据所述训练区边界从所述样本sar影像中选取训练样本，将所述训练样本的所述样本特征输入随机森林分类器，训练得到所述分类算法在其中一个实施例中，所述根据所述样本光学影像，确定训练区边界，包括：确定所述样本光学影像的混合水体指数灰度图像；采用最大类间方差法，将所述混合水体指数灰度图像转化为二值影像，所述二值影像包括表征水体部分和陆地部分的像元；将所述二值影像中的所述水体部分矢量化，得到训练区边界。
10.在其中一个实施例中，所述通过分类算法对所述sar影像序列进行分类，根据分类结果确定所述sar影像序列中的水域像元，包括：根据所述sar影像序列得到sar影像中各像元的特征向量，所述特征向量包括垂直-垂直的后向散射系数、垂直-水平的后向散射系数、滑动平均后的垂直-垂直的后向散射系数、滑动平均后的垂直-水平的后向散射系数，高程值和坡度值；将所述各像元的所述特征向量输入分类算法，得到所述各像元的分类结果；根据所述分类结果确定所述sar影像序列中的水域像元。
11.在其中一个实施例中，所述根据水位-水储量关系曲线和所述目标水位序列，得到所述目标水库的水储量序列之前，所述方法还包括：从激光测高卫星获取所述目标水库中最高水位以上的激光点云高程数据；根据所述激光点云高程数据对数字高程模型进行校正；从校正后的所述数字高程模型中获取计算范围内每个栅格点的高程值，所述计算范围根据所述目标水库的最大水面范围得到；根据目标水位、所述计算范围内所述栅格点的数目、所述计算范围内每个所述栅格点的高程值，确定所述目标水位对应的目标水储量，进而得到所述目标水库的所述水位-水储量关系曲线。
12.第二方面，本技术还提供了一种水库水储量反演装置，所述装置包括：影像获取模块，用于获取目标水库中目标局部水域的合成孔径雷达sar影像序列；面积计算模块，用于根据所述sar影像序列，确定所述目标局部水域的水域面积序列；关系计算模块，用于获取所述目标水库中水位与局部水域面积对应的第一关系；水位计算模块，用于根据所述第一关系将所述水域面积序列转化为目标水位序列；水储量计算模块，用于根据水位-水储量关系曲线和所述目标水位序列，得到所述
目标水库的水储量序列。
13.在其中一个实施例中，所述面积计算模块，还用于通过分类算法对所述sar影像序列进行分类，根据分类结果确定所述sar影像序列中的水域像元，所述水域像元为类别为水体的像元；根据所述sar影像序列中的各所述水域像元，确定所述目标局部水域的水域面积序列。
14.在其中一个实施例中，所述关系计算模块，还用于根据激光测高卫星和/或雷达测高卫星获取所述目标水库的初始水位序列；根据所述初始水位序列对应的时间信息，从所述水域面积序列中获取所述初始水位序列对应的初始水域面积序列；通过多项式回归对所述目标水库的所述初始水位序列和所述初始水域面积序列进行处理，得到水位与局部水域面积的第一关系。
15.在其中一个实施例中，所述水库水储量反演装置还包括算法训练模块，用于获取所述目标局部水域的多个样本影像对，所述样本影像对包括样本光学影像和样本sar影像；根据所述样本光学影像，确定训练区边界，所述训练区边界为所述目标局部水域中水体和陆地的边界；根据所述样本sar影像得到样本特征，所述样本特征包括垂直-垂直的后向散射系数、垂直-水平的后向散射系数、滑动平均后的垂直-垂直的后向散射系数、滑动平均后的垂直-水平的后向散射系数、高程值和坡度值；根据所述训练区边界从所述样本sar影像中选取训练样本，将所述训练样本的所述样本特征输入随机森林分类器，训练得到所述分类算法。
16.在其中一个实施例中，所述算法训练模块，还用于确定所述样本光学影像的混合水体指数灰度图像；采用最大类间方差法，将所述混合水体指数灰度图像转化为二值影像，所述二值影像包括表征水体部分和陆地部分的像元；将所述二值影像中的所述水体部分矢量化，确定训练区边界。
17.在其中一个实施例中，所述面积计算模块，还用于根据所述sar影像序列得到sar影像中各像元的特征向量，所述特征向量包括垂直-垂直的后向散射系数、垂直-水平的后向散射系数、滑动平均后的垂直-垂直的后向散射系数、滑动平均后的垂直-水平的后向散射系数，高程值和坡度值；将所述各像元的所述特征向量输入分类算法，得到所述各像元的分类结果；根据所述分类结果确定所述sar影像序列中的水域像元。
18.在其中一个实施例中，所述水库水储量反演装置还包括曲线获取模块，用于从激光测高卫星获取所述目标水库中最高水位以上的激光点云高程数据；根据所述激光点云高程数据对数字高程模型进行校正；从校正后的所述数字高程模型中获取计算范围内每个栅格点的高程值，所述计算范围根据所述目标水库的最大水面范围得到；根据目标水位、所述计算范围内所述栅格点的数目、所述计算范围内每个所述栅格点的高程值，确定所述目标水位对应的目标水储量，进而得到所述目标水库的所述水位-水储量关系曲线。
19.第三方面，本技术还提供了一种计算机设备，所述计算机设备包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现上述各方法实施例中的步骤。
20.第四方面，本技术还提供了一种计算机可读存储介质，所述计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述各方法实施例中的步骤。
21.上述水库水储量反演方法、装置、计算机设备和计算机可读存储介质，获取目标水库中目标局部水域的合成孔径雷达sar影像序列根据sar影像序列，确定目标局部水域的水域面积序列；获取目标水库中水位与局部水域面积对应的第一关系；根据第一关系将水域面积序列转化为目标水位序列；根据水位-水储量关系曲线和目标水位序列，得到目标水库的水储量序列。传统的基于光学影像提取完整的水库水域面积，基于激光测高卫星反演的水位，从而借助水库的水位-水储量关系或面积-水储量关系来计算水库的水储量的方法，需要较大运算量，增大了计算成本且引入更多不确定性，时间分辨率和计算精度低。相比于传统方法中提取完整水库面积，本技术获取的是局部水域的水域面积序列，数据量更小，也即降低了计算的数据量，计算更为简单，进而降低了计算资源消耗和计算误差，且由于sar影像从sentinel-1卫星获取，sentinel-1卫星时间分辨率较高，故而本技术最终转化得到的目标水位序列的时间分辨率较高，从而提高了水储量的反演精度。
附图说明
22.图1为一实施例中水库水储量反演方法的流程示意图；图2为一实施例中步骤104的流程示意图；图3为一实施例中步骤106的流程示意图；图4为一实施例中水库水储量反演方法的流程示意图；图5为一实施例中步骤404的流程示意图；图6为一实施例中步骤202的流程示意图；图7为一实施例中水库水储量反演方法的流程示意图；图8为一实施例中水库水储量遥感反演算法的流程图；图9为一实施例中小湾水库的示意图；图10为一实施例中多源遥感反演的小湾水库水储量序列与死库容和总库容示意图；图11为一实施例中小湾水库水位多源遥感反演值和实测值对比示意图；图12为一个实施例中水库水储量反演装置的结构框图；图13为一个实施例中计算机设备的内部结构图。
具体实施方式
23.为了使本技术的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本技术进行进一步详细说明。应当理解，此处描述的具体实施例仅仅用以解释本技术，并不用于限定本技术。
24.在光学遥感、雷达遥感、测高技术、地理信息系统、水文水资源等技术领域中，水库在储存地表水资源方面发挥着至关重要的作用。在过去几十年中，全球已经建造了大量的水库和大坝，用于防洪、发电和灌溉。水库会对流域的径流产生重大影响，并影响地表水资源的时空分布。模型和卫星测高数据显示，水库的季节性储水变化占地表水变化的一半以上。一些水文模型考虑了水库对河流径流的影响，通常使用概念模型模拟水库运行过程，但这可能与实际情况不同。准确监测水库水位和水储量信息有助于认识水库在径流调节和水资源管理中的作用。然而，水库水位和水储量的原位监测数据十分有限，或由于信息保密难
以获取，因此，卫星监测获得水库水量变化是一种有效的方法。
25.水库水储量变化的遥感反演方法主要包括基于水域面积和水库水位的方法，获得完整的水库水域面积或水位，并借助水库的水位-水储量关系或面积-水储量关系来计算水库的水储量。提取水库水域面积主要依靠光学或sar（synthetic aperture radar，合成孔径雷达）影像。光学影像提取的水面面积精度较高，但极易受到云雾污染导致无法获得有效观测，时间分辨率大为降低。提取水库的水位则主要依靠雷达或激光测高卫星。雷达测高卫星反演的水位精度较低，在地形复杂地区的小水库可能不适用，激光测高卫星反演的水位精度较高，但是重访期较长、时间分辨率较低。
26.因此，申请人研究发现基于不受云雾影响的sar影像提取水库水域面积结合测高卫星数据的方法，可以更好地获得高精度、高时间分辨率水库水储量序列。本技术提出的水库水储量反演方法就是基于这一理论框架，水库水储量反演方法解决的核心问题是：（1）如何从sar影像中精确地识别水面范围；（2）如何利用测高卫星数据修正dem以获得高精度水位-库容曲线。
27.目前其他水储量反演算法存在较多缺陷，例如：通过sar影像提取水库的水域面积时，由于sar影像噪声较大，且易受地形影响，难以获得精确的水域范围；直接从sar影像中提取整个水库水面面积需要较大运算量，增大了计算成本其会引入更多不确定性。
28.基于此，本技术实施例提供了一种水库水储量反演方法，以解决上述问题，克服传统水库水储量反演算法缺乏普适性、实施条件苛刻、反演结果精度低等多方面的缺陷，实现低成本、大范围、高效率的水库水储量遥感监测。
29.在一个实施例中，如图1所示，提供了一种水库水储量反演方法，本实施例以该方法应用于服务器进行举例说明，可以理解的是，该方法也可以应用于终端，还可以应用于包括终端和服务器的系统，并通过终端和服务器的交互实现。本实施例中，该方法包括以下步骤：步骤102，获取目标水库中目标局部水域的合成孔径雷达sar影像序列。
30.本技术实施例中，目标水库为待反演水储量的水库，目标局部水域为针对目标水库选取的合适的局部水面面积提取范围内的水域。sar影像序列为随时间变化的多个sar影像。在确定目标局部水域后，可以将经过目标局部水域的多个sar影像组成sar影像序列。
31.其中，举例来说，目标水库为小湾水库这一类较为狭长的水库时，可以在公开数据集（如joint research centre global surface water, jrc gsw或global reservoir and dam, grand）中选择小湾水库较宽阔的区段并加缓冲区作为局部水域面积提取范围（roi，region of interest，感兴趣区域）。如果目标水库的矢量边界包含在grand数据集中，可以直接从grand数据集下载使用，如果不包含在grand数据集中，可以下载jrc gsw数据中目标水库的历史最大水面范围，将历史最大水面范围的栅格文件在gis（地理信息系统，geographic information system或 geo－information system）软件中矢量化后作为目标水库的矢量边界。矢量边界为水库水面范围的矢量文件，常用的文件格式是.shp。sar影像可以从sentinel-1卫星获取。
32.步骤104，根据sar影像序列，确定目标局部水域的水域面积序列。
33.本技术实施例中，水域面积序列为目标水库中目标局部水域的水域面积随时间的变化。可以根据sar影像序列中随时间变化的多个sar影像中目标局部水域的水域面积得到
水域面积序列。
34.步骤106，获取目标水库中水位与局部水域面积对应的第一关系。
35.本技术实施例中，第一关系可以表征目标水库的水位与目标水库处于该水位时，目标局部水域的局部水域面积之间的关系。可以先获取目标水库某些天的水位，和同日期的目标局部水域的局部水域面积，来确定目标水库中水位与局部水域面积对应的第一关系。
36.步骤108，根据第一关系将水域面积序列转化为目标水位序列。
37.本技术实施例中，目标水位序列为目标水库的水位随时间的变化。可以将水域面积序列代入第一关系中，以通过第一关系将水域面积序列中的局部水域面积换算为水位，从而将水域面积序列转化为目标水位序列。目标水位序列与水域面积序列的时间分辨率相同。
38.步骤110，根据水位-水储量关系曲线和目标水位序列，得到目标水库的水储量序列。
39.本技术实施例中，水位-水储量关系曲线为目标水库的水位与水位对应的水储量之间的关系曲线。水位-水储量关系曲线可以由dem（digital elevation model，数字高程模型）计算得到。将目标水位序列代入表征水位-水储量关系曲线的关系公式中，即可得到目标水库的水储量序列，实现目标水库水储量的反演。
40.上述水库水储量反演方法，获取目标水库中目标局部水域的合成孔径雷达sar影像序列根据sar影像序列，确定目标局部水域的水域面积序列；获取目标水库中水位与局部水域面积对应的第一关系；根据第一关系将水域面积序列转化为目标水位序列；根据水位-水储量关系曲线和目标水位序列，得到目标水库的水储量序列。相比于传统方法中提取完整水库面积，本技术获取的是局部水域的水域面积序列，数据量更小，也即降低了计算的数据量，计算更为简单，进而降低了计算资源消耗和计算误差，且由于sar影像从sentinel-1卫星获取，sentinel-1卫星时间分辨率较高，故而本技术最终转化得到的目标水位序列的时间分辨率较高，从而提高了水储量的反演精度。
41.在一个实施例中，如图2所述，步骤104中，根据sar影像序列，确定目标局部水域的水域面积序列，可以包括：步骤202，通过分类算法对sar影像序列进行分类，根据分类结果确定sar影像序列中的水域像元，水域像元为类别为水体的像元。
42.本技术实施例中，分类算法为rf（random forest，随机森林）算法，通过分类算法可以对sar影像序列中每个sar影像进行水体分类，sar影像中每个像元都可以分为水体和陆地两种类型，分类结果可以用来表征像元是否为水体。水域像元为类别为水体的像元。
43.步骤204，根据sar影像序列中的各水域像元，确定目标局部水域的水域面积序列。
44.本技术实施例中，水域面积序列为目标水库中目标局部水域的水域面积随时间的变化。在确定sar影像序列中每个sar影像的水域像元后，可以根据水域像元的个数以及sar影像的分辨率计算得到每个sar影像中的水域面积，从而获得目标局部水域的水域面积序列。
45.其中，各sar影像为从sentinel-1卫星获取，因此水域面积序列的时间分辨率由sentinel-1的重访周期决定，sentinel-1由两颗相同的卫星组成，重访周期不固定，约为7
天。目标局部水域的水域面积序列的时间分辨率为7天。上述通过分类算法对sar影像序列进行分类，以及确定目标局部水域的水域面积序列的过程都在google earth engine（gee）云计算平台完成，可以极大减少本地计算量。
46.本公开实施例，可以通过对sar影响序列分类获取目标局部水域的水域面积序列，相比于传统的基于水库整体面积的反演方法，降低了水储量计算结果的误差和计算资源消耗，显著提升了水库水储量反演的时空分辨率和精度。
47.在一个实施例中，如图3所述，步骤106中，获取目标水库中水位与局部水域面积对应的第一关系，可以包括：步骤302，根据激光测高卫星和/或雷达测高卫星获取目标水库的初始水位序列。
48.本技术实施例中，初始水位序列为根据卫星测高数据获得的目标水库的水位随时间的变化。对于激光测高卫星icesat-2，可以直接从其atl 13（内陆水体高程）数据集中提取目标水库的水位，对目标水库同一天获得的水位，筛除3倍标准差之外的离群值并取剩余水位的中值作为当天的水位，进而获得初始水位序列。如果icesat-2中目标水库的数据量较少，可以补充雷达测高卫星的数据。以jason-3卫星为例，对jason-3卫星波形数据进行阈值法重采样，并进行其他校正以反演目标水库的水位，采用同上方式筛选同一天的水位值，补充进从icesat-2获得的水位数据中，组成初始水位序列。
49.步骤304，根据初始水位序列对应的时间信息，从水域面积序列中获取初始水位序列对应的初始水域面积序列。
50.其中，在获取初始水位序列后，可以根据初始水位序列中包含的时间信息，从水域面积序列中筛选出这些时间下的水域面积，组成初始水域面积序列。初始水域面积序列的时间分辨率与初始水位序列相同。由于激光测高卫星重访期较长，其反演的初始水位时间分辨率较低，而水域面积序列的时间分辨率较高，可以从水域面积序列中获取与初始水位序列中时间相同的局部水域面积，组成初始水位序列对应的初始水域面积序列。
51.步骤306，通过多项式回归对目标水库的初始水位序列和初始水域面积序列进行处理，得到水位与局部水域面积的第一关系。
52.其中，水位与局部水域面积的第一关系用于表征相同时间下水位与局部水域面积的对应关系。本技术实施例并不对多项式回归的具体计算方法进行限定，只要能得到水位与局部水域面积的第一关系即可。
53.本公开实施例，通过目标水库的初始水位序列结合水域面积序列，构建水位-局部水域面积关系，进而便于将水域面积序列转化为目标水位序列，提高水位序列的时间分辨率，以在进行水储量反演时提高水储量序列的时间分辨率。
54.在一个实施例中，如图4所述，水库水储量反演方法还可以包括：步骤402，获取目标局部水域的多个样本影像对，样本影像对包括样本光学影像和样本sar影像。
55.其中，样本光学影像可以从sentinel-2卫星获取，样本sar影像可以从sentinel-1卫星获取。在确定目标局部水域后，可以筛选经过目标局部水域的时间接近的12对样本影像，例如每个样本影像对中样本光学影像和样本sar影像的时间间隔在5日以内。12个样本影像对应尽可能覆盖目标局部水域的水域面积最大/最小的时间点，增加了训练和验证分类算法的样本数目，提升分类算法的稳健性。
56.步骤404，根据样本光学影像，确定训练区边界，其中训练区边界为目标局部水域中水体和陆地的边界。
57.其中，由于光学影像提取的水面面积精度较高，但极易受到云雾污染导致无法获得有效观测，在确定训练区边界前，可以先对样本光学影像进行筛选，得到有效观测像元比例较高的样本光学影像。有效观测像元即为没有被云覆盖的像元，例如可以选择目标局部水域内云覆盖率低于20％的样本光学影像。sentinel-1获取的样本sar影像不受云雾影响往往可以完全覆盖目标局部水域。使用筛选后的样本光学影像可以确定训练区边界，训练区边界即为目标局部水域中水域范围的边界，可以将训练区边界输入google earth engine（gee）云计算平台，用于分类算法的训练。由于sentinel-2卫星对应波段的空间分辨率为10-20 m，故使用样本光学影像确定的训练区边界的空间分辨率较高。
58.步骤406，根据样本sar影像得到样本特征，样本特征包括垂直-垂直的后向散射系数、垂直-水平的后向散射系数、滑动平均后的垂直-垂直的后向散射系数、滑动平均后的垂直-水平的后向散射系数，高程值和坡度值。
59.其中，样本特征用于输入随机森林分类器以训练得到分类算法。sentinel-1 获取的样本sar影像分辨率为10 m，样本sar影像包含目标局部水域内水面和陆地在垂直-垂直（vv，vertical-vertical）和垂直-水平（vh，vertical-horizontal）极化通道的后向散射系数。sentinel-1 sar卫星发射极化后的雷达波，并且在接收回波时也进行一次极化，例如vv指vertical-vertical极化，即发射和接收都采用垂直极化，vh指vertical-horizontal，即发射和接收分别采用垂直和水平极化。样本sar影像中每个像元都具有垂直-垂直的后向散射系数、垂直-水平的后向散射系数，由于水面的后向散射系数往往低于陆地，可以基于此实现水体和陆地的像元分类。
60.然而样本sar影像容易受到地形影响，故需要对样本sar影像进行地形校正，以一定程度消除地形产生的收缩、叠掩和阴影。传统的处理方法往往是对sar影像进行低通滤波。本技术实施例中直接引入了5
×
5滑动平均后的vv和vh后向散射系数，以及30 m分辨率数字高程模型nasadem（national aeronautics and space administration，美国国家航空航天局；digital elevation model，数字高程模型），以消除地形对样本sar影像的影响。5
×
5滑动平均即为使用某个像元及其周围24个像元组成的5x5正方形窗口作为滑动平均窗口，将此窗口内像元后向散射系数值的均值赋值给这个中心像元，从而减小高频噪声。nasadem垂直高度的分辨率为1m，可以体现水库地形变化，nasadem经过处理后可以得到水库库区地形坡度的分布，可以从nasadem获取每个像元的高程值和坡度值，与垂直-垂直的后向散射系数、垂直-水平的后向散射系数、滑动平均后的垂直-垂直的后向散射系数、滑动平均后的垂直-水平的后向散射系数，一起组成样本特征。
61.步骤408，根据训练区边界从样本sar影像中选取训练样本，将训练样本的样本特征输入随机森林分类器，训练得到分类算法。
62.其中，训练样本为从样本sar影像根据训练区边界，在水面范围和水面以外随机选取的像元。在确定目标局部水域的训练区边界后，可以对训练区边界进行缓冲区处理，以保证选取训练样本的区域中水面面积和非水面面积的比例为1：3。缓冲区处理即为将训练区边界向外扩展一定距离，防止某些时刻目标局部水域的水面扩展到了我们之前确认的水面范围之外，导致水域面积提取结果有偏差。选取的训练样本中水体像元与陆地像元的个数
比例为1：3。例如在分别水面范围内外选择5000和15000个像元点作为训练样本，将这些像元点的样本特征作为自变量，作为随机森林分类器（rf）的训练集，以样本特征对应的像元的类型作为标注信息，训练随机森林分类器，得到分类算法。随机森林分类器可以由50棵决策树组成。
63.示例性的，从12组样本sar影像中选取12组训练样本后，使用12折检验方法训练随机森林分类器得到分类算法。12折检验方法（k-fold，k=12）为机器学习算法的常用检验方式，每次选择11组训练样本做训练，剩余1组训练样本做验证，重复12次，获得12个独立的精度结果，综合起来可以用于分类算法的精度评价。
64.本技术实施例中，分类算法的训练过程均可以在gee中完成，极大减少本地计算量。利用样本光学影像确定训练区边界，进而根据训练区边界从样本sar影像中选取训练样本，通过训练样本的6个样本特征训练得到分类算法，消除了地形对sar影像的影像，提高分类算法的精准性。
65.在一个实施例中，如图5所述，步骤404中，根据样本光学影像，确定训练区边界，可以包括：步骤502，确定样本光学影像的混合水体指数灰度图像。
66.本技术实施例中，混合水体指数（mwi，mixed water index）将多个波段数据转化为一个波段，因此mwi的分布就是一个灰度图像，而不是常见的彩色卫星影像。在对样本光学影像进行筛选，得到有效观测像元比例较高的样本光学影像后，样本光学影像的mwi满足下列公式（一）、（二）和（三）。
67.公式（一）公式（二）公式（三）其中，re3、re4、blue、green、nir、swir1、swir2分别代表sentinel-2卫星影像红边3号、红边4号、蓝色、绿色、近红外、短波红外1号和短波红外2号的波段反射率，ndmi（normalized differencemud index）和awei
sh
（automated water extraction index）分别代表归一化差异淤泥指数和自动水体提取指数。
68.步骤504，采用最大类间方差法，将混合水体指数灰度图像转化为二值影像，二值影像包括表征水体部分和陆地部分的像元。
69.其中，使用最大类间方差法将混合水体指数灰度影像转化的二值影像为水体/陆地二值影像，水体部分值为1，陆地部分为0。
70.步骤506，将二值影像中的水体部分矢量化，得到训练区边界。
71.其中，可以在qgis（开源地理信息系统软件 ，quantum geographic information system ）等地理信息系统软件中将二值影像中的水体部分矢量化。同时可以结合二值影像对应的样本光学影像目视调整，以获得高精度的水面范围，作为训练区边界。得到的训练区边界可以输入gee云计算平台。
72.本公开实施例，通过将样本光学影像转化为二值影像以获得训练区边界，以便于后续步骤中分类算法的建立，由于样本光学影像从sentinel-2获取，sentinel-2对应波段的空间分辨率为10-20 m，故本技术实施例中获取的训练区边界的空间分辨率较高，进而可
以提高分类算法的精度，反演更精确的局部水域面积。
73.在一个实施例中，如图6所述，步骤202中，通过分类算法对sar影像序列进行分类，根据分类结果确定sar影像序列中的水域像元，可以包括：步骤602，根据sar影像序列得到sar影像中各像元的特征向量，特征向量包括垂直-垂直的后向散射系数、垂直-水平的后向散射系数、滑动平均后的垂直-垂直的后向散射系数、滑动平均后的垂直-水平的后向散射系数，高程值和坡度值。
74.其中，可以对sar影像序列中的各sar影像和nasadem（national aeronautics and space administration，美国国家航空航天局；digital elevation model，数字高程模型）进行预处理，得到各sar影像中每个像元vv，vh后向散射值，滑动平均后的vv，vh后向散射值、nasadem导出的高程值和坡度值。
75.步骤604，将各像元的特征向量输入分类算法，得到各像元的分类结果。
76.其中，vv，vh后向散射值，滑动平均后的vv，vh后向散射值、nasadem导出的高程值和坡度值共同组成的特征向量可以作为一个像元的自变量，输入分类算法。分类算法输出的分类结果可以表征该像元是否为水体，如果该像元为水体，输出的分类结果为0，否则输出的分类结果为1。
77.步骤606，根据分类结果确定sar影像序列中的水域像元。
78.其中，在获得每个像元的分类结果后，可以根据分类结果确定每个sar影像中为水体的像元个数，即水域像元的个数。
79.本技术实施例中，通过垂直-垂直的后向散射系数、垂直-水平的后向散射系数、滑动平均后的垂直-垂直的后向散射系数、滑动平均后的垂直-水平的后向散射系数，高程值和坡度值6个参数组成的特征向量作为分类算法的输入，以此实现对sar影像序列的分类，提高了分类结果的精确性。
80.在一个实施例中，如图7所述，步骤110中，根据水位-水储量关系曲线和目标水位序列，得到目标水库的水储量序列之前，可以包括：步骤702，从激光测高卫星获取目标水库中最高水位以上的激光点云高程数据。
81.本技术实施例中，水储量的反演计算过程不仅需要水位的输入，也需要目标水库的水位-水储量关系曲线。水位-水储量关系曲线需要通过dem计算。常用的dem是srtm （航天飞机雷达地形测绘使命，shuttle radar topography mission）dem，而在计算目标水库的水位-水储量关系前，首先需要利用激光测高卫星数据对dem数据的系统偏差进行校正。
82.示例性的，首先，可以利用icesat-2卫星的atl 03沿轨道光子高程数据，以及atl 08沿轨道地表和冠层高度数据和phoreal (photon research and engineering analysis library)软件，获得目标水库边界1000 m缓冲区（从目标水库边界向外扩展1000m）内地表和水体表面的激光点云高程数据，该激光点云高程数据的精度高于srtm dem数据。由于激光卫星发射激光光子的量很大，这些光子和地表或水面的交点很多，在三维空间中构成了所谓的点云，这里的每个点有对应的高程，即为激光点云高程数据。从上述内地表和水体表面的激光点云高程数据中提取目标水库最高水位以上的激光点云高程数据。
83.步骤704，根据激光点云高程数据对数字高程模型进行校正。
84.其中，每个激光点所在的空间位置有对应的srtm dem高程（即海拔），该高程数据可能和icesat-2测得的激光点云高程数据（icesat-2测得的高程值更精确）不同。可以将提
取后的激光点云高程数据与数字高程模型srtm dem中对应的高程数据进行比较，在目标水库周围，这个高程差值存在一个均值，通过在srtm dem每个数据点上减去该均值，可以消除srtm dem的整体误差。
85.步骤706，从校正后的数字高程模型中获取计算范围内每个栅格点的高程值，计算范围根据目标水库的最大水面范围得到。
86.其中，srtm dem是30 mx 30 m的栅格数据，每个栅格点（即像素）有一个高程值（可以理解成海拔值，地表到水准面的距离），因此srtm dem描述了目标水库周边的地形。计算范围可以为目标水库的最大水面范围叠加缓冲区（如200 m），以保证不漏算栅格点。目标水库的最大水面范围可以从公开数据集jrc gsw或grand中获取。
87.步骤708，根据目标水位、计算范围内栅格点的数目、计算范围内每个栅格点的高程值，确定目标水位对应的目标水储量，进而得到目标水库的水位-水储量关系曲线。
88.其中，目标水位为选取的需要计算的一个水位，目标水储量为此时目标水库处于目标水位时，目标水库的水储量。
89.示例性的，srtm dem获取于2000年2月，因此对于2000年2月后蓄水的水库，可以直接计算水库的水位-水储量关系，水位-水储量关系曲线计算方法如公式（四）所示。
90.公式（四）其中，h为目标水位，s(h)为目标水位对应的目标水储量，hi为第i个栅格点的srtm dem高程值，n为计算范围内srtm dem栅格点的数目。
91.对于2000年2月之前已经蓄水的水库，可以基于srtm dem计算水面以上的水位面积关系，计算方法如公式（五）所示。
92.公式（五）其中h为目标水位，a(h)为目标水位对应的水面面积，hi为第i个栅格点的srtm dem高程值，n为计算范围内srtm dem栅格点的数目，sgn函数是符号函数，定义如公式（六）所示。
93.公式（六）获得水面以上的水位面积关系后，将a(h)进行多项式拟合并延伸至水位以下，对水位面积关系进行积分即可获得目标水库的水位-水储量关系关系，如公式（七）所示。
94.公式（七）其中，h0满足a(h0)=0，s(h)为目标水位对应的目标水储量。
95.本公开实施例，通过激光测高卫星对数字高程模型进行校正，从校正的数字高程模型中获取的计算范围内每个栅格点的高程值，获得目标水库的水位-水储量关系曲线，以获得高精度的水位-水储量关系曲线，进而提高目标水库的水储量的反演精度。
96.应该理解的是，虽然如上所述的各实施例所涉及的流程图中的各个步骤按照箭头
的指示依次显示，但是这些步骤并不是必然按照箭头指示的顺序依次执行。除非本文中有明确的说明，这些步骤的执行并没有严格的顺序限制，这些步骤可以以其它的顺序执行。而且，如上所述的各实施例所涉及的流程图中的至少一部分步骤可以包括多个步骤或者多个阶段，这些步骤或者阶段并不必然是在同一时刻执行完成，而是可以在不同的时刻执行，这些步骤或者阶段的执行顺序也不必然是依次进行，而是可以与其它步骤或者其它步骤中的步骤或者阶段的至少一部分轮流或者交替地执行。
97.为了便于本技术水库水储量反演方法的进一步理解，参见图8，本技术在此提供水库水储量遥感反演算法的流程图。首先，选择目标水库的局部水面（对狭长的水库，最好选择水面较宽的河段）作为研究区，选取时间接近（往往5天以内间隔）的sentinel-2光学影像和sentinel-1 sar影像对，依据有效观测像元的比例保留sentinel-2成像质量好，云量少的影像。使用混合水体指数和最大类间方差法提取光学影像的局部水域面积作为参考，基于sar影像中vv和vh极化的后向散射系数和数字高程模型（dem），训练随机森林分类器得到分类算法（rf算法），使用分类算法和sar影像提取局部水域面积序列。结合局部水域面积序列，以及雷达和激光测高卫星反演的水库水位，构建水库水位-局部水域面积关系，将水域面积序列转化为目标水位序列。使用激光测高卫星数据修正dem，基于修正后的dem计算水库的水位-水储量关系，结合目标水位序列计算目标水库的水储量时间序列。
98.本公开实施例针对卫星遥感反演水库水储量变化时间分辨率和精度低的问题，提出了一种融合光学、sar遥感影像和雷达、激光卫星测高数据的水库水储量反演方法，可有效利用sar影像提取水库的局部水域面积，相比于传统的基于水库整体面积的反演方法，本技术实施例降低了水储量计算结果的误差和计算资源消耗，显著提升了水库水储量反演的时空分辨率和精度。
99.可以通过反演的水位与水库实际的水位对比，进行水位精度验证。以小湾水库为例，小湾水库位于东经100度，北纬25度，是澜沧江干流梯级水库中第二大的水库，总库容为14.65 km3，死库容为4.75 km3。小湾水库调节方式为年调节，约每年6至11月蓄水，12月至次年5月放水。sentinel-1/2和icesat-2卫星均经过小湾水库，为水储量反演提供了数据基础。用于验证的实测数据是水库水位计数据，可以提供2019年9月以来每日精确水库水位距平值，小湾水库地形和roi的选择参见图9。本技术提供的水库水储量反演方法在小湾水库的水储量反演结果如图10、11所示，将实测水位和遥感反演水位变化对比，得到均方根误差为2.72 m（对应的水储量误差约为0.38 km3），r2约为0.973，拟合直线斜率0.995，证明遥感反演结果基本没有系统性偏差。从图10、11中可以看出，除个别日期外，遥感反演结果对水库水储量变化的捕捉十分准确，具有很高的应用价值。
100.因此，本技术提供的种基于光学和雷达遥感影像及测高卫星的水库水储量反演方法，解决了复杂地形条件下水库的水储量监测，可服务水库调度、河流管理等，并为水库调节径流情况下的缺资料流域水文模拟提供技术基础。本技术实施例的实施基于sentinel-2卫星光学影像，sentinel-1卫星合成孔径雷达影像，激光（icesat-2）和雷达（jason-3）测高卫星数据，以无云光学影像中的局部水域面积为参考，训练sar影像水体分类算法，提取周时间尺度分辨率局部水面面积信息，结合卫星测高数据反演的水位数据和数字高程模型，建立水库水位-水储量关系，计算水库周时间分辨率的水储量。相比于传统的基于光学影像或测高卫星的水储量反演算法，本技术实施例反演的水库水储量具有更高的时间分辨率和
更高的反演精度，基于小湾水库实测水位验证结果显示：水位遥感反演的均方根误差为2.72 m，拟合优度r2达到0.987。本技术实施例适用于各类型水库，但由于sar影像具有一定穿透性且存在噪声，在应用于水位或面积变化极小的水库时存在一定不确定性。水库水储量反演的时间分辨率由sentinel-1的重访周期决定，由于sentinel-1由两颗相同的卫星组成，重访周期不固定，约为7天。
101.本技术实施例提供的水库水储量反演方法主要涉及光学遥感、雷达遥感、测高技术、地理信息系统、水文水资源，可实现低成本、大范围、高效率的水库水储量监测。在掌握水库的水位和水储量相关变化规律的基础上，利用光学，sar影像的水体分类结果，结合雷达和激光测高卫星的水位反演结果及水库的水位-水储量关系，反演水库的水储量。通过卫星遥感、大地测量、水文学等多学科的联合创新，提供可靠的水库水储量估计值，服务于水库监测和调度、流域水文模拟、水资源管理等，可以有效降低水库实地监测的成本。本技术实施例利用多源卫星数据可以获得较为理想的水库水储量反演效果，随着对地观测卫星的更新换代，如2022年即将发射的swot卫星等，本技术的实用性和可靠性有望得到进一步提高，因此本技术具有较高的应用潜力。
102.基于同样的发明构思，本技术实施例还提供了一种用于实现上述所涉及的水库水储量反演方法的水库水储量反演装置。该装置所提供的解决问题的实现方案与上述方法中所记载的实现方案相似，故下面所提供的一个或多个水库水储量反演装置实施例中的具体限定可以参见上文中对于水库水储量反演方法的限定，在此不再赘述。
103.在一个实施例中，参见图12，提供了一种水库水储量反演装置1200。水库水储量反演装置1200包括：影像获取模块1202，用于获取目标水库中目标局部水域的合成孔径雷达sar影像序列；面积计算模块1204，用于根据sar影像序列，确定目标局部水域的水域面积序列；关系计算模块1206，用于获取目标水库中水位与局部水域面积对应的第一关系；水位计算模块1208，用于根据第一关系将水域面积序列转化为目标水位序列；水储量计算模块1210，用于根据水位-水储量关系曲线和目标水位序列，得到目标水库的水储量序列。
104.本技术提供的水库水储量反演装置，获取目标水库中目标局部水域的合成孔径雷达sar影像序列根据sar影像序列，确定目标局部水域的水域面积序列；获取目标水库中水位与局部水域面积对应的第一关系；根据第一关系将水域面积序列转化为目标水位序列；根据水位-水储量关系曲线和目标水位序列，得到目标水库的水储量序列。相比于传统方法中提取完整水库面积，本技术获取的是局部水域的水域面积序列，数据量更小，也即降低了计算的数据量，计算更为简单，进而降低了计算资源消耗和计算误差，且由于sar影像从sentinel-1卫星获取，sentinel-1卫星时间分辨率较高，故而本技术最终转化得到的目标水位序列的时间分辨率较高，从而提高了水储量的反演精度。
105.在一个实施例中，面积计算模块1204还用于通过分类算法对sar影像序列进行分类，根据分类结果确定sar影像序列中的水域像元，水域像元为类别为水体的像元；根据sar影像序列中的各水域像元，确定目标局部水域的水域面积序列。
106.在一个实施例中，关系计算模块1206还用于根据激光测高卫星和/或雷达测高卫
星获取目标水库的初始水位序列；根据初始水位序列对应的时间信息，从水域面积序列中获取初始水位序列对应的初始水域面积序列；通过多项式回归对目标水库的初始水位序列和初始水域面积序列进行处理，得到水位与局部水域面积的第一关系。
107.在一个实施例中，水库水储量反演装置1200还包括算法训练模块，用于获取目标局部水域的多个样本影像对，样本影像对包括样本光学影像和样本合成孔径雷达sar影像；根据样本光学影像，确定训练区边界，训练区边界为目标局部水域中水体和陆地的边界；根据样本sar影像得到样本特征，样本特征包括垂直-垂直的后向散射系数、垂直-水平的后向散射系数、滑动平均后的垂直-垂直的后向散射系数、滑动平均后的垂直-水平的后向散射系数、高程值和坡度值；根据训练区边界从样本sar影像中选取训练样本，将训练样本的样本特征输入随机森林分类器，训练得到分类算法。
108.在一个实施例中，算法训练模块还用于确定样本光学影像的混合水体指数灰度图像；采用最大类间方差法，将混合水体指数灰度图像转化为二值影像，二值影像包括表征水体部分和陆地部分的像元；将二值影像中的水体部分矢量化，确定训练区边界。
109.在一个实施例中，面积计算模块1204还用于根据sar影像序列得到sar影像中各像元的特征向量，特征向量包括垂直-垂直的后向散射系数、垂直-水平的后向散射系数、滑动平均后的垂直-垂直的后向散射系数、滑动平均后的垂直-水平的后向散射系数，高程值和坡度值；将各像元的特征向量输入分类算法，得到各像元的分类结果；根据分类结果确定sar影像序列中的水域像元。
110.在一个实施例中，水库水储量反演装置还包括曲线获取模块，用于从激光测高卫星获取目标水库中最高水位以上的激光点云高程数据；根据激光点云高程数据对数字高程模型进行校正；从校正后的数字高程模型中获取计算范围内每个栅格点的高程值，计算范围根据目标水库的最大水面范围得到；根据目标水位、计算范围内所述栅格点的数目、计算范围内每个栅格点的高程值，确定目标水位对应的目标水储量，进而得到目标水库的所述水位-水储量关系曲线。
111.上述水库水储量反演装置中的各个模块可全部或部分通过软件、硬件及其组合来实现。上述各模块可以硬件形式内嵌于或独立于计算机设备中的处理器中，也可以以软件形式存储于计算机设备中的存储器中，以便于处理器调用执行以上各个模块对应的操作。
112.在一个实施例中，提供了一种计算机设备，该计算机设备可以是服务器，其内部结构图可以如图13所示。该计算机设备包括通过系统总线连接的处理器、存储器和网络接口。其中，该计算机设备的处理器用于提供计算和控制能力。该计算机设备的存储器包括非易失性存储介质和内存储器。该非易失性存储介质存储有操作系统、计算机程序和数据库。该内存储器为非易失性存储介质中的操作系统和计算机程序的运行提供环境。该计算机设备的网络接口用于与外部的终端通过网络连接通信。该计算机程序被处理器执行时以实现一种水库水储量反演方法。
113.本领域技术人员可以理解，图13中示出的结构，仅仅是与本技术方案相关的部分结构的框图，并不构成对本技术方案所应用于其上的计算机设备的限定，具体的计算机设备可以包括比图中所示更多或更少的部件，或者组合某些部件，或者具有不同的部件布置。
114.在一个实施例中，提供了一种计算机设备，包括存储器和处理器，存储器中存储有计算机程序，该处理器执行计算机程序时实现上述各方法实施例中的步骤。
115.在一个实施例中，提供了一种计算机可读存储介质，其上存储有计算机程序，计算机程序被处理器执行时实现上述各方法实施例中的步骤。
116.需要说明的是，本技术所涉及的用户信息（包括但不限于用户设备信息、用户个人信息等）和数据（包括但不限于用于分析的数据、存储的数据、展示的数据等），均为经用户授权或者经过各方充分授权的信息和数据。
117.本领域普通技术人员可以理解实现上述实施例方法中的全部或部分流程，是可以通过计算机程序来指令相关的硬件来完成，所述的计算机程序可存储于一非易失性计算机可读取存储介质中，该计算机程序在执行时，可包括如上述各方法的实施例的流程。其中，本技术所提供的各实施例中所使用的对存储器、数据库或其它介质的任何引用，均可包括非易失性和易失性存储器中的至少一种。非易失性存储器可包括只读存储器（read-only memory，rom）、磁带、软盘、闪存、光存储器、高密度嵌入式非易失性存储器、阻变存储器（reram）、磁变存储器（magnetoresistive random access memory，mram）、铁电存储器（ferroelectric random access memory，fram）、相变存储器（phase change memory，pcm）、石墨烯存储器等。易失性存储器可包括随机存取存储器（random access memory，ram）或外部高速缓冲存储器等。作为说明而非局限，ram可以是多种形式，比如静态随机存取存储器（static random access memory，sram）或动态随机存取存储器（dynamic random access memory，dram）等。本技术所提供的各实施例中所涉及的数据库可包括关系型数据库和非关系型数据库中至少一种。非关系型数据库可包括基于区块链的分布式数据库等，不限于此。本技术所提供的各实施例中所涉及的处理器可为通用处理器、中央处理器、图形处理器、数字信号处理器、可编程逻辑器、基于量子计算的数据处理逻辑器等，不限于此。
118.以上所述实施例的各技术特征可以进行任意的组合，为使描述简洁，未对上述实施例中的各个技术特征所有可能的组合都进行描述，然而，只要这些技术特征的组合不存在矛盾，都应当认为是本说明书记载的范围。
119.以上所述实施例仅表达了本技术的几种实施方式，但并不能因此而理解为对本技术专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本技术构思的前提下，还可以做出若干变形和改进，这些都属于本技术的保护范围。因此，本技术专利的保护范围应以所附权利要求为准。