本發(fā)明涉及一種基于集合卡爾曼濾波的地下水估算方法，屬于水文地質(zhì)。

背景技術(shù)：

1、水資源短缺的地區(qū)，其水資源的分布與利用面臨著嚴峻的挑戰(zhàn)，地下水資源不僅是支撐當?shù)厣鐣婧桶l(fā)展的關(guān)鍵，也是緩解水資源短缺、促進農(nóng)業(yè)生產(chǎn)和經(jīng)濟發(fā)展的重要支柱。然而，隨著全球氣候變化的影響加劇，加之人口增長、工業(yè)化進程加速、以及不合理的水資源管理，水資源短缺地區(qū)的地下水系統(tǒng)正遭受前所未有的壓力。氣候變化導(dǎo)致的降水模式變化、溫度升高和極端天氣事件頻發(fā)，直接影響了地下水的補給和儲存。同時，過度抽取地下水以滿足不斷增長的人口需求和農(nóng)業(yè)生產(chǎn)，導(dǎo)致了地下水位下降，水源枯竭的風(fēng)險日益增加。這種情況不僅威脅到地區(qū)居民的生活用水安全，還對農(nóng)業(yè)產(chǎn)量、生態(tài)系統(tǒng)健康和區(qū)域經(jīng)濟發(fā)展構(gòu)成了重大風(fēng)險。然而傳統(tǒng)的地下水儲量估算方法仍存在許多的局限性，例如不能準確反映整體區(qū)域水文動態(tài)變化、地下水補給的時間特征，缺乏準確的地下水狀態(tài)信息也極大限制了水資源規(guī)劃和管理。

技術(shù)實現(xiàn)思路

1、本發(fā)明所要解決的技術(shù)問題是提供一種基于集合卡爾曼濾波的地下水估算方法，能夠準確估算地下水儲量的時空變化，為水資源管理決策提供支持。

2、本發(fā)明為了解決上述技術(shù)問題采用以下技術(shù)方案：本發(fā)明設(shè)計了一種基于集合卡爾曼濾波的地下水估算方法，實現(xiàn)目標區(qū)域地下水儲量的估計，包括如下步驟：

3、步驟a.獲得目標區(qū)域關(guān)于預(yù)設(shè)目標空間網(wǎng)格分辨率對應(yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的陸地水儲量網(wǎng)格分布、土壤含水量網(wǎng)格分布、雪水當量網(wǎng)格分布、冠層水儲量網(wǎng)格分布、降水量網(wǎng)格分布、蒸散量網(wǎng)格分布，以及獲得目標區(qū)域中各水體區(qū)分別對應(yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的水位值、邊界信息，然后進入步驟b；

4、步驟b.基于目標區(qū)域中各水體區(qū)分別對應(yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的邊界信息，根據(jù)各水體區(qū)分別在目標空間網(wǎng)格分辨率下相應(yīng)各網(wǎng)格中對應(yīng)各預(yù)設(shè)歷史時間點的面積，結(jié)合各水體區(qū)分別對應(yīng)各預(yù)設(shè)歷史時間點的水位值，計算獲得各水體區(qū)分別在目標空間網(wǎng)格分辨率下相應(yīng)各網(wǎng)格中對應(yīng)各預(yù)設(shè)歷史時間點的地表水儲量，然后聯(lián)系定義目標空間網(wǎng)格分辨率下其他各網(wǎng)格分別對應(yīng)各預(yù)設(shè)歷史時間點的地表水儲量等于0，構(gòu)建目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率對應(yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的地表水儲量網(wǎng)格分布，然后進入步驟c；

5、步驟c.基于集合卡爾曼濾波，計算目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率對應(yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的全部網(wǎng)格點的地下水儲量估計值，實現(xiàn)目標區(qū)域地下水儲量的估計。

6、作為本發(fā)明的一種優(yōu)選技術(shù)方案：還包括步驟d至步驟f如下，執(zhí)行完步驟c之后，進入步驟d；

7、步驟d.根據(jù)目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率對應(yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點t的網(wǎng)格點地下水儲量估計值按如下公式：

8、

9、獲得目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率對應(yīng)歷史時間點t相較上一歷史時間點t-1的網(wǎng)格點地下水儲量波動估計值其中，1,…,i,…,i，i表示目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率的網(wǎng)格數(shù)，表示目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率第i個網(wǎng)格對應(yīng)歷史時間點t的網(wǎng)格點地下儲水量估計值，表示目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率第i個網(wǎng)格對應(yīng)歷史時間點t-1的網(wǎng)格點地下儲水量估計值，然后進入步驟e；

10、步驟e.基于歷史分析時間段l，按如下公式：

11、

12、獲得目標區(qū)域?qū)?yīng)歷史分析時間段l的地下儲水量潛在補給量rl；

13、同時，按如下公式：

14、

15、獲得目標區(qū)域?qū)?yīng)歷史分析時間段l的地下儲水量潛在排放量dl，||表示取絕對值；

16、然后進入步驟f；

17、步驟f.按如下公式：

18、rnet,l＝rl-dl

19、獲得目標區(qū)域?qū)?yīng)歷史分析時間段l的地下儲水量的凈補給量rnet,l。

20、作為本發(fā)明的一種優(yōu)選技術(shù)方案：所述步驟a包括如下步驟a1至步驟a4；

21、步驟a1.獲得目標區(qū)域?qū)?yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的陸地水儲量網(wǎng)格分布、土壤含水量網(wǎng)格分布、雪水當量網(wǎng)格分布、冠層水儲量網(wǎng)格分布、降水量網(wǎng)格分布、蒸散量網(wǎng)格分布，以及獲得目標區(qū)域中各水體區(qū)分別對應(yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的水位值、邊界信息，然后進入步驟a2；

22、步驟a2.判斷目標區(qū)域中各水體區(qū)分別對應(yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的水位值中是否存在缺失，是則針對該各個缺失的水位值進行填補；否則進入步驟a3；

23、步驟a3.分別針對目標區(qū)域?qū)?yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的土壤含水量網(wǎng)格分布、雪水當量網(wǎng)格分布、冠層水儲量網(wǎng)格分布、降水量網(wǎng)格分布、蒸散量網(wǎng)格分布、陸地水儲量網(wǎng)格分布，更新統(tǒng)一到預(yù)設(shè)目標空間網(wǎng)格分辨率，然后進入步驟a4；

24、步驟a4.判斷目標區(qū)域?qū)?yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的陸地水儲量網(wǎng)格分布中是否存在缺失，是則針對該各個缺失的陸地水儲量進行填補；否則進入步驟b。

25、作為本發(fā)明的一種優(yōu)選技術(shù)方案：所述步驟a2中，判斷目標區(qū)域中各水體區(qū)分別對應(yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的水位值中存在缺失時，則執(zhí)行如下步驟a2-1至步驟a2-4，針對該各個缺失的水位值進行填補；

26、步驟a2-1.分別針對各個未缺失的水位值，執(zhí)行如下：

27、根據(jù)未缺失水位值所對應(yīng)水體區(qū)關(guān)于相應(yīng)歷史時間點的邊界信息，獲得該水體區(qū)基于該邊界信息所涉及降水量網(wǎng)格分布下的各個網(wǎng)格，并獲得該未缺失水位值所對應(yīng)歷史時間點下該各網(wǎng)格降水量的平均值，構(gòu)成該未缺失水位值所對應(yīng)水體區(qū)所在目標區(qū)域的局部區(qū)域關(guān)于相應(yīng)歷史時間點的降水量；

28、同時，根據(jù)該未缺失水位值所對應(yīng)水體區(qū)關(guān)于相應(yīng)歷史時間點的邊界信息，獲得該水體區(qū)基于該邊界信息所涉及蒸散量網(wǎng)格分布下的各個網(wǎng)格，并獲得該未缺失水位值所對應(yīng)歷史時間點下該各網(wǎng)格蒸散量的平均值，構(gòu)成該未缺失水位值所對應(yīng)水體區(qū)所在目標區(qū)域的局部區(qū)域關(guān)于相應(yīng)歷史時間點的蒸散量；

29、即獲得該未缺失水位值所對應(yīng)的降水量與蒸散量，進而獲得各個未缺失水位值分別對應(yīng)的降水量與蒸散量，以未缺失水位值與相對應(yīng)的降水量、蒸散量構(gòu)成第一樣本，獲得各個第一樣本，然后進入步驟a2-2；

30、步驟a2-2.基于各個第一樣本，以第一樣本中降水量、蒸散量為輸入，第一樣本中未缺失水位值為輸出，針對lstm神經(jīng)網(wǎng)絡(luò)進行訓(xùn)練，獲得水位值神經(jīng)網(wǎng)絡(luò)模型，然后進入步驟a2-3；

31、步驟a2-3.分別針對各個缺失的水位值，執(zhí)行如下：

32、根據(jù)缺失水位值所對應(yīng)水體區(qū)關(guān)于相應(yīng)歷史時間點的邊界信息，獲得該水體區(qū)基于該邊界信息所涉及降水量網(wǎng)格分布下的各個網(wǎng)格，并獲得該缺失水位值所對應(yīng)歷史時間點下該各網(wǎng)格降水量的平均值，構(gòu)成該缺失水位值所對應(yīng)水體區(qū)所在目標區(qū)域的局部區(qū)域關(guān)于相應(yīng)歷史時間點的降水量；

33、同時，根據(jù)該缺失水位值所對應(yīng)水體區(qū)關(guān)于相應(yīng)歷史時間點的邊界信息，獲得該水體區(qū)基于該邊界信息所涉及蒸散量網(wǎng)格分布下的各個網(wǎng)格，并獲得該缺失水位值所對應(yīng)歷史時間點下該各網(wǎng)格蒸散量的平均值，構(gòu)成該缺失水位值所對應(yīng)水體區(qū)所在目標區(qū)域的局部區(qū)域關(guān)于相應(yīng)歷史時間點的蒸散量；

34、即獲得該缺失水位值所對應(yīng)的降水量與蒸散量，進而獲得各個缺失水位值分別對應(yīng)的降水量與蒸散量，然后進入步驟a2-4；

35、步驟a2-4.分別針對各個缺失水位值，以缺失水位值所對應(yīng)的降水量與蒸散量為輸入，應(yīng)用水位值神經(jīng)網(wǎng)絡(luò)模型進行處理，獲得相對應(yīng)水位值，對該缺失的水位值進行填補，進而對各個缺失的水位值進行填補。

36、作為本發(fā)明的一種優(yōu)選技術(shù)方案：所述步驟a3中，分別針對目標區(qū)域?qū)?yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的土壤含水量網(wǎng)格分布、雪水當量網(wǎng)格分布、冠層水儲量網(wǎng)格分布、降水量網(wǎng)格分布、蒸散量網(wǎng)格分布、陸地水儲量網(wǎng)格分布，執(zhí)行球諧展開，獲得球諧模型系數(shù)，進而對網(wǎng)格分布進行空間分辨率調(diào)整，更新統(tǒng)一到預(yù)設(shè)目標空間網(wǎng)格分辨率。

37、作為本發(fā)明的一種優(yōu)選技術(shù)方案：所述步驟a3中，執(zhí)行球諧展開的過程中，使用稀疏最小二乘方法使原始場和模擬場之間的殘差最小化，削弱球諧展開中的吉布斯效應(yīng)。

38、作為本發(fā)明的一種優(yōu)選技術(shù)方案：所述步驟a4中，判斷目標區(qū)域?qū)?yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點的陸地水儲量網(wǎng)格分布中存在缺失時，則執(zhí)行如下步驟a4-1至步驟a4-4，針對該各個缺失的陸地水儲量網(wǎng)格值進行填補；

39、步驟a4-1.針對各個未缺失的陸地水儲量網(wǎng)格值，以未缺失陸地水儲量網(wǎng)格值，結(jié)合其所在網(wǎng)格對應(yīng)相應(yīng)歷史時間點的土壤含水量網(wǎng)格值、降水量網(wǎng)格值、蒸散量網(wǎng)格值，構(gòu)成第二樣本，獲得各個第二樣本，然后進入步驟a4-2；

40、步驟a4-2.基于各個第二樣本，以第二樣本中土壤含水量網(wǎng)格值、降水量網(wǎng)格值、蒸散量網(wǎng)格值為輸入，第二樣本中未缺失的陸地水儲量網(wǎng)格值為輸出，針對lstm神經(jīng)網(wǎng)絡(luò)進行訓(xùn)練，獲得陸地水儲量網(wǎng)格值神經(jīng)網(wǎng)絡(luò)模型，然后進入步驟a4-3；

41、步驟a4-3.分別針對各個缺失陸地水儲量網(wǎng)格值，以缺失陸地水儲量網(wǎng)格值所在網(wǎng)格對應(yīng)相應(yīng)歷史時間點的土壤含水量網(wǎng)格值、降水量網(wǎng)格值、蒸散量網(wǎng)格值為輸入，應(yīng)用陸地水儲量網(wǎng)格值神經(jīng)網(wǎng)絡(luò)模型進行處理，獲得相對應(yīng)陸地水儲量網(wǎng)格值，對該缺失的陸地水儲量網(wǎng)格值進行填補，進而對該各個缺失的陸地水儲量網(wǎng)格值進行填補。

42、作為本發(fā)明的一種優(yōu)選技術(shù)方案：所述步驟c包括如下步驟c1至步驟c8；

43、步驟c1.按如下公式：

44、

45、

46、定義目標區(qū)域?qū)?yīng)歷史時間點t的狀態(tài)向量xt，以及目標區(qū)域?qū)?yīng)歷史時間點t的觀測向量yt，其中，1,…,i,…,i，i表示目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率的網(wǎng)格數(shù)，表示目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率第i個網(wǎng)格對應(yīng)歷史時間點t的聯(lián)合土壤含水量網(wǎng)格值、雪水當量網(wǎng)格值、冠層水儲量網(wǎng)格值、地表水儲量網(wǎng)格值、陸地水儲量網(wǎng)格值計算的地下儲水量網(wǎng)格值，表示目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率第i個網(wǎng)格對應(yīng)歷史時間點t的陸地水儲量網(wǎng)格值，[·]t表示轉(zhuǎn)置，然后進入步驟c2；

47、步驟c2.構(gòu)建狀態(tài)向量xt與觀測向量yt之間的相關(guān)方程如下：

48、yt＝htxt+vt

49、其中，ht表示對應(yīng)歷史時間點t與狀態(tài)向量xt同維度的系數(shù)矩陣，vt表示對應(yīng)歷史時間點t與狀態(tài)向量xt同維度的觀測噪聲，然后進入步驟c3；

50、步驟c3.基于預(yù)設(shè)e個互不相同的擾動，針對聯(lián)合土壤含水量、雪水當量、冠層水儲量、地表水儲量，陸地水儲量計算的狀態(tài)向量xt，分別添加各個擾動，構(gòu)建對應(yīng)歷史時間點t的集合xt如下：

51、

52、其中，1,…,e,…,e，表示狀態(tài)向量xt添加第e個擾動的初始量，然后進入步驟c4；

53、步驟c4.分別針對狀態(tài)向量xt-1添加各擾動的初始量執(zhí)行時間正向傳播，獲得各分別對應(yīng)的更新值然后進入步驟c5；

54、步驟c5.按如下公式：

55、

56、

57、構(gòu)建集合矩陣xt|t-1和狀態(tài)協(xié)方差更新值矩陣pt|t-1，其中，表示集合矩陣xt|t-1中各元素的平均值，然后進入步驟c6；

58、步驟c6.按如下公式：

59、

60、計算對應(yīng)歷史時間點t的增益矩陣kt，其中，rt表示對應(yīng)歷史時間點t的觀測向量的協(xié)方差矩陣，然后進入步驟c7；

61、步驟c7.按如下公式：

62、

63、獲得各個的校正值其中，表示觀測向量yt添加第e個擾動的初始量，然后進入步驟c8；

64、步驟c8.按如下公式：

65、

66、獲得對應(yīng)歷史時間點t的地下水儲量估計值即目標區(qū)域關(guān)于目標空間網(wǎng)格分辨率對應(yīng)預(yù)設(shè)歷史時間段內(nèi)各預(yù)設(shè)歷史時間點t的全部網(wǎng)格點的地下水儲量估計值，實現(xiàn)目標區(qū)域地下水儲量的估計。

67、本發(fā)明所述一種基于集合卡爾曼濾波的地下水估算方法，采用以上技術(shù)方案與現(xiàn)有技術(shù)相比，具有以下技術(shù)效果：

68、(1)本發(fā)明所設(shè)計一種基于集合卡爾曼濾波的地下水估算方法，基于目標區(qū)域預(yù)設(shè)各屬性數(shù)據(jù)分布，首先獲得各水體區(qū)分別對應(yīng)各預(yù)設(shè)歷史時間點的水位值、邊界信息，然后構(gòu)建目標區(qū)域?qū)?yīng)各預(yù)設(shè)歷史時間點的地表水儲量網(wǎng)格分布，最后基于集合卡爾曼濾波，計算目標區(qū)域?qū)?yīng)各預(yù)設(shè)歷史時間點的全部網(wǎng)格點的地下水儲量估計值，實現(xiàn)目標區(qū)域地下水儲量的估計；設(shè)計方案能夠準確估算地下水儲量變化，有助于更好地理解氣候變率和人類活動影響下水資源短缺地區(qū)地下水儲量的時空特征，能為可持續(xù)地下水管理、水資源規(guī)劃和分配決策提供支持，并有利于確定因氣候變化或過度開采，可能面臨地下水資源風(fēng)險的潛在地區(qū)；

69、(2)本發(fā)明所設(shè)計一種基于集合卡爾曼濾波的地下水估算方法中，采用了球諧函數(shù)統(tǒng)一了水文數(shù)據(jù)的空間分辨率，并采用稀疏最小二乘法有效削弱了吉布斯效應(yīng)，并通過建立lstm神經(jīng)網(wǎng)絡(luò)重建了缺失的陸地水儲量和水位數(shù)據(jù)，保障了地下水儲量估算結(jié)果的完整性，此外，利用集合卡爾曼濾波估算地下水儲量的方法，納入了觀測數(shù)據(jù)和模型中的不確定性，不再依賴傳統(tǒng)估算方法的先驗假設(shè)和分解流程，提高了估算結(jié)果的精度和可靠性。