基于穩健回歸-去趨勢波動分析法的山前平原地下水轉換關系研究

劉曉芮，王 清，陳植華，胡 成

(1.中國地質大學(武漢)環境學院，湖北 武漢 430074；2.中國地質調查局武漢地質調查中心，湖北 武漢 430205)

山前平原是山區向平原的過渡地段，既包括基巖山區，也包括平原河谷階地。大別山-云應盆地山前平原過渡區內，東北部地勢最高處為中元古界青白口系變質巖，向西南地勢降低，為新生界古近系砂巖與變質巖不整合接觸，伏于第四系松散堆積物之下，是研究區主要開采含水層。由于區域內第四系上更新統上部沉積物為相對隔水的亞黏土，其孔隙含水層與砂巖含水層相對封閉，主要開采層的地下水補給來源不明確，各類型含水層中地下水的轉換關系缺乏較深入的研究。因此，查明區域內含水層中地下水之間的相互轉換關系對地下水資源的合理開發和保護具有重要意義。

地下水動態變化特征的研究是目前地下水資源研究中具有挑戰性的領域和熱點之一[1-3]。地下水動態既反映出地區地下水資源的相關信息，也可對地下水之間的相互轉換關系進行指示。地下水在補徑排過程中，水位發生波動，水位過程線產生分形，由于地表、土壤、含水層等的阻滯作用，使各含水層中地下水水位時間序列的標度指數發生變化[4]，因此可通過對不同類型含水層中地下水的分形標度進行量化，來分析探討不同類型含水層中地下水之間的轉換關系。

對于時間序列的分形標度行為的分析方法眾多，去趨勢波動分析(Detrended Fluctuation Analysis,DFA)方法能夠有效檢測嵌入非平穩時間序列中的長程相關性，避免由外部影響產生的虛假效果，被認為是一種可靠的非平穩信號的分析方法[5-7]，在降雨、河流等水文時間序列分形方面的應用較為廣泛，但在地下水時間序列分形方面的應用較少。Tu等[8]運用DFA法對承壓含水層中地下水水位的單分形進行了量化，結果表明地下水水位時間序列存在分形，得到的地下水水位標度指數與時間序列的長度和時間間隔相關；Li等[4]運用DFA法對降雨-徑流-地下水-基流的標度指數進行分析，認為地下水水位波動為分形布朗運動，由于土壤及含水層的抑制作用，降雨-地下水-基流的標度指數逐漸增大;Yu等[9]運用DFA法對波多黎不同類型含水層中地下水分形尺度進行了研究，結果表明沖積河谷含水層中地下水單分形指數高于巖溶含水層中地下水單分形指數。運用DFA法對單分形信號進行分析時，由于交叉點的確定是主觀的，故會影響單分形量化的可靠性。為了克服此缺點并使量化過程自動化，Habib等[10]提出了穩健回歸-去趨勢波動分析(r-DFA)方法，并運用該方法對英國瓦林福德地區地下水及河流水位、地下水及河流水溫、降雨和氣溫的分形尺度進行了量化，客觀地確定交叉點，使其具有統計學意義，并在時間域和分形域對該地區降雨、空氣、地下水、河流之間的相互聯系進行分析。

1 研究方法

穩健回歸-去趨勢波動(robust Detrended Fluctuation Analysis,r-DFA)方法是包括DFA和統計模型的過程，能夠對信號進行單分形量化，并自動確定交叉點位置。本文使用r-DFA方法的代碼由Habib于2016年在Mathworks網站上提供(文件ID為60026)。

1.1 去趨勢波動分析(DFA)過程

DFA法最早由Peng等[11]分析DNA相關性時提出，其計算步驟如下：

(1)

(2) 將Y(ti)分成m個時間長短為L的非重疊段，使m=int(N/L)，每一段累積離差都標注為Yj,k，其中j表示分段長度，j= 1,2,…,L；k表示分段個數，k=1,2,…,m。

(2)

(4) 計算長度L中所有段的均方根F(L)：

(3)

(5) 對不同的L值重復步驟(1)至(4)，然后在對數軸上繪制F(L)與L的關系圖，以確定標度指數α，即最佳擬合直線的斜率：

F(L)≈Lα

(4)

采用穩健回歸(bisquare估計)線性擬合得到的α為r-DFA的全局標度指數，是忽略標度指數的任何局部變化而確定的直線斜率，可以確保其基于預定義范圍內的殘差。當α=0.5時，時間序列為白噪聲，無長程相關性；α=1.5時，時間序列為分形布朗運動；當1.5<α<2時，一個時間序列是持續性的；當α<1.5時，時間序列是非持續性的。持續性表明，如果時間序列在一段時間內持續增加(或減少)，則在類似的一段時間內，時間序列將繼續增加(或減少)；反之，非持續性則表明時間序列在類似時間內，時間序列與原趨勢相反。

1.2 確定交叉點

全局標度指數α的變化指示出時間序列出現單分形行為，其表現為F(L)與L雙對數關系圖的斜率發生變化，發生變化的時間(L)為交叉點[4]。

在通過穩健回歸擬合確定α后，還需要進行分段回歸，通過最小化數據與擬合直線之間的最小二乘誤差來優化交叉點位置；然后通過對DFA結果進行協方差分析(ANCOVA)和多重比較，確定擬合數據交叉點的數量，保證在95%顯著性水平的基礎上能產生顯著差異的斜率的最大交叉點；將交叉點的數量逐漸增加，直到該方法不能產生任何進一步的顯著標度指數。

1.3 去趨勢波動分析(DFA)修正

(5)

2 研究區概況與數據來源

2.1 研究區水文地質概況

研究區位于云應盆地北部澴河河谷階地，屬于大別山山前平原區，氣候為亞熱帶大陸性氣候，雨熱同期，地勢北高南低，澴河從北向南縱穿，沿河兩側發育河漫灘與多級階地,見圖1。

受地形地貌及含水層結構特征的影響，研究區地下水接受大氣降水補給后，總體上由北往南、由東西兩側往中部澴水方向徑流及排泄(見圖1)。

2.2 數據來源

圖1 研究區水文地質平面簡圖及地下水水位長期觀測孔分布

圖2 澴河東側水文地質剖面示意圖

表1 研究區地下水水位長期觀測孔基本信息統計表

地下水水位動態監測數據全部采用Levelogger solinst自動水位計進行采集，采集頻率為1次/4 h；同時，收集了國家氣象站孝感站2017年5月—2019年2月的4小時降雨量數據。

2.3 地下水水位動態數據的變化特征

研究區2017年5月—2019年2月降雨量和各類型含水層中地下水水位的動態監測數據，見圖3。

圖3 研究區2017年5月—2019年2月降雨量和各類型含水層中地下水水位的動態監測數據

由圖3可以看出：

(1) 每次降雨前后，研究區西北部Ey砂巖孔隙-裂隙水含水層中地下水水位上升約0.05～0.1 m，波動微小；豐水期，Ey砂巖含水層中地下水水位于降雨增加約30 d后呈上升趨勢，地下水水位上升幅度約0.5 m，持續上升約一個月后緩慢下降，地下水水位下降約0.7 m[見圖3(c)]。結果表明:在降雨后研究區北部地區地下水得到補給，并向南部Ey砂巖含水層徑流；枯水期北部補給區地下水接受降雨補給減少，Ey砂巖含水層中地下水水位開始緩慢下降。

(2) QbW變質巖類風化裂隙水含水層的富水性弱，對降雨的響應快，每次降雨后1 d內QbW變質巖含水層中地下水水位上升約0.3～0.5 m；2018年7—10月受到人為開采的影響，QbW變質巖含水層中地下水水位持續下降約2 m[見圖3(b)]；在豐水期，QbW變質巖含水層受到降雨持續補給，含水層中地下水水位上升約0.5 m，從峰值變化滯后降雨約30 d；在枯水期，地下水水位下降約為1 m?？傮w上看，青白口系變質巖含水層對降雨較為敏感，雨水轉化為QbW變質巖類風化裂隙水。

在對昌馬灌區近十年地下水水位年內變化分析時，分別以8眼監測井地下水水位2002—2010年期間逐月平均值（見圖2）和2002年與 2010年年內地下水水位埋深作為指標進行分析，如圖5所示。從圖5中可以看出，以8眼監測井地下水水位年平均埋深作為指標來看，在研究時段內，Ⅰ-1、Ⅰ-2和Ⅰ-3監測井的水位年內隨季節變化幅度較大，而其余監測井地下水水位年內隨季節變化幅度不顯著。同時，可以看出，以2002年和2010年為典型年的年內變化研究中可以看出，除Ⅰ-1、Ⅱ-1和Ⅱ-2監測井外，在兩個典型年內，其余監測井地下水水位埋深的年內變化基本相同，且與多年年內地下水水位變化趨勢基本一致。

3 結果與分析

3.1 不同階數r-DFA結果比較

本文先對降雨量和6個地下水水位觀測點的地下水水位動態數據進行DFA計算，再利用公式(5)來解決分形行為中的偏差，打亂順序重新配置的次數選擇為100，得到r-DFA1～r-DFA6，并采用穩健回歸來確定所有階的DFA的全局標度指數α；然后使用分段線性回歸確定交叉點，并使用ANCOVA分析結果，再通過多重比較程序評估，最終得到研究區降雨r-DFA1～r-DFA6的F(n)(L)與L雙對數關系圖，見圖4。

圖4 研究區降雨r-DFA1～r-DFA6的F(n)(L)與L雙對數關系圖

由圖4可見，經過修正后的高階r-DFA擬合直線偏離共線趨勢較小，與低階r-DFA1、r-DFA2擬合結果相近；在高階r-DFA擬合結果中，r-DFA4和r-DFA6原始數據與擬合直線在小時間尺度上擬合較差，r-DFA5擬合直線與其他階擬合直線存在偏差，出現負斜率線段，而r-DFA3濾去了r-DFA1的線性趨勢和r-DFA2的二階趨勢，與直線擬合最為接近，具有典型性和代表性，因此本文選擇r-DFA3的擬合結果對斜率α值和各分段點斜率及交叉點橫坐標進行分析。

3.2 降雨和地下水水位r-DFA3標度指數結果

研究區降雨r-DFA3的F(L)與L雙對數關系圖，見圖5。

圖5 研究區降雨r-DFA3的F(L)與L雙對數關系圖

由圖5可見，研究區降雨的全局標度指數α值為0.607 6，類似于α=0.5的白噪聲；研究區降雨序列在4.8 d左右發生交叉，當時間從小尺度上升到中等尺度時，降雨的α值從0.747 2下降到0.600 7。

研究區不同類型含水層中地下水水位r-DFA3的F(L)與L雙對數關系圖，見圖6。

由圖6可以看出：

(2) ZK06觀測孔的QbW變質巖風化裂隙水含水層中地下水水位的α值為1.614 7，在19.5 d發生交叉[見圖6(b)]；在小時間尺度上，QbW變質巖含水層中地下水水位波動趨勢是持續性的(α=1.641 3)，在中等時間尺度上，其波動趨勢是非持續性的(α=1.461 8)。

圖6 研究區各觀測孔地下水水位r-DFA3的F(L)與L雙對數關系圖

3.3 結果分析

研究區降雨和不同類型含水層中地下水水位r-DFA3的計算結果匯總于表2。

研究區降雨的全局標度指數α值為0.607 6，與Li等[4]和Matsoukas等[7]的研究結果一致。Matsoukas等研究了美國9個地區的降雨序列，在5～10 d左右均出現交叉，推測交叉點的出現與氣象和氣候系統的分離有關。

表2 研究區降雨和不同類型含水層中地下水水位的r-DFA3計算結果匯總

研究區不同類型含水層中地下水水位的全局標度指數α值表明，地下水水位波動更傾向于分形布朗運動(α=1.5)，與前人的研究結果相同[4,8-9]。Li等[4]報道了美國愛荷華州地下水水位分別在6～30 d和50～365 d之間存在交叉點，其6～30 d的交叉結果與本文在較小時間尺度(4～10 d)的交叉結果相似，表明本研究區地下水分形特征與前人的研究結果一致，具有較高的可信度。

在所有的水文信號中，降雨的波動最大，在降雨入滲的過程中，由于地面、包氣帶和含水層的阻滯作用，將時間序列中較為復雜的波動濾掉，留下較為穩定的連續波動，其標度指數增大，通過比較標度指數，可體現含水層性質以及各類型含水層中地下水之間的轉換關系。

研究區Ey砂巖含水層中觀測孔ZK08與ZK04-1在第一交叉點前地下水水位波動趨勢為非持續性，點后為持續性，第二交叉點后地下水水位波動趨勢為非持續性，但兩交叉點橫坐標不同，且ZK08觀測孔無第三交叉點，表明地下水分形具有場地特異性。

研究區QbW變質巖含水層的滲透系數在0.016～0.14 m/d之間，滲透性差且受到人為活動干擾最小，其地下水水位標度指數最大。

4 結 論

本文結合研究區降雨量和地下水水位數據的動態特征及r-DFA分析結果，得到如下結論：

(1) 研究區降雨波動類似于白噪聲，地下水水位波動更傾向于分形布朗運動，且地下水分形存在場地特異性。

(3) r-DFA法通過對地下水水位的分形研究，得到其全局標度指數α，可定量地研究不同類型含水層中地下水之間的相互轉換關系，本文的研究結果與定性分析結果一致，表明r-DFA法在分析各類型含水層中地下水之間的轉換關系具有可行性和廣闊的應用前景。