基于數據同化的地下水埋深插值研究

馬 歡 岳德鵬 YANG Di 于 強 張啟斌 黃 元

(1.北京林業大學精準林業北京市重點實驗室， 北京 100083； 2.佛羅里達大學地理系， 蓋恩斯維爾 32611)

基于數據同化的地下水埋深插值研究

馬 歡1岳德鵬1YANG Di2于 強1張啟斌1黃 元1

(1.北京林業大學精準林業北京市重點實驗室， 北京 100083； 2.佛羅里達大學地理系， 蓋恩斯維爾 32611)

以西北干旱區典型縣域磴口縣為研究區，以2015年8月份40個地下水采樣點的樣品數據為基礎，引入集合卡爾曼濾波(EnKF)數據同化將其優化作為主變量，以蒸散發量反演結果以及歸一化植被指數(NDVI)數據為協變量，進行協同克里金插值，同時與未采用同化的協同克里金插值結果以及經同化采用普通克里金插值結果進行交叉驗證。結果表明：三者在較大空間尺度上對地下水埋深空間分布趨勢的模擬基本一致，南部沙漠地區整體較高，在空間分布上表現為明顯的地理規律性。同化后的數據進行協同克里金插值的結果改善最顯著，平均誤差、均方根誤差、平均標準誤差均優于未同化插值結果，其中平均誤差僅為0.270 5 m。與普通克里金插值方法相比，協同克里金插值考慮蒸散發與NDVI的協同作用，精度明顯提高，平均誤差減小0.409 7 m，均方根誤差減小0.078 4 m，平均標準誤差減小1.016 7 m。

地下水埋深； 干旱區； 數據同化； 協同克里金法； 插值

引言

地下水是能夠影響地表生態環境系統的一個重要影響因子。地下水的變化往往會引起整個生態環境系統天然平衡狀態的失衡和破壞。過量開采地下水，造成地下水水位下降，在世界范圍內總的來說會誘發4方面的環境地質問題：區域性地下水降落漏斗、地面沉降、海水入侵和水質惡化等。在干旱地區淺層地下水水位大幅度下降，疏干了原有的沼澤濕地，原有的綠洲會變成沙漠，原有的景觀會遭到破壞[1-2]。因此，對西北干旱地區而言，地下水作為水環境和生態環境狀況的重要指標，實時掌握地下水分布變化有著重要的意義。

空間插值方法是研究地下水屬性的空間變異特征和分布的有效工具之一，并已廣泛地應用在地下水資源領域[3-5]。由于通過分布在各處的監測站點只能獲取有限的點尺度上地下水埋深信息，空間上是間斷的，不便于在區域面尺度上開展細節研究。利用空間插值方法對已有站點數據進行處理可以獲得空間連續的地下水埋深數據，但插值結果受制于能獲取資料站點的密度，在站點分布稀疏的區域，地下水數據具有不確定性、離散性高，很難得到精確的插值效果。在多數研究中，往往不重視進行空間數據分析、插值方法的篩選和插值參數的優化。不適合的插值方法和插值參數會造成嚴重扭曲，使數據具有的某種空間變化趨勢難以辨識[6]。

插值方法的比較與參數優化已成為空間插值技術在各領域應用研究的熱點，實際工程中許多變量之間存在一定的相關性，其中一些有較豐富的資料，另一些變量的資料比較少，有的變量很容易測量，有的變量則難以測量或者費時費力。能否利用資料多、容易測量的變量資料來提高資料少、難以測量的變量估計精度，協同克里金法就可以解決此類問題，利用多變量進行高精度估值[7-9]。

地下水埋深的空間分布特征由于存在取樣數據之間自相關和互相關的特性，未知點的埋深屬性可通過特定的方法內插或外推，經典的數理統計方法難以確定周圍的數據樣本點以及相應的插值系數。通常情況下，地下水埋深與蒸散發和NDVI數據在統計量及物理意義上存在某種聯系[10-11]，由于蒸散發和NDVI這2種協同變量的獲取相對簡單，所以可以利用大范圍蒸散發數據、NDVI及部分范圍的地下水埋深數據，利用協同克里金方法估計大范圍的地下水埋深。

本文以西北干旱區典型城市磴口縣為研究區，以2015年8月份的地下水埋深監測數據為基礎數據，采用數據同化技術中的集合卡爾曼濾波算法模擬修正地下水埋深數據，以此為主變量，以蒸散發反演結果以及NDVI作為協變量，進行小樣本地下水埋深協同克里金插值，分析地下水的空間分布規律，以期為干旱區水資源均衡管理和分水政策的制定提供依據。

1 材料與方法

1.1 研究區概況

磴口縣位于內蒙古巴彥淖爾市西南部(106°9′～107°10′E、40°9′～40°57′N)，地處內蒙古河套平原源頭，黃河中上游，背靠狼山山脈，西鄰烏蘭布和沙漠。全縣東西長約92 km，南北寬約65 km，總面積約為4 167 km2。境內海拔高度為1 030～2 046 m，整個地形除山區外東南高西北低，總體漸向鄂爾多斯高原傾斜。磴口縣屬溫帶大陸性季風氣候，歷年平均風速3.0 m/s，瞬間最大風速28 m/s，多年平均降水量143.9 mm，多年平均蒸發量2 327 mm，多年平均氣溫7.6℃。黃河流經磴口縣52 km，年徑流量310億m3，共有水域面積2 406萬m2。河套地區地下水埋深0.5～3 m，沙區地下水埋深3～10 m，山前洪積扇地下水埋深3～30 m，地表水資源較豐富，主要得益于黃河來水，但巴彥淖爾市總體水資源緊缺，時空分布差異大。

1.2 數據來源與處理

地下水埋深樣點共40個，分別采集于磴口縣水務局、中國林科院沙林中心、內蒙古巴彥淖爾市黃灌局設立的監測井，如圖1所示。這些井多數從2006年開始觀測，有較為完整的資料。本文選取2015年8月份地下水埋深數據作為空間插值的原始數據。選取磴口縣2015年8月份Landsat TM影像，利用ENVI5.1軟件對影像進行波段合成、圖像增強和幾何校正處理，通過建立解譯標志與實地調查得到磴口縣2015年土地利用現狀圖以及NDVI分布圖。利用SEBAL模型反演出磴口縣2015年8月份蒸散發數據。其他地表溫度、水文、植被及社會經濟數據來自巴彥淖爾市磴口縣水土保持局、國土資源局、統計局等。選用ArcGlS 10.2作為數據統計分析和插值平臺，建立地理信息數據庫，將數字化底圖和監測井地下水埋深數據導入數據庫。

圖1 磴口縣地下水監測井位置分布Fig.1 Location of groundwater monitoring wells in Dengkou County

1.3 co-Kriging模型

協同克里金法(co-Kriging)建立在協同區域化變量理論基礎之上，利用多個區域化變量之間的互相關性，通過建立交叉協方差函數和交叉變異函數模型，用易于觀測和控制的變量對不易觀測的變量進行局部估計[12]。在本研究中以不易觀測的地下水變量篩選蒸散發以及NDVI作為協同變量，建立交叉變異函數，與普通克里金法相比，可更有效改進地下水數據的估計精度和采樣效率。

取其中一個主要變量作為主變量，即目標變量，本文為地下水埋深數據。同時選取另外的一個或者多個變量作為次要變量，本文選取相關性較高的蒸散發數據以及NDVI數據。主變量的自相關性和其他類型變量之間的交叉相關性都有助于對結果做出更好的預測。這主要因為它不只是對各變量做自相關性預測，還包括對所有變量之間交叉相關性的估計[13]。co-Kriging插值模型步驟為：

(1)

(2)

式中ai、bj——協同克里金權重系數

(1)無偏性條件

(3)

式中μu——蒸散發的平均值μv——NDVI的平均值

(2)最優性條件

在滿足無偏條件下，協同克里金估計方差為

(4)

(5)

式中μ1、μ2——拉格朗日乘子

對F求偏導數并令其為零，得協同克里金線性方程組

(6)

整理此協同克里金線性方程組為

(7)

式(7)是一個n+m+2階線性方程組，解該方程組可得協同克里金權重系數ai、bj，再代入協同克里金線性估計式(2)，得到協同克里金線性無偏最優估計量。

根據協同克里金方程組，地下水協同克里金方差為

(8)

若有多個變量，則求解λak的協同克里金方程組為

(9)

(10)

1.4EnKF-co-Kriging模型

集合卡爾曼濾波(Ensemble Kalman filter，EnKF)通過引入觀測數據利用集合卡爾曼濾波同化算法提高模型模擬精度，首先通過擾動生成模型預報的狀態變量集合，通過卡爾曼濾波進行分析，同時考慮模型預測誤差和觀測值的誤差，利用卡爾曼增益(Kalman lain)對模擬結果進行校正，以每一時刻狀態變量集合的平均值為該時刻狀態變量的最優估計[14-15]。通過向前積分，計算不同時刻的概率密度函數所對應的統計特性(如均值和協方差)。

本文基于EnKF方法的數據同化步驟為：

(1)初始化。給定Q個符合高斯分布的狀態變量，在本文中為各測站點的地下水實測數據。

(11)

式中q——模型狀態變量的個數Q——集合的個數l——時間刻度

(12)

本研究將M()引入協同克里金內插半方差函數模型，通過交叉驗證進行篩選，可以更加精確的模擬地下水埋深。

(3)計算l+1時刻的卡爾曼增益矩陣Ll+1。

(13)

其中

(vl+1～N(0,R))

(14)

(15)

(16)

式中Yl+1——l+1時刻的狀態變量的觀測值vl+1——期望為0、方差為R的高斯白噪聲

(5)判斷是否到結束時刻，即目標年份。如果未到結束時刻，那么返回步驟(2)，否則結束。

1.5 精度驗證

采用交叉驗證方法可以對插值方法的空間插值精度進行比較，在上述研究基礎上，以站點實測值為檢驗點，采用交叉驗證的方法對比3種插值方法的插值精度。交叉檢驗的方法很多，包括相關系數、均值、相對誤差、均方根誤差等。本文將平均誤差(Mean error)、均方根誤差(Root mean square error)、平均標準誤差(Mean standard error)作為評估插值效果的標準[16-19]。

2 結果與分析

2.1 變量空間變異

圖2 蒸散發、NDVI與地下水埋深相關性分析Fig.2 Correlation analysis of evapotranspiration, NDVI and groundwater depth

地下水埋深是地下水環境要素的重要指標，能較好地反映地下水的情況。為了合理利用地下水，需對地下水埋深及其影響因素有較好的認識。植被分布情況、環境和土地利用類型等因素都會對地下水埋深產生影響。由于氣候變化和人類活動導致地下水水位頻繁、快速的波動，進而誘發依賴地下水的植被蒸騰量變化，是全球半干旱地區廣泛存在的一種現象。在依賴地下水植被廣泛分布的地區，植被蒸騰同蒸散發一樣，是地下水排泄不可忽略的水均衡要素。同時，蒸散發量也是判斷植被水分脅迫響應及植被對地下水依賴程度的重要指標。許多學者對此進行了研究，發現半干旱區植被覆蓋狀況及蒸散發對地下水變化的響應作用較大[20-22]，蒸散發數據利用SEBAL模型的反演結果與NDVI均以像元為單位，因此在整個研究區域有大量數據點作為協變量，且均勻分布于整個研究區域。

將40個地下水埋深取樣點的數據與蒸散發數據、NDVI數據進行相關性分析(圖2)，地下水埋深與蒸散發兩者擬合的決定系數R2=0.810 6，地下水埋深與NDVI兩者擬合的系數R2=0.722 5，可以看出，蒸散發、NDVI與地下水埋深均具有強相關性，研究區地下水埋深與蒸散發以及NDVI在統計量及物理意義上存在某種聯系。對三者進行單樣本柯爾莫戈洛夫-斯米爾諾夫(Kolmogorov-Smirnov)檢驗，地下水埋深與蒸散發、NDVI的檢驗均滿足正態分布[23]。基于以上分析，可進一步利用協同克里金理論，將蒸散發數據以及NDVI數據作為本次研究的協變量，對地下水埋深進行空間分析。

利用蒸散發數據、NDVI作為本次研究的協變量，利用基于EnKF同化的40個取樣點的地下水埋深優化數據作為主變量，對地下水埋深進行空間插值研究，利用較容易獲取的大范圍數據補充地下水埋深數據的不足。變異函數應用于分析區域化變量的變異特征及結構性狀，但由于樣本數量有限，通常變異函數值是由有限的實測樣本值計算所得。因此，當定量地描述變量在整個區域的特征時，必須通過變異函數進行推斷[24]。由于實際工作中區域化變量的變異通常很復雜，它可能在不同方向上有不同的變異性，或者在同一方向上包含著不同尺度多層次的變異性，因此在計算出實驗變異函數后，常采用多種理論模型進行擬合[25]。

為了增加方法的準確性，本文采用11種半方差函數模型對地下水埋深插值的結果進行實證對比，精度評價見表1，發現穩定模型是研究區地下水埋深插值的最優模型，說明不同的指標應根據參數設置比較才能確定最優的協同克里金插值。

表1 不同半方差函數模型空間內插結果的精度評價

由于地下水埋深研究是對中尺度范圍的空間內插進行分析，因此，平均誤差和均方根誤差是主要的精度評價指標和模型選擇的依據[26]。在ArcGIS軟件中的交叉驗證是利用參加空間插值的數據組進行的驗證，包含5個評價指標：平均誤差反映估計值誤差的大小，其值越接近0越好；均方根誤差、平均標準誤差、標準化平均誤差三者反映插值模型的反演靈敏度、極值效應，其值越小越好；標準化均方根誤差越接近1則越好。綜合考慮平均標準誤差、標準化平均誤差、標準化均方根誤差，交叉驗證分析數據表明，地下水埋深的半方差函數模型在選擇穩定模型進行時內插精度最高，其經過擬合后的穩定半方差數學方程表示為：

地下水埋深-地下水埋深模型

h=0.001 114 3C0+1.114 3S(0.136 44,2)

地下水埋深-蒸散發模型

h=-0.418 61S(0.136 44,2)

蒸散發-蒸散發模型

h=0.001 045 3C0+1.045 3S(0.136 44,2)

式中h——半方差函數值S——穩定模型半方差數學方程C0——塊金值，反映最小抽樣尺度以下變量的變異性及測量誤差

模型結果見圖3。

圖3 穩定模型模擬協方差分析Fig.3 Covariance analysis of stable model

2.2 地下水埋深插值

圖4 2015年8月磴口縣地下水埋深空間分布圖(單位：m)Fig.4 Spatial distributions of groundwater depth of Dengkou County in August 2015

利用2015年8月1日—8月30日的地下水埋深數據作為同化的初始狀態變量，同化間隔為24 h，并只在表層加誤差(包括初始誤差、模型誤差、觀測誤差)，其他各層誤差均設為0。利用EnKF融合協同克里金插值模型，得到優化的地下水埋深數據，以此作為插值主變量，以蒸散發數據與NDVI數據為協變量，對磴口縣2015年地下水埋深數據進行空間插值，利用Jenks分類方法重分成8類得到2015年地下水埋深空間分布結果(圖4a)。同時分別將未利用數據同化的地下水數據采用協同克里金插值模型的模擬結果(圖4b)、經同化采用普通克里金插值模型的模擬結果(圖4c)，與協同克里金插值模型結果(圖4a)進行對比，并從物理分布規律和統計學誤差方面進行分析。

各插值結果可以直觀地反映該地區地下水因子的整體連續變化情況，總體上來看，三者在較大空間尺度上對地下水埋深空間分布趨勢的模擬基本一致，大范圍地下水埋深數據均比較符合一般意義上的生態埋深標準范圍(0～10 m)，河套地區地下水埋深0.5～3 m，沙區地下水埋深3～10 m。地下水埋深最高值區域出現在南部沙漠地區，其值均大于5 m；最低值區域出現在東部黃河流域附近區域，其值范圍為0.5～2 m，流域附近水資源較為豐富，地下水埋深較淺；東南部磴口縣城區地下水埋深較大，表明隨著城市化發展，人類活動用水較多，導致地下水埋深較深。從圖4中可以看出3種方法插值結果的主要差別在磴口縣西北部，EnKF-co-Kriging模型的結果中西北部的地下水埋深最大，最高值達到了10.83 m，西北部山前監測點較少，由于EnKF-co-Kriging模型考慮NDVI因素，山前植被覆蓋率較高，蒸散發作用較強，使得其插值結果中地下水埋深較深。通過小樣本插值后可更加了解地下水埋深分布范圍和分布面積，對比離散型監測數據可直觀地反映區域地下水分布變化規律。

從較小的空間尺度觀察插值結果，發現EnKF-co-Kriging生成的插值結果圖更為平滑，反映出該方法具有較高的耐抗性，能使等值線的平滑度與精確度之間達到較好的平衡，利于總體趨勢與局部趨勢的良好表現。而另外2種方法產生的等值線會圍繞埋深極大值或極小值采樣點發生大曲率的彎曲，甚至閉合形成小等值線圈，插值結果較為破碎，說明其對極端值敏感，耐抗性低，對總體趨勢的反映產生較大干擾。

2.3 交叉驗證

為了比較上述3種方法的插值精度，需要分別計算地下水埋深40個樣點插值估算值和實測值之間的平均誤差、均方根誤差和平均標準誤差，結果見表2。

表2 3種插值方法交叉驗證結果

從表2中可以看出，在平均誤差和均方根誤差方面，EnKF-co-Kriging方法均小于co-Kriging插值方法，在進行數據同化后，2種插值方法的平均誤差和均方根誤差都有所減小，可見數據同化方法可以大大提高插值精度；同時，EnKF-co-Kriging方法的平均誤差和均方根誤差均小于EnKF-Kriging插值方法，說明蒸散發、NDVI的協同作用對地下水埋深影響較大，對地下水埋深數據的插值精度來說，考慮蒸散發和NDVI的協同作用進行插值在精度上要優于普通的克里金插值法。因此，基于EnKF的協同克里金方法整體要優于另外2種插值方法得到的結果。

綜上所述，在采取一系列插值優化措施后，3種插值方法均能在較大空間尺度上，比較客觀地模擬出地下水埋深的空間分布趨勢，在較小空間尺度上，EnKF-co-Kriging對局部趨勢的模擬最為準確和合理。結合交叉驗證的結果綜合評定，EnKF-co-Kriging是最適合研究區地下水埋深的插值方法。

3 結論

(1) 磴口縣地下水埋深數據利用集合卡爾曼濾波進行同化，將同化的修正值作為協同克里金插值的主變量，依據蒸散發數據以及NDVI數據的協同作用進行協同克里金插值，得到地下水的空間分布變化。同化后的數據進行協同克里金插值的模擬結果改善最顯著，平均誤差小于未同化協同克里金插值0.057 9 m，均方根誤差減小0.048 1 m，平均標準誤差減小0.461 6 m。

(2) 地下水埋深的空間分布特征由于存在取樣數據之間自相關和互相關的特性，利用地下水埋深與蒸散發以及NDVI在統計量及物理意義上存在的某種聯系，將蒸散發數據及NDVI數據作為協變量進行協同克里金插值，可以彌補主變量(地下水埋深)樣本點少、變異方向不明顯等缺點，與普通克里金插值方法相比，協同克里金插值的精度明顯提高，平均誤差減小0.409 7 m，均方根誤差減小0.078 4 m，平均標準誤差減小1.016 7 m。

(3) 運用協同克里金插值方法對地下水數據分析的結果表明，三者在較大空間尺度上對地下水埋深空間分布趨勢的模擬基本一致，大范圍地下水埋深數據均比較符合一般意義上的生態埋深標準范圍(0～10 m)，空間分布規律與實地情況基本一致，南部沙漠地區整體較高，在空間分布上表現為明顯的地理規律性。

1 鄭麗萍. 地下水資源的保護與管理[J].農業與技術，2010(3)：58-61. ZHENG Liping. Protection and management of groundwater resources[J].Agriculture and Technology,2010(3):58-61.(in Chinese)

2 樊自立，陳亞寧，李和平，等. 中國西北干旱區生態地下水埋深適宜深度的確定[J].干旱區資源與環境，2008,22(2)：1-5. FAN Zili, CHEN Yaning, LI Heping, et al. Determination of suitable ecological groundwater depth in arid areas in northwest part of China[J]. Journal of Arid Land Resources and Environment,2008,22(2):1-5.(in Chinese)

3 蔣慶. 地下水時空變化及監測網多目標優化研究[D].武漢：華中科技大學，2008. JIANG Qing. Research on spatio-temporal changes of the groundwater and monitoring network multi-objective optimization[D]. Wuhan:Huazhong University of Science and Technology,2008.(in Chinese)

4 王春梅，孟慶巖，占玉林，等. 西北干旱區農田春小麥蒸散量的空間插值方法[J/OL].農業機械學報，2014，45(11)：172-179.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20141127&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2014.11.027. WANG Chunmei, MENG Qingyan, ZHAN Yulin, et al. Spatial interpolation analysis of spring wheat evapotranspiration measurements in Northwest China[J/OL].Transactions of the Chinese Society for Agricultural Machinery,2014,45(11):172-179. (in Chinese)

5 何紅艷，郭志華，肖文發. 降水空間插值技術的研究進展[J].生態學雜志，2005,24(10)：1187-1191. HE Hongyan, GUO Zhihua, XIAO Wenfa. Review on spatial interpolation techniques of rainfall[J]. Chinese Journal of Ecology,2005,24(10):1187-1191.(in Chinese)

6 文雯，周寶同，汪亞峰，等. 基于輔助環境變量的土壤有機碳空間插值——以黃土丘陵區小流域為例[J].生態學報，2013,33(19)：6389-6397. WEN Wen, ZHOU Baotong, WANG Yafeng, et al. Soil organic carbon interpolation based on auxiliary environmental covariates:a case study at small watershed scale in loess hilly region[J]. Acta Ecologica Sinica,2013,33(19):6389-6397.(in Chinese)

7 王長虹，朱合華，錢七虎. 克里金算法與多重分形理論在巖土參數隨機場分析中的應用[J].巖土力學，2014，35(增刊2)：386-392. WANG Changhong, ZHU Hehua, QIAN Qihu. Application of Kriging methods and multi-fractal theory to estimate of geotechnical parameters spatial distribution[J].Rock and Soil Mechanics, 2014,35(Supp.2):386-392.(in Chinese)

8 DAI H Y, REN L Y, WANG M, et al. Water distribution extracted from mining subsidence area using Kriging interpolation algorithm[J]. Transactions of Nonferrous Metals Society of China,2011，21(Supp.3):723-726.

9 PANG Su, LI Tingxuan, WANG Yongdong, et al. Spatial interpolation and sample size optimization for soil copper(Cu) investigation in cropland soil at county scale using CoKriging[J]. Agricultural Sciences in China,2009,8(11):1369-1377.

10 郝興明，陳亞寧，李衛紅，等. 塔里木河中下游荒漠河岸林植被對地下水埋深變化的響應[J].地理學報，2008，63(11)：1123-1130. HAO Xingming, CHEN Yaning, LI Weihong, et al. Response of desert riparian forest vegetation to groundwater depth changes in the middle and lower Tarim River[J]. Acta Geographica Sinica,2008,63(11):1123-1130.(in Chinese)

11 趙文智，劉鵠. 荒漠區植被對地下水埋深響應研究進展[J].生態學報，2006，26(8)：2702-2708. ZHAO Wenzhi, LIU Hu. Recent advances in desert vegetation response to groundwater table changes[J]. Acta Ecologica Sinica,2006,26(8):2702-2708. (in Chinese)

12 于正軍，董冬冬，宋維琪，等. 相帶控制下協克里金方法孔隙度預測[J].地球物理學進展，2012，27(4)：1581-1587. YU Zhengjun, DONG Dongdong, SONG Weiqi, et al. Porosity prediction with co-Kriging method controlled by sedimentary fades[J]. Progress in Geophysics,2012,27(4):1581-1587.(in Chinese)

13 俞超鋒. 基于線性矩的協克里金插值在降雨頻率分析中的應用[C]∥中國水利學會2013學術年會論文集——S5河口治理與保護，2013：1773-1780.

14 劉昭，周艷蓮，居為民，等. 基于集合卡爾曼濾波同化方法的農田土壤水分模擬[J].應用生態學報，2011，22(11)：2943-2953. LIU Zhao, ZHOU Yanlian, JU Weimin, et al.Simulation of cropland soil moisture based on an ensemble Kalman filter[J]. Chinese Journal of Applied Ecology,2011,22(11):2943-2953. (in Chinese)

15 于強，岳德鵬，YANG Di，等. 基于EnKF-MCRP模型的生態用地擴張模擬研究[J/OL]. 農業機械學報，2016，47(9)：285-293.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160939&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2016.09.039. YU Qiang, YUE Depeng, YANG Di, et al. Simulation on ecological land use expansion based on EnKF-MCRP model[J/OL]. Transactions of the Chinese Society for Agricultural Machinery,2016,47(9):285-293. (in Chinese)

16 田雷，董德明，魏強，等. 3種空間插值方法在道路塵中Pb監測數據統計處理中的應用比較[J].吉林大學學報:理學版，2011，49(5)：964-968. TIAN Lei, DONG Deming, WEI Qiang, et al. Comparison of three spatial interpolation methods statistical processing of lead monitoring data in road dust[J].Journal of Jilin University:Science Edition, 2011,49(5):964-968. (in Chinese)

17 鄭小波，羅宇翔，于飛，等.西南復雜山地農業氣候要素空間插值方法比較[J].中國農業氣象，2008，29(4)：458-462. ZHENG Xiaobo, LUO Yuxiang, YU Fei, et al. Comparisons of spatial interpolation methods for agro-climinate factors in complex mountain areas of Southwest China[J]. Chinese Journal of Agrometeorology,2008,29(4):458-462. (in Chinese)

18 李巍，范文義，毛學剛，等. 降雨量空間插值方法比較研究[J].安徽農業科學，2014，42(12)：3667-3669. LI Wei, FAN Wenyi, MAO Xuegang, et al. A comparative study of spatial interpolation methods for precipitation[J]. Journal of Anhui Agricultural Sciences,2014,42(12):3667-3669. (in Chinese)

19 云安萍，鞠正山，胡克林，等. 基于距離反比法的土壤鹽分三維空間插值研究[J/OL].農業機械學報，2015，46(12)：148-156.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20151221&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2015.12.021. YUN Anping, JU Zhengshan, HU Kelin, et al. Three limensional spatial interpolation of soil salinity based on inverse distance weighting method[J/OL]. Transactions of the Chinese Society for Agricultural Machinery,2015,46(12):148-156. (in Chinese)

20 王利書，悅琳琳，唐澤軍，等. 氣候變化和農業發展對石羊河流域地下水位的影響[J/OL].農業機械學報，2014，45(1)：121-128.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140120&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2014.01.020. WANG Lishu, YUE Linlin, TANG Zejun, et al. Influence of climate change and agricultural development on groundwater level in Shiyang River Basin[J/OL].Transactions of the Chinese Society for Agricultural Machinery,2014,45(1):121-128. (in Chinese)

21 黃金廷. 半干旱區蒸散發對地下水變化響應機制研究[D].西安：長安大學，2013. HUANG Jinting. The responses of evapotranpiration to the groundwater changes in the semi-arid area[D]. Xi’an: Chang’an University,2013. (in Chinese)

22 馬玉蕾，王德，劉俊民，等. 黃河三角洲典型植被與地下水埋深和土壤鹽分的關系[J].應用生態學報，2013，24(9)：2423-2430. MA Yulei, WANG De, LIU Junmin, et al. Relationships between typical vegetations, soil salinity and groundwater depth in the Yellow River Delta of China [J]. Chinese Journal of Applied Ecology,2013,24(9):2423-2430. (in Chinese)

23 徐馳，曾文治，黃介生，等. 基于高光譜與協同克里金的土壤耕作層含水率反演[J].農業工程學報，2014，30(13)：94-103. XU Chi, ZENG Wenzhi, HUANG Jiesheng,et al. Modelling soil water content in plow layer using co-Kriging interpolation based on hyperspectral data[J].Transactions of the CSAE,2014,30(13):94-103. (in Chinese)

24 葉回春，黃珊瑜，張世文，等. 土壤有機碳空間變異性對采樣密度的響應研究[J/OL].農業機械學報，2014，45(12)：215-223.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20141232&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2014.12.032. YE Huichun, HUANG Shanyu, ZHANG Shiwen, et al. Spatial variability response of soil organic carbon to sampling density change[J/OL]. Transactions of the Chinese Society for Agricultural Machinery,2014,45(12):215-223. (in Chinese)

25 張朔，魯學軍，趙英俊，等. 空間變異函數等效應橢圓套合方法及其應用[J].地球信息科學學報，2009，11(3)：342-348.ZHANG Shuo, LU Xuejun, ZHAO Yingjun, et al. Equbalent ellihtic nested structure and its application in variogram nesting[J]. Journal of Geo-information Science,2009,11(3):342-348. (in Chinese)

26 王文，劉永偉，寇小華，等. 基于集合卡爾曼濾波和HYDRUS-1D模型的土壤剖面含水量同化試驗[J].水利學報，2012，43(11)：1302-1311. WANG Wen, LIU Yongwei, KOU Xiaohua, et al. EnKF and HYDRUS-1D based data assimilation experiments for improving soil moisture profile prediction[J]. Journal of Hydraulic Engineering,2012,43(11):1302-1311. (in Chinese)

Interpolation of Groundwater Depth Based on Data Assimilation

MA Huan1YUE Depeng1YANG Di2YU Qiang1ZHANG Qibin1HUANG Yuan1

(1.BeijingKeyLaboratoryofPrecisionForestry,BeijingForestryUniversity,Beijing100083,China2.DepartmentofGeography,UniversityofFlorida,GainesvilleFL32611,USA)

Groundwater monitoring is limited by practical conditions, and only limited monitoring results can be obtained when it is observed. As a kind of geostatistical interpolation method, cooperative Kriging (co-Kriging) method can effectively represent the transformation of discrete point-like information to planar continuous information. Dengkou County, a typical county in the arid region of Northwest China, was selected as the study area. The sampled data from 40 groundwater sampling sites in 2015 was selected as the main variable. And this data optimized by EnKF was used as the basic data of co-Kriging interpolation. The evapotranspiration results and NDVI data were selected as the covariates. Co-Kriging interpolation was carried out by using the sampled data from 40 groundwater sampling sites in August, 2015, as the main variable, which were optimized by EnKF, and the evapotranspiration results and NDVI data were used as the covariates. Meanwhile, the results of co-Kriging interpolation without using EnKF model and Kriging interpolation optimized by EnKF model were used to verify the accuracy. The results showed that the spatial distribution trend of groundwater depth was basically the same at large scale, the value in the southern desert region was higher, and the spatial distribution showed obvious geography regularity. The most significant improvement was achieved with EnKF model. Based on this improvement, the mean error, root mean square error and mean standard error were all better than those without assimilation, with the mean error of 0.270 5 m. Compared with the ordinary Kriging interpolation method, co-Kriging model took the synergistic effect of evapotranspiration and NDVI into consideration, and the precision was obviously improved. The mean error was decreased by 0.409 7 m, the root mean square error was decreased by 0.078 4 m and the mean standard error was decreased by 1.016 7 m. This study can provide a scientific basis for spatial visualization simulation and reasonable management of water resources in arid areas.

groundwater depth; arid region; data assimilation; co-Kriging; interpolation

10.6041/j.issn.1000-1298.2017.04.027

2017-01-11

2017-02-11

國家自然科學基金項目(41371189)和“十二五”國家科技支撐計劃項目(2012BAD16B00)

馬歡(1992—)，女，博士生，主要從事3S技術在生態環境中的應用研究，E-mail： xiaohuan27@163.com

岳德鵬(1963—)，男，教授，博士生導師，主要從事景觀生態學和土地評價研究，E-mail： yuedepeng@126.com

K903

A

1000-1298(2017)04-0206-09