基于最小響應單元的地下水埋深時空分異特征研究

牛 騰 岳德鵬 于佳鑫 于 強 王朋沖 胡雅慧

(1.北京林業大學精準林業北京市重點實驗室， 北京 100083； 2.北京林業大學林學院， 北京 100083)

0 引言

地下水埋深狀態是水文地質要素的綜合反映，它不僅反映了地下水補給、徑流、排泄的隨機變化過程，同時也是判斷地下水是否出現環境問題及其嚴重程度的重要指標。地下水的變化會引起整個生態環境系統天然平衡狀態的失衡和破壞。過量開采地下水將造成地下水埋深下降。淺層地下水埋深大幅度下降，原有的綠洲會變成沙漠，原有的景觀會遭到破壞。地下水埋深是地下水資源量最直接的表現形式，是地下水管理最重要的控制性指標[1]。

翁牛特旗是內蒙古自治區嚴重缺水地區之一，缺水面積占全旗的30%以上。由于特定的自然地理和水文地質條件，地下水一直是該地區的重要水源，為地區經濟的發展發揮了重要的作用[2-3]。但是，由于嚴重缺水，加之長期超量開采，地下水資源匱乏已成為制約該地區經濟發展的生態“瓶頸”。改進水文響應單元是根據流域內土壤、降水、植被等因素劃分的具有相同水文特性的最小水文單元，基于最小水文響應單元的研究方法可在水資源的合理調配、環境風險評估等領域發揮重要作用。將最小水文響應單元理念運用于地下水埋深的模擬，能將柵格單元轉化為矢量單元，將局部的數據模擬至全局，并能較為有效地反映實際情況[4-7]。空間自相關法以空間關聯測度為核心，研究地下水埋深的空間分布特征，發現其空間關聯模式，能較好地從區域角度分析土壤中地下水埋深的空間自相關特征。

采用協同克里金插值方法，使插值結果不僅基于地下水埋深數據的空間地理特性，也對影響地下水埋深的多種因素進行擬合，并對結果進行綜合模擬，但協同克里金插值結果具有一定的范圍局限性[8]。利用水文響應單元理論將空間內的柵格單元轉化為矢量單元，使之更加貼近實際、符合地下水埋深的擴散和匯流狀態。本文采用協同克里金插值方法，以2005—2017年每隔4 a、共4 a地下水埋深數據為主變量，對應年份的NDVI(歸一化植被指數)、降水量和河網密度為協變量，將地下水埋深在研究區空間范圍內呈現。利用空間自相關的方法將賦予地下水埋深插值結果的最小響應單元進行空間分析[9-10]，研究翁牛特旗時間、空間雙尺度上的分異規律及特性，以期對地下水資源的合理開發利用、當地土地規劃利用和農業發展規劃提供參考。

1 材料與方法

1.1 研究區概況

翁牛特旗地處內蒙古自治區赤峰市中部，南與敖漢旗和赤峰市松山區相連，西接克什克騰旗，北與林西縣、巴林右旗、阿魯科爾沁旗為鄰，東與通遼市奈曼旗接壤。位于東經117°49′～120°43′、北緯42°26′～43°25′之間，西部是高寒臺地，中部是起伏平緩、山川交錯的丘陵區；東部為科爾沁沙地。全旗既有沙海和草甸相間的低洼沼澤，也有老哈河與西拉沐淪河交匯地帶流水沖擊形成的平原，構成全旗“五沙四山一分田”的復雜地形。全旗總面積11 889 km2。研究區屬于溫帶半干旱大陸性季風氣候。翁牛特旗歷年平均降水量在300～400 mm之間。西部山區降水量為400 mm，中部丘陵區降水量360 mm，東部荒漠區降水量300 mm。區內地下水時空分布極不均勻，主要通過大氣降水、農田灌溉、干渠入滲和水庫滲漏等進行補給[11]。

1.2 數據來源與處理

地下水埋深樣點共40個，分布于翁牛特旗監測井，如圖1所示。監測井均從2005年開始觀測，觀測至2017年12月，有較為完整的資料。本文選取2005、2009、2013、2017年夏季地下水埋深數據作為空間插值的原始數據。選取翁牛特旗夏季且少云的2005年7月和2009年7月Landsat TM 影像和2013年7月和2017年7月的Landsat OLI影像，利用ENVI 5.1軟件對影像進行波段合成、圖像增強和幾何校正處理，通過建立解譯標志與實地調查得到翁牛特旗2005、2009、2013、2017年土地利用現狀圖以及NDVI分布圖。從空間地理數據云下載翁牛特旗DEM，用于水文分析中水文響應單元的劃分。區域降水量根據巴林左旗、林西、翁牛特旗、寶國圖和赤峰5個降水測站點的數據插值獲取。根據全國分級水網數據，通過ArcGIS軟件密度分析工具，構建翁牛特旗水網密度。

圖1 研究區區劃及地下水測站分布圖Fig.1 Division of study area and distribution map of groundwater stations

1.3 協同克里金插值模型

如圖1所示，地下水埋深是以點數據的形式分布于研究區內，為研究整個空間范圍內地下水埋深的時空分異規律，并對其進行定性定量分析，首先需對地下水埋深數據在空間上進行呈現。因此，選取協同克里金插值，相對于其他插值方式，協同克里金插值(co-Kriging)在區域范圍內進行插值的主變量與多個屬性變量之間存在很強的空間相關性，協同克里金插值是將區域空間估計值的參考變量在一個主變量的基礎上增加了多個區域化協同變量的方法，通過模型的擬合，分析主變量與協同變量之間的相關性[12-14]，有利于更好地估算整個區域范圍內的因變量分布情況。

本文結合影響地下水埋深的靜態特征(如植被覆蓋程度、降水量和水網分布等因子)作為協變量，模擬地下水埋深在研究區范圍內的分布情況。地表植被覆蓋情況和地形的分布特點對于地下水埋深的擴散范圍、擴散方式與擴散程度有重要的影響，植被覆蓋情況能反映土壤的固水能力，植被覆蓋率高的地方普遍地下水埋深較為穩定，多年間變化不大，影響地下水埋深的時間變化特征，且輻射能力較弱，地形則影響地下水的匯流過程，是地下水擴散的主要途徑，降水是地下水徑流的主要來源，水網是地下水擴散的重要載體，二者共同決定地下水埋深的空間變化特征。因此，本文選取NDVI(歸一化植被指數)、區域降水量和水網密度作為協變量，對近10年間地下水埋深進行協同克里金插值。

協同區域化的協變量K所對應的影響因子數據{Zk(x)，k=1，2，…，x}，k0為主變量地下水埋深變化速率，在各地下水測站點上埋深變化速率Zk0(xi)(i=1,2,…,n)的數學期望存在并為一個定值mk0。

ZVk0對中心點x0在待估測區域V(x0)內區域化變量地下水埋深變化速率k0估測的平均值為

(1)

承載Vαk上的平均值為

(2)

(3)

式中μαk——協同克里金插值各因素權重系數

Zαk——協同化區域變量

協同化區域變量在地下水測站點{x1，x2，…,xn}的范圍內符合二階平穩假設，二階平穩假設是指區域化變量Z(x)的任一n維分布函數不因空間點x發生位移而改變，具有數學期望存在且平穩，方差和協方差存在且平穩的性質。在本文中協同化區域變量為降水量、河網密度和NDVI，分別用ui、vj、ps表示，即

(4)

式中ai、bj、cs——降水量、河網密度和NDVI的權重系數

1.4 最小水文響應單元

依據協同克里金插值對地下水埋深的處理結果在空間上呈現為柵格數據，但在實際情況下，因多重因素空間差異和土地利用類型的不規則化，使得地下水埋深的劃分單元不是類似于柵格的規則形狀。因此，利用劃分水文響應單元的方法，構建最小水文響應單元[15-17]，使柵格數據轉換為矢量數據，使后面的空間分析更貼合于研究區的實際情況。

最小水文響應單元基于水文響應單元(Hydrological response units，HRU)，原定義為SWAT模型中模擬的基本單元，SWAT模型中的基本單元HRU通過閾值劃分由土地利用、土壤類型和坡度共同定義。結合適應指數法，構建指數閾值體系，能在DEM的基礎上通過水文分析更好地劃分研究區的子流域。在生成子流域的基礎上，依據土地利用數據作為劃定標準，確定水文響應單元，并以此尺度開展相關研究[18-21]。地下水的來源主要有降水入滲、灌溉水入滲、地表水入滲補給、越流補給和人工補給。綜合分析，主要由降水、地表徑流和所處環境決定。地表徑流入滲是地下水補給的重要途徑，但本文不僅依據地表徑流和地形劃分的水文響應單元來作為分析單元，而是在此基礎上依據研究區的地類信息，考慮不同土壤特性，不同土地利用情況的因素，并以此構建最小響應單元。因地下水的響應單元劃分無法完全實地獲取，相對于柵格單元和單純地表匯流劃分的水文響應單元，最小響應單元更具備科學性。

流域計算根據ArcGIS軟件中的水文分析功能，利用D8單流向算法原理對DEM每個柵格單元的周圍8個方向進行定義，并假定雨水降落在地形中某一個格子上，該格子的水流將會流向周圍8個格子中地形最低的格子。如果多個像元格子的最大下降方向都相同，則會擴大相鄰像元范圍，直到找到最陡下降方向為止。由于DEM數據并不均勻分布，存在因周邊像元過高產生不合理的凹陷點而帶來基礎數據的誤差，故將不合理的凹陷點填充至合理值[22-24]。

因此，用水文響應單元來表達地下水埋深的變化規律具備可行性，將協同克里金插值結果在ArcGIS軟件中利用水文分析工具對水文響應單元進行劃定，水文響應單元依據研究區的地形進行劃分，對基于DEM劃分的流域進行二次劃分，根據DEM劃分的流域在空間內僅能代表地形一種特點，在同一流域內也存在多種地類，存蓄的地下水含量也各不相同。因此，利用土地覆蓋數據對流域進行二次劃分，并對邊界的小斑塊進行合并，得到最終的最小水文響應單元，每個最小水文響應單元的土地覆蓋與流域劃分各不相同，具備一定的獨特性，以此作為地下水埋深評價的最小水文響應單元能夠最大限度地分析地下水埋深的時空分異規律。

對劃分好的全新最小水文響應單元進行賦值，因協同克里金插值結果是在整個研究區的柵格數據，而水文響應單元是空間數據，因此，通過分析工具中的zonal模塊對插值結果統計空間范圍內各水文響應單元的均值，并通過屬性連接在空間上進行呈現和分析。

1.5 空間自相關

定量研究地下水的時空分異特征，需要將劃分好的水文響應單元和地下水埋深協克里金插值結果耦合分析，將地下水埋深數值通過水文響應單元進行呈現，對其結果進行空間自相關分析。空間自相關是指一些變量在同一個分布區內的觀測數據之間潛在的相互依賴性。空間自相關就是研究空間中，某空間單元與其周圍單元間，就某種特征值，通過統計方法，進行空間自相關性程度的計算，以分析這些空間單元在空間上分布現象的特性[25-28]。描述空間自相關的參數包括全局自相關和局部自相關，本文主要通過局部自相關來表達地下水的擴散及匯流過程。

局部空間自相關是用來度量局部空間單元對于整體研究范圍空間自相關的影響程度的相關性指數，代表特定水文響應單元與鄰近水文響應單元的地下水埋深的相關程度，計算公式為

(5)

式中m——與水文響應單元i相鄰接的水文響應單元個數

Wij——n階矩陣的元素，代表第i個單元格與第j個單元格的空間拓撲關系

xi——第i個單元格地下水埋深

xj——第j個單元格地下水埋深

結合多年地下水埋深數據，根據最小水文響應單元劃分的區域進行空間自相關分析。可得到近20年間地下水埋深的變化趨勢，并以空間地塊的方式進行呈現，每個最小水文響應單元分別具備不同的景觀類型，分析實際情況下研究區內地下水埋深分布的聚類程度和相關情況。具備空間高聚集特性的區域，地下水匯流過程較為通暢，由結果反推過程，根據明顯的自相關特性分析得出在范圍內地下水流通的過程不具備較大的阻礙。

2 結果與分析

2.1 地下水埋深協同克里金插值

基于ArcGIS軟件中的地統計向導工具，以2005、2009、2013、2017年地下水埋深為主變量，以對應年份夏季的NDVI、降水量和研究區內水網密度為協變量，進行協同克里金插值分析。在協同克里金插值中選擇不同插值方式和參數應用，確定最好的插值效果圖像。對比普通克里金插值、簡單克里金插值、泛克里金插值和析取克里金插值4種插值方法，普通克里金插值多用于單個變量的無偏估計，地下水埋深結果具有較強的聚類特點，所以普通克里金插值不適用于地下水埋深研究；泛克里金插值需要知道整個插值空間的整體變化趨勢，在本研究中，插值前并沒有地下水變化趨勢的整體估計；析取克里金插值是非線性的插值模型，在協同克里金插值結合協變量的基礎上非線性插值增加了插值結果的不確定性，故此，普通克里金插值、泛克里金插值和析取克里金插值并不適用于本文的研究，研究地表形變及其影響因子的耦合關系分析，確定簡單克里金插值效果最好。選擇不同協方差函數進行函數擬合精度分析，環形模型、球面模型、三球模型、K-貝塞爾模型、J-貝塞爾模型和穩定模型等11個函數模型的對比如表1所示，在對協方差模型的精度驗證中，將標準平均值作用選擇最優模型的標準，在11個模型中穩定模型的效果最好。

表1 不同半方差函數模型插值精度評價Tab.1 Evaluation of interpolation accuracy of different semivariance function models m

圖2 普通克里金插值結果Fig.2 Ordinary Kriging interpolation results

圖3 協同克里金插值結果Fig.3 Co-Kriging interpolation results

對比用單一的地下水埋深作為插值因素的克里金插值(圖2)，協同克里金插值結果(圖3)沒有突兀的點狀聚集特性，插值結果在空間上更平滑，能更好體現地下水在空間上的分布規律。在整個研究區范圍內，插值結果分布空間分異較為明顯，其中，東西兩側的地下水埋深差距明顯，且在所有年份中地下水埋深受土地利用類型影響很大。在4 a時間里，東側地下水埋深較大，主要分布在西拉木倫河下游以勝利村為中心的農田和杜勒本吉附近的沙丘荒地；翁牛特旗中部及東部沙丘主要是平原沙丘，主要為沙漠，零星有荒草等植被覆蓋，連年的干旱缺水環境和無植被根系固水的情況下，在這個范圍地下水埋深較低；而大片的農田地帶雖然緊鄰西拉木倫河，且存在一定的農作物覆蓋，但在農作物密集種植的情況下，其需水量巨大，當地居民普遍在農田周圍附近打井就近灌溉，每年巨大的需水量導致在大范圍農田分布區地下水埋深也普遍較深，老哈河下游農田集中區也與其類似，河流下游的農業開發活動導致地下水埋深情況不佳。地下水埋深較好的區域主要分布在西部山區，綜合4 a插值結果來看，西部山區植被覆蓋程度較高，高大喬灌木較多，也更為密集，植物的多樣化也有利于土壤環境的維持與改善，發達的根系和良好的土壤環境有利于土壤水分的維持，每次降雨有大量的水分被植物固定在土壤中，因此，研究區東部普遍地下水埋深較淺；此外，河流附近的地下水埋深與周圍相比也呈現出較好的特點，其中尤以研究區北部西拉木倫河和南部羊腸子河、紅山水庫最為明顯，地下水埋深都是通過就近匯入河流進行擴散，在一些大型河流的附近由于小河流的匯入和水分的擴散，致使地下水埋深較高。

對比2005、2009、2013、2017年的地下水埋深插值結果，以4 a為界限，分析地下水埋深在時間上的變化規律。2005年地下水埋深在2.731 67～4.003 64 m之間，整體地下水埋深分布較為均勻，地下水埋深較大值主要分布在研究區的東北部沙丘地帶，且研究區內大部分區域地下水埋深在3 m以上，說明當年地下水擴散和匯流較為順暢，沒有太多阻礙因素，地下水埋深低值沿河流分布特征明顯。2009年地下水埋深在2.603 1～4.107 38 m之間，2009年地下水空間分異不明顯，所有地下水埋深高值區域集中分布在東北部和南部農田范圍內，且地下水埋深達到3 m以上，說明該年地下水埋深受環境影響較大，且降水較為稀少，河流區與農田區地下水埋深較為接近，農田需要過多抽取地下水進行農業灌溉。2013年地下水埋深在3.086 97～4.643 16 m之間，地下水埋深差異性很小，地下水在土壤中擴散與匯流過程通暢，不受太多阻礙，但整體地下水埋深較低，整體分布趨勢與2009年類似，降水較為平均但不夠充沛，地面徑流起到主要調節作用。2017年地下水埋深在2.324 57～4.520 02 m之間，地下水埋深在研究區內分布較為均勻，在西拉木倫河、老哈河和羊腸子河附近地下水埋深較低，河流對于地下水的匯集有突出作用。綜合以上4 a結果可得出，短期內地下水埋深的平均值不會有太大波動，但局部地區，如：東北部農田、中東部沙丘等，變化明顯，且地下水埋深受當年降水等因素的影響較大，在時間序列上，地下水埋深受河流影響越來越大，地下水在土壤中的匯流過程起到了越來越重要的作用。因此，對翁牛特旗地下水資源的開發、利用與保護應充分利用這一點，以改善水資源現狀，維持水環境的平衡。

2.2 最小響應單元構建

依據ArcGIS軟件中的水文分析，根據DEM劃分得到水文響應單元6 614個，如圖4所示，但傳統的水文響應單元在邊界處存在大量的細碎小斑塊，通過消除工具對小于75 hm2的細碎小斑塊進行消除，得到水文響應單元83個。但此水文響應單元僅僅是根據研究區內地形劃分的水文響應單元，地下水埋深在空間上的分區劃分還受土地利用類型的影響，例如耕地、林地和建筑用地對于地下水擴散和匯流的影響各不相同，不同的土地類型具備不同的土壤類型、植被覆蓋程度，土壤含水率和用水程度等多種特征的差別。因此，將依據地形劃分的水文響應單元結合土地利用情況進行二次劃分，劃分后進行小斑消除，最終得到最小響應單元551個。最小響應單元代表不同的土地利用類型和地下水匯流方向，以此為載體在空間上呈現地下水埋深的擴散和匯流情況具有重要意義。

圖4 水文響應單元分布圖Fig.4 Distribution map of hydrological response units

研究區內最小響應單元的分布如圖5所示，不同類型的最小響應單元具備一定的分布規律。其中林地的最小響應單元50個，均分布在翁牛特旗的西部和西南部，并存在明顯的以國道為界限，國道以西以林地為主，且高程普遍在1 000 m以上，林地最小響應單元分布較為集中，且面積較大，在林地范圍內，最小響應單元的聚集特點突出，范圍內的地下水埋深擴散與匯流特性較為一致。研究區內水體分布也較為集中，最小響應單元主要分布在北部西拉木倫河流域以及老哈河下游的紅山水庫，共26個，面積均較小。耕地最小響應單元158個，面積普遍不大且分布較為零散，具體特點不顯著，分布在翁牛特旗東北角，西拉木倫河下游平原區、西部林地中河流兩側、中部林草-沙漠交界帶以及老哈河下游沖擊平原。建筑用地最小響應單元1個，位于翁牛特旗市中心。草地最小響應單元230個，為6種地類中最多的，也是分布最為零散的一類最小響應單元，其特點是環繞沙漠地區分布。裸地最小響應單元86個，主要分布在研究區的中部和東部，其中有5塊面積較大，較為集中，地下水在這類最小響應單元中擴散的規律最為明顯。

圖5 最小響應單元分布圖Fig.5 Distribution map of minimum response unit

2.3 空間自相關分析

協同克里金插值的結果在空間上呈現為柵格數據，完成最小響應單元構建后，將協同克里金插值結果賦給最小響應單元，再利用空間統計工具中的空間自相關下的局部自相關分析4 a地下水埋深的分布及演變規律。空間自相關結果如圖6所示。

圖6 地下水埋深空間自相關特征Fig.6 Spatial autocorrelation characteristics of groundwater depth

由圖6可知，研究區東北部的135個最小響應單元的地下水埋深結果在2005—2017年均呈現出高聚集特點，其大部分為耕地，說明在長時間序列上，耕地的地下水埋深均處在一個水位較低的區間內，且聚集特點明顯，說明在大片農田范圍內，地下水的擴散機制較為通暢，農作物的需水量使得這部分連片區域呈現為高值聚集區；研究區中東部的沙丘平原區也呈現出大范圍的高值聚集區，在個別年份呈現出非較大值聚集特點，說明平原沙丘區的地下水埋深空間分布很不規律，不夠穩定，空間內地下水埋深較低，但是存在個別斑塊地下水埋深較高，地下水交互存在障礙，聚集特征不如耕地明顯；研究區西部有114個最小響應單元在4 a中均為低值聚集區，西部依據林地和地形劃分的最小響應單元偏大，同樣面積內，低值聚集區的最小響應單元個數會少很多，由低值聚集區的分布可明顯發現其主要分布于林地的范圍內，而不是在水體等地下水埋深較高的最小響應單元，這說明西部山區范圍內地下水埋深呈現大范圍的聚集，森林中樹木發達的根系能夠促進地下水的產流與匯流。對比4 a地下水埋深的空間分異特性，研究區西部山區林地最小響應單元內的低值聚集區變化差異不大，但在微弱的年際變化間，低值聚集區呈現明顯的聚集特點，部分低值次聚集區和低值分散區逐漸減少甚至消失，致使研究區西部的聚集特點逐年增加；在研究區東部，2005—2017年整體分布和變化規律與西部類似，但除2013年外東部的變化更加微弱，耕地的聚集特性較為穩定。這說明近年來整體地下水埋深呈現聚集特點，地下水的匯流過程呈現一種逐漸有序的狀態。

3 結論

(1)以2005—2017年4期地下水埋深為主變量，以NDVI、降水量和河網密度等在空間上影響地下水埋深的數據作為協變量，采用協同克里金插值方法將地下水埋深數據在空間呈現。在4 a插值結果中，地下水埋深平均變化不大，地下水埋深整體在2～5 m之間浮動，但在中東部沙丘區變化較為明顯，整體呈現出向河流方向匯集的趨勢。

(2)運用水文響應單元模式，在空間上將協同克里金插值后的地下水埋深柵格數據轉換為矢量數據。通過水文分析劃定水文響應單元，結合研究區土地利用類型，對其重新劃分為551個最小水文響應單元。在對最小水文響應單元賦值的基礎上，進行空間自相關分析，4 a數據中，持續高值聚集區135個，低值聚集區114個，整體地下水埋深呈現聚集特點，地下水的匯流過程呈現一種逐漸有序的狀態。

(3)對研究區內地下水埋深時空分異特性分析表明，在空間上，地下水埋深東西分異規律明顯，這主要源于西部山區森林帶、東部農田和沙丘帶的土地利用以及地形地貌的差異，地下水位呈現西高、東低的態勢，且地下水埋深情況較好的區域逐漸沿翁牛特旗的主要河流分布，受降水和河流影響逐漸變大；在時間序列上，地下水埋深平均變化不大，但逐漸呈現聚集的趨勢，這說明阻礙地下水擴散的障礙正在逐漸消除。因此，建議對于翁牛特旗地下水資源的開發與利用應結合以上特點，合理規劃、重點開發、有效利用。