基于穩(wěn)健回歸-去趨勢(shì)波動(dòng)分析法的山前平原地下水轉(zhuǎn)換關(guān)系研究

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

(1.中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院，湖北 武漢 430074；2.中國(guó)地質(zhì)調(diào)查局武漢地質(zhì)調(diào)查中心，湖北 武漢 430205)

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

地下水動(dòng)態(tài)變化特征的研究是目前地下水資源研究中具有挑戰(zhàn)性的領(lǐng)域和熱點(diǎn)之一[1-3]。地下水動(dòng)態(tài)既反映出地區(qū)地下水資源的相關(guān)信息，也可對(duì)地下水之間的相互轉(zhuǎn)換關(guān)系進(jìn)行指示。地下水在補(bǔ)徑排過程中，水位發(fā)生波動(dòng)，水位過程線產(chǎn)生分形，由于地表、土壤、含水層等的阻滯作用，使各含水層中地下水水位時(shí)間序列的標(biāo)度指數(shù)發(fā)生變化[4]，因此可通過對(duì)不同類型含水層中地下水的分形標(biāo)度進(jìn)行量化，來分析探討不同類型含水層中地下水之間的轉(zhuǎn)換關(guān)系。

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

1 研究方法

穩(wěn)健回歸-去趨勢(shì)波動(dòng)(robust Detrended Fluctuation Analysis,r-DFA)方法是包括DFA和統(tǒng)計(jì)模型的過程，能夠?qū)π盘?hào)進(jìn)行單分形量化，并自動(dòng)確定交叉點(diǎn)位置。本文使用r-DFA方法的代碼由Habib于2016年在Mathworks網(wǎng)站上提供(文件ID為60026)。

1.1 去趨勢(shì)波動(dòng)分析(DFA)過程

DFA法最早由Peng等[11]分析DNA相關(guān)性時(shí)提出，其計(jì)算步驟如下：

(1)

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

(2)

(4) 計(jì)算長(zhǎng)度L中所有段的均方根F(L)：

(3)

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

F(L)≈Lα

(4)

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

1.2 確定交叉點(diǎn)

全局標(biāo)度指數(shù)α的變化指示出時(shí)間序列出現(xiàn)單分形行為，其表現(xiàn)為F(L)與L雙對(duì)數(shù)關(guān)系圖的斜率發(fā)生變化，發(fā)生變化的時(shí)間(L)為交叉點(diǎn)[4]。

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

1.3 去趨勢(shì)波動(dòng)分析(DFA)修正

(5)

2 研究區(qū)概況與數(shù)據(jù)來源

2.1 研究區(qū)水文地質(zhì)概況

研究區(qū)位于云應(yīng)盆地北部澴河河谷階地，屬于大別山山前平原區(qū)，氣候?yàn)閬啛釒Т箨懶詺夂?，雨熱同期，地?shì)北高南低，澴河從北向南縱穿，沿河兩側(cè)發(fā)育河漫灘與多級(jí)階地,見圖1。

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

2.2 數(shù)據(jù)來源

圖1 研究區(qū)水文地質(zhì)平面簡(jiǎn)圖及地下水水位長(zhǎng)期觀測(cè)孔分布

圖2 澴河?xùn)|側(cè)水文地質(zhì)剖面示意圖

表1 研究區(qū)地下水水位長(zhǎng)期觀測(cè)孔基本信息統(tǒng)計(jì)表

地下水水位動(dòng)態(tài)監(jiān)測(cè)數(shù)據(jù)全部采用Levelogger solinst自動(dòng)水位計(jì)進(jìn)行采集，采集頻率為1次/4 h；同時(shí)，收集了國(guó)家氣象站孝感站2017年5月—2019年2月的4小時(shí)降雨量數(shù)據(jù)。

2.3 地下水水位動(dòng)態(tài)數(shù)據(jù)的變化特征

研究區(qū)2017年5月—2019年2月降雨量和各類型含水層中地下水水位的動(dòng)態(tài)監(jiān)測(cè)數(shù)據(jù)，見圖3。

圖3 研究區(qū)2017年5月—2019年2月降雨量和各類型含水層中地下水水位的動(dòng)態(tài)監(jiān)測(cè)數(shù)據(jù)

由圖3可以看出：

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

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

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

3 結(jié)果與分析

3.1 不同階數(shù)r-DFA結(jié)果比較

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

圖4 研究區(qū)降雨r-DFA1～r-DFA6的F(n)(L)與L雙對(duì)數(shù)關(guān)系圖

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

3.2 降雨和地下水水位r-DFA3標(biāo)度指數(shù)結(jié)果

研究區(qū)降雨r-DFA3的F(L)與L雙對(duì)數(shù)關(guān)系圖，見圖5。

圖5 研究區(qū)降雨r-DFA3的F(L)與L雙對(duì)數(shù)關(guān)系圖

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

研究區(qū)不同類型含水層中地下水水位r-DFA3的F(L)與L雙對(duì)數(shù)關(guān)系圖，見圖6。

由圖6可以看出：

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

圖6 研究區(qū)各觀測(cè)孔地下水水位r-DFA3的F(L)與L雙對(duì)數(shù)關(guān)系圖

3.3 結(jié)果分析

研究區(qū)降雨和不同類型含水層中地下水水位r-DFA3的計(jì)算結(jié)果匯總于表2。

研究區(qū)降雨的全局標(biāo)度指數(shù)α值為0.607 6，與Li等[4]和Matsoukas等[7]的研究結(jié)果一致。Matsoukas等研究了美國(guó)9個(gè)地區(qū)的降雨序列，在5～10 d左右均出現(xiàn)交叉，推測(cè)交叉點(diǎn)的出現(xiàn)與氣象和氣候系統(tǒng)的分離有關(guān)。

表2 研究區(qū)降雨和不同類型含水層中地下水水位的r-DFA3計(jì)算結(jié)果匯總

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

在所有的水文信號(hào)中，降雨的波動(dòng)最大，在降雨入滲的過程中，由于地面、包氣帶和含水層的阻滯作用，將時(shí)間序列中較為復(fù)雜的波動(dòng)濾掉，留下較為穩(wěn)定的連續(xù)波動(dòng)，其標(biāo)度指數(shù)增大，通過比較標(biāo)度指數(shù)，可體現(xiàn)含水層性質(zhì)以及各類型含水層中地下水之間的轉(zhuǎn)換關(guān)系。

研究區(qū)Ey砂巖含水層中觀測(cè)孔ZK08與ZK04-1在第一交叉點(diǎn)前地下水水位波動(dòng)趨勢(shì)為非持續(xù)性，點(diǎn)后為持續(xù)性，第二交叉點(diǎn)后地下水水位波動(dòng)趨勢(shì)為非持續(xù)性，但兩交叉點(diǎn)橫坐標(biāo)不同，且ZK08觀測(cè)孔無(wú)第三交叉點(diǎn)，表明地下水分形具有場(chǎng)地特異性。

研究區(qū)QbW變質(zhì)巖含水層的滲透系數(shù)在0.016～0.14 m/d之間，滲透性差且受到人為活動(dòng)干擾最小，其地下水水位標(biāo)度指數(shù)最大。

4 結(jié) 論

本文結(jié)合研究區(qū)降雨量和地下水水位數(shù)據(jù)的動(dòng)態(tài)特征及r-DFA分析結(jié)果，得到如下結(jié)論：

(1) 研究區(qū)降雨波動(dòng)類似于白噪聲，地下水水位波動(dòng)更傾向于分形布朗運(yùn)動(dòng)，且地下水分形存在場(chǎng)地特異性。

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