一種水文序列變異診斷耦合模型及其應用

錢龍霞，李漢霖，汪 騰，王藝瑋，王紅瑞

(1.南京郵電大學 理學院，江蘇 南京 210023；2.中國水利水電科學研究院 流域水循環(huán)模擬與調(diào)控國家重點實驗室，北京 100038；3.北京師范大學 水科學研究院，北京 100875)

根據(jù)IPCC（政府間氣候變化專門委員會）第四次及第五次報告[1–3]，全球氣候急劇變化，且變暖的趨勢不斷加強，人類活動是這一變化的重要因素。在人類活動和自然變異的雙重因素影響下，世界各地水文數(shù)據(jù)的統(tǒng)計規(guī)律和混沌程度已經(jīng)發(fā)生了顯著變化[4–5]。水文序列變異診斷對防洪減災、生態(tài)環(huán)境保護具有重要的科學意義和應用價值[6–8]。

目前，水文變異診斷方法主要用于均值變異的檢測，包括Mann–Kendall法、Pettitt檢驗法[9–10]、Spearman秩相關檢驗法[11]、BG（啟發(fā)式分割）算法等。國外相關的研究進展如下：蘭甜[1]把常用的均值變異診斷方法劃分成了5種，其中參數(shù)檢驗方法（如T檢驗和F檢驗）及非參數(shù)檢驗方法（如Mann–Kendall法和秩次相關分析法），是國內(nèi)外學者采用較多的兩類方法。參數(shù)檢驗方法對數(shù)據(jù)序列要求較高，使得非參數(shù)檢驗方法應用較多。Bernaola–Galván等[12]首次提出啟發(fā)式分割算法（又稱BG算法），適用于非平穩(wěn)、非線性序列的均值變異診斷。此后，F(xiàn)ukuda等[13]指出BG算法對于不同統(tǒng)計特性的理想樣本均具有良好的效果?？偟膩碚f，BG方法劃分得到的各子序列能在不同的物理背景下表征其特性，避免了傳統(tǒng)檢測方法對數(shù)據(jù)序列平穩(wěn)性和線性的嚴苛要求。國內(nèi)相關的應用研究有：王生雄等[14]采用Mann–Kendall法對渭河華縣站1956—2000年徑流序列進行了變異診斷分析；張蔚等[15]采用Pettitt檢驗法和Mann–Kendall法對20世紀50年代以來珠江三角洲河網(wǎng)內(nèi)17個主要控制水文站的年際潮差變化趨勢進行了分析，所得結果與實際情況符合較好；朱恒峰等[16]采用Spearman秩相關檢驗法結合有序聚類分析法，推估1960—1980年代末延河流域輸沙時變過程的顯著干擾點，其結論與現(xiàn)實擬合得很好；劉秀花等[17]以Spearman秩相關系數(shù)衡量黃河寧夏段1996—2004年污染變化趨勢，得到了污染程度下降的結論。雖然Mann–Kendall法、Pettitt檢驗法和Spearman秩相關檢驗法，在水文序列均值變異診斷中效果較好，但均存在明顯的不足，如：Mann–Kendall法檢測時會出現(xiàn)虛假變異點、Pettitt檢驗法和Spearman秩相關檢驗法僅適用于單點變異情形[18–20]。實際的水文序列可能存在多個變異點[21–22]，BG算法能夠檢測出水文序列中存在的所有變異點，如陳廣才等[20]采用BG算法對珠江三角洲順德河網(wǎng)三水站1900—2004年的年最大洪峰流量序列進行了變異診斷，分析了BG算法在水文變異診斷的可行性，并與傳統(tǒng)水文變異檢測方法進行對比，結果顯示BG算法對均值變異的敏感性和穩(wěn)定性較好。

然而，除均值變異之外，水文序列還存在動力學結構變異，傳統(tǒng)均值變異診斷方法無法刻畫這種變異特征。由于水文系統(tǒng)是一個動態(tài)的非線性復雜系統(tǒng)，通過熵值的演變規(guī)律可以揭示水文要素序列演變的內(nèi)在機理。何文平等[23]提出基于近似熵的變異檢測新方法，即滑動移除近似熵法，能夠找到動力學結構的變異區(qū)間。該方法的本質(zhì)是基于近似熵曲線進行人為判斷后得到變異診斷結果，具有很強的主觀性，而且滑動步長（即窗口長度）L在一定范圍內(nèi)對近似熵曲線的形成有著重要影響，主要存在以下問題：1）較大的L易于凸顯變異點附近的差異，并減弱整體的波動性，雖然增加了主觀判斷的可操作性，但是較大的L使得變異區(qū)間變大，變異診斷結果不精確；2）較小的L使得變異點附近的差異變小，并使整體波動性增強，增加了主觀判斷的難度。因此，如何對水文序列的內(nèi)在演變機理進行客觀診斷并給出相應的參數(shù)范圍是水文序列變異診斷模型需要解決的關鍵問題。

綜上，本文提出一種基于BG和近似熵的水文變異診斷耦合模型。首先，利用近似熵和數(shù)據(jù)滑動技術揭示水文序列內(nèi)在演變機理；再利用BG算法對近似熵序列進行分割，通過置信度檢驗檢測出水文序列存在的所有變異點；最后，通過敏感性實驗給出耦合模型參數(shù)的最優(yōu)取值范圍。

1 水文序列變異診斷耦合模型

1.1 BG算法

目前，水文序列診斷通常采用傳統(tǒng)的統(tǒng)計檢驗方法，如滑動T檢驗、滑動F檢驗、游程檢驗、秩和檢驗等，上述方法都有較為嚴格的假設，即要求序列具有線性和平穩(wěn)性的特征。而水文序列普遍具有非線性，同時，由于人類活動和自然變異的雙重影響，水文序列一般還具有非平穩(wěn)性。所以，傳統(tǒng)的方法對變異水文序列的檢測能力有限，檢測結果存在一定的偏差。

對序列的假設較為寬松，能準確檢測非線性、非平穩(wěn)的均值變異點。對于一個由N個點構成的時間序列X(t)，BG算法進行變異診斷的建模步驟如下[12]：

第1步，從左到右，依次計算每個點左邊部分和右邊部分的平均值和標準差，分別是： μ1(i)、s1(i)和μ2(i)、s2(i)。

于是，i點的合并偏差sD(i)為：

式中，N1和N2分別為i點左邊和右邊部分的點數(shù)，i=1,2,···,N。

第2步，用T檢驗的統(tǒng)計值T(i) 量化i點左右兩部分均值的差異：

對X(t)中的每一個點重復上述計算過程，得到與X(t)一 一對應的檢驗統(tǒng)計值序列T(t)，T越大，表示該點左右兩部分的均值相差越大。

第3步，計算T(t) 中 的最大值Tmax的統(tǒng)計顯著性P(Tmax)：式中，P(Tmax) 表示在隨機過程中取到T值小于等于Tmax的概率。一般情況下P(Tmax)可近似表示為：

第2步，計算任意向量X(i) 與其余向量X(j)之間的相對歐式距離d[X(i),X(j)]。

第3步，設定容許偏差為r，統(tǒng)計每個向量X(i)的d[X(i),X(j)]小 于r的數(shù)目，求出該數(shù)目與向量總數(shù)N?m+1之 比(r)。

第5步，維數(shù)m增加1，重復第1～4步，得+1(r)和?m+1(r)。

第6步，近似熵的估計值定義為ApEn(m,r)=?m(r)??m+1(r)。 顯然，ApEn值與m、r的取值有關，Pincus[25]建議r取0 .10σ ～0 .20σ （σ是原始序列的標準偏差）。本文中，取參數(shù)m=2 ，最初默認取r=0.15σ。

由蒙特卡洛模擬實驗得η=4.19lnN?11.54,δ=0.40,v=N?2，N為時間序列X(t) 總 長。Iv/(v+Tm2ax)(δv,δ)為不完全β函數(shù)，其中下標v/(v+ax) 為不完全β函數(shù)的上限。

第4步，設定置信度P0和最小分割尺度l0， 對T(t)序列的Tmax進行診斷。若其概率值大于P0，則于該點將X(t)分割成兩段均值有一定差異的子序列；否則，不分割。

第5步，對新得到的兩個子序列分別重復上述4個步驟，新子序列滿足分割條件的同時，如果子序列與其左右相鄰的子序列間的均值差異程度均滿足上述條件，則繼續(xù)分割；否則，不分割。如此重復直至所有的子序列都不可分割為止。此外，為保證統(tǒng)計的有效性，當子序列的長度小于等于l0（l0為最小分割尺度）時，不再對其進行分割。通過上述操作，作者將原序列分割為若干個不同均值的子序列，每個分割點即為均值變異點。一般而言，l0≥25,P0∈[0.50,0.95]，本文中l(wèi)0取滑動移除近似熵法處理后的數(shù)據(jù)樣本量的一半，P0取0.95。根據(jù)文獻[13]可知，l0和P0的取值，對檢測結果的影響較小。在解決不同樣本量的問題時，通過對兩者進行調(diào)整，可對序列實現(xiàn)不同尺度的變異檢測。

1.2 近似熵

近似熵是一個理想的非線性動力學指數(shù)，基于熵理論可以刻畫出系統(tǒng)的混亂程度，同時能反映水文序列在模式上的自相似程度。需要的數(shù)據(jù)點較少，抗噪能力較好，對序列的假設比較寬松，適用于多種信號分析。所以，近似熵在生物、醫(yī)學、機械故障診斷、通信信號分析等領域都有廣泛應用。近似熵方法的主要思想及建模步驟如下[23–24]：

第1步，對時間序列 {u(i),i=1,2,···,N}進行相空間重構，重構維數(shù)為m，據(jù)此可構造一組維數(shù)為m的新向量X(1),X(2),···,X(N?m+1)，其中，X(i)=(u(i),u(i+1),···,u(i+m?1)),i=1,2,···,N?m+1。

1.3 數(shù)據(jù)滑動技術

數(shù)據(jù)滑動技術用于剔除水文序列的隨機性，其建模步驟如下：第1步，選擇滑動移除數(shù)據(jù)的窗口長度L；第2步，從待分析時間序列的第i（i=1,2,···,N?L+1，N為時間序列的總記錄個數(shù)）個數(shù)據(jù)開始連續(xù)移除L個數(shù)據(jù)，再將剩余N?L個數(shù)據(jù)直接連在一起得到一個新的時間序列；第3步，利用近似熵方法計算新序列的近似熵值；最后，固定滑動步長L不斷重復第2步和第3步直至序列結束，即可得到一個隨窗口移動而變化的近似熵序列。

1.4 耦合模型

本文構建一種水文變異診斷耦合模型，建模步驟如下：首先，以近似熵刻畫水文序列各點的動力學狀態(tài)，基于數(shù)據(jù)滑動技術充分利用水文序列的各點信息，得到隨著滑動窗口移動而變化的近似熵值。其次，運用BG算法依次對各點進行一分為二的迭代計算、對比。再次，根據(jù)提前設定的置信度檢測出所有的變異點，以變異點對近似熵序列進行動力學狀態(tài)的劃分，由此揭示水文序列的內(nèi)在演變規(guī)律。本文提出的耦合模型，具有客觀性、穩(wěn)定性和直觀性。對一個由N個點組成的水文序列X(t)來說，耦合模型建模步驟如下：

第1步，利用近似熵和數(shù)據(jù)滑動技術對X(t)進行處理，從而得到一個待檢測的近似熵序列X′(t)，初步得到了刻畫水文要素內(nèi)在演變機理的量化值。

第2步，對近似熵序列各點逐一進行一分為二的檢測分析，從而得到近似熵序列各點T檢驗統(tǒng)計值序列和置信概率序列，統(tǒng)計值和置信概率的計算參見式（2）和（4）。

第3步，根據(jù)第2步得到的分割置信概率序列對原序列X′(t)進行水文動力學狀態(tài)劃分；若其置信概率大于原先設定的P0，則以此為不同狀態(tài)的分割點，也即變異點；反之，則不進行分割。

第4步，畫出X′(t)的序列圖和不同動力學狀態(tài)劃分結果，初步判斷近似熵值驟升或驟降的變異點位置，再根據(jù)置信概率計算結果進一步確定變異點的位置和數(shù)目。

綜上，水文變異診斷耦合模型的算法流程和建模步驟如圖1所示。

圖1 耦合模型的算法流程和建模步驟Fig. 1 Algorithm flow and modeling steps of the coupling model

2 模擬實驗

2.1 耦合模型在線性時間序列中的性能測試

2.1.1 線性時間序列的構建

為比較耦合模型和僅使用滑動移除近似熵的檢測結果的優(yōu)劣，本文沿用的線性理想時間序列是在王啟光等[26]構造的方程（5）基礎上，修改一些參數(shù)所得的方程（6），以檢驗耦合模型的變異診斷結果。需要說明的是，對方程（5）修改的目的是降低演化曲線的周期數(shù)和均值差異，使之存在弱性的動力學結構變異[26]。

圖2（a）和（b）分別是式（5）和（6）所構造的理想時間序列IS1和IS2的演化，由式（5）和（6）可知，IS1和IS2在t=1 001處都發(fā)生了動力學結構變異，其控制方程由簡單的正弦函數(shù)變異為正弦函數(shù)和余弦函數(shù)所組成的新方程。

通過觀察圖2，能直觀感受到IS1和IS2的不同，IS1時間序列代表的水文系統(tǒng)在1～2 000樣本點內(nèi)包含更多的周期數(shù)，在t=1000變異點前后，表現(xiàn)出更復雜的變化；IS2則反之。因此IS1和IS2分別代表了強趨勢動力學結構變異和弱趨勢動力學結構變異。

圖2 線性時間序列的演化Fig. 2 Evolution curves of linear time series

2.1.2 耦合模型對強趨勢動力學結構變異（IS1）的性能測試

令r=0.15σ,L=10、20、50、100情況下，利用耦合模型對IS1序列進行變異診斷，結果如圖3所示。

由圖3可知，不同滑動步長情形下，耦合模型均得到IS1序列在t=1000左右發(fā)生動力學結構變異，與滑動移除近似熵的變異診斷結果一致。

圖3 不同的L下使用耦合模型的IS1檢測結果Fig. 3 Results of coupling model detection for IS1 under different L

耦合模型變異診斷的置信曲線如圖4所示。

圖4表明不同滑動步長下耦合模型的變異置信度均超過了95%。因此，對于強趨勢動力學結構變異，耦合模型不僅能準確給出變異點的位置，而且能給出變異點的置信水平。而滑動移除近似熵方法只能通過直觀觀察近似熵圖像得到變異點的大概位置，具有很大的主觀性。

圖4 圖3中對應4幅圖位于t =1 000附近的變異點置信曲線Fig. 4 Confidence images of variation points near t =1 000 for the corresponding four images in Fig. 3

2.1.3 耦合模型對弱趨勢動力學結構變異（IS2）的性能測試

IS1所刻畫的動力學結構變異是一種理想的強趨勢變異；當發(fā)生的動力學結構變異較弱時，如式（6）所構造的IS2（圖2（b）），僅使用滑動移除近似熵方法無法對變異情況進行主觀判斷。

1）滑動移除近似熵的診斷結果

利用滑動移除近似熵方法對IS2進行變異診斷，結果如圖5所示。由圖5可以看出，弱變異情況下，L的取值對檢測得到的ApEn曲線形狀的影響很大，說明滑動移除近似熵方法在弱趨勢變異時對L的取值很敏感。觀察圖5，得到以下判斷結果：L=10時，能觀察到在t=1 000附近有較為明顯的變異；L=20時，無法觀察到明顯的變異；L=50時，在t=1 000觀察到變異，同時在t=1 000～1 400之間也存在一個類似的變異，無法給出確定的變異點位置和數(shù)目；L=100時，觀察不到變異。綜上，當滑動窗口長度L≥20時，滑動移除近似熵方法無法確定變異點的數(shù)目和位置。

圖5 不同L設置下的滑動移除近似熵法的IS2變異診斷結果Fig. 5 Results of moving cut approximate entropy for IS2 under different L

2）耦合模型診斷結果

在給出耦合模型結果之前，先介紹變異區(qū)間的概念。變異區(qū)間由算法偵測到的變異點所決定。當對N=2 000樣本量的序列進行診斷時，如果設滑動窗口L=l，檢測到的變異點為Y，那么變異區(qū)間為[Y?l+1,Y]。 當設定了滑動窗口的數(shù)值l后，就會把該值之前（包括該值在內(nèi)）的l個數(shù)值移除；然后，對剩余序列進行分析，得到一個近似熵值；接著，會從2l開始往前取到l+1， 由此往復直至取完nl=N。因此，變異區(qū)間的含義是：某次移除了l個數(shù)值，則剩余序列得到的近似熵值發(fā)生了異常，認為變異點就出現(xiàn)在這l個數(shù)值中。

利用耦合模型對IS2進行變異診斷，結果如圖6所示。

圖6 不同L設置下的耦合模型IS2變異診斷結果Fig. 6 Result of the the coupled model for IS2 under different values of L

由圖6可以看出：L=10時，得出的變異點為t= 990和1 020，均值曲線在該點附近呈凹谷狀，t=990處均值減弱，t=1 020處均值變異增強，在形式上雖不如階梯狀直觀，但大概能找出預先設置的變異點位置；L=20時，得出變異點為t=980，t=[961,980]為變異區(qū)間，由于構造的方程為線性的，因此在右側邊界附近檢測到一個虛假變異點；L=50時，變異點為t=1 010，t=[961,1 010]為變異區(qū)間，另外檢測到了一個虛假變異點t=600；L=100時，變異點為t=1 900，t=[1 801,1 900]為變異區(qū)間。綜合比較滑動移除近似熵方法的結果和耦合模型的診斷結果如表1所示。

從表1看出：大部分情形下滑動移除近似熵方法不能檢測出實驗設置的變異點，且即使可檢測到，也只能基于觀察給出模糊的位置；而耦合模型可以檢測出變異點及變異區(qū)間；雖然耦合模型的精度在某些滑動步長下存在一定的偏差，但都在可接受范圍內(nèi)，如L=20時有一個多余變異點，是由式（5）的線性特點影響的，但并不影響應用；在L=50情況下有一個虛假變異點，這說明耦合模型對L是敏感的，需要進一步探索出參數(shù)的合理取值以提高模型性能。綜上，在弱趨勢變異下，樣本量為2 000時，耦合模型取L=20、r=0.15σ時可以得到良好的結果。

表1 耦合模型與僅使用滑動移除近似熵方法的檢測效果對比Tab. 1 Comparison of the detection effects of the coupling model and the sliding removal approximate entropy method only

通過上述實驗，關于耦合模型可得出以下結論：對強趨勢和弱趨勢變異均有較好的檢測表現(xiàn)。對動力學結構弱趨勢變異，在特定L情況下具有良好性能，但具體L應該如何取值使耦合模型的性能更好，僅憑目前的實驗結果還無從得知，需要在不同樣本量下對r和L的取值進一步研究才能提高模型的檢測精度。

2.2 弱趨勢變異情況下耦合模型對r 、L的取值評價

為進一步提高模型的精度，對模型中的參數(shù)進行敏感性分析。耦合模型精度主要受3個參數(shù)的影響，即樣本量N、滑動窗口長度L和容許偏差r。

2.2.1r、L對近似熵序列的影響

r的影響：

1）取較大值，可以接受時間序列數(shù)據(jù)點向近似熵值轉(zhuǎn)化過程中更大的偏差，近似熵曲線波動振幅增強，趨勢變化因受到各點偏差的影響而被掩蓋；

2）取較小值，可以強化精度的要求，近似熵曲線的波動振幅會減弱，趨勢變化對于耦合模型會更清晰。

圖7能進一步說明以下規(guī)律：在樣本量N都為2 000，且L=20 時，容許偏差r從0.15σ 到0.10σ變化，近似熵序列波動振幅減弱，變異點附近的變化趨勢更明顯；從0.15σ到 0.60σ是r取值變大的過程，近似熵序列波動振幅增大明顯，在變異點附近的趨勢直接被掩蓋掉了。

圖7 N=2 000，L=20時取不同的r值的近似熵序列Fig. 7 Sequence of approximate entropy with different values of r when N=2 000 and L=20

L的影響：

由圖3可知：L越大，近似熵序列波動性越小，變化趨勢越明顯；L越小，則反之。在同一樣本量N下，L和r對診斷結果的影響并不是相互獨立的，而是交織在一起的，因此，需要構建搜索模型，通過依次生成(r,L)數(shù)對，在2 000樣本量下探索出耦合模型具有良好性能時的(r,L)取值。

2.2.2 搜索模型

確立一個搜索模型對 (r,L)進行搜索，其核心方法是遍歷算法，由于要篩選出符合精度的取值，先設定兩個指標：檢測誤差Γ （即L與 樣本量N之比的百分數(shù)）和誤差度Ψ（即耦合模型得出的變異區(qū)間的上限b或 下限a距 離理想變異點Y的最小值z與 樣本量N之比的百分數(shù)）。

其中，z=min{|b?t|,|a?t|},a=Y?L+1,b=Y。

在此基礎上，定義權重誤差ξ，ξ= 0.5Γ+0.5Ψ，為檢測誤差和誤差度的加權組合，權重均取0.5，如圖8所示。

圖8 權重誤差闡釋Fig. 8 Weight error illustration

結合上述權重誤差公式與圖8，權重的選取基于以下考慮：1）在樣本量固定時，檢測誤差Γ來源于滑動窗口L大小，控制變異區(qū)間，以防止過大。若區(qū)間過大，即變異點Y被變異區(qū)間包括，并不能僅以Y被檢測入變異區(qū)間為評判方法好壞的標準。2）為了使變異點Y不被包含在變異區(qū)間，若Y處于變異區(qū)間中，則Ψ值恒為0，誤差主要來自Γ ；若Y不在變異區(qū)間中，則產(chǎn)生了兩部分誤差，其中，一部分是Γ，另一部分通過誤差度Ψ進行量化。綜上，視兩部分誤差影響處于均勢的地位，權重各取0.5。可以在不同應用場景下對權重進行調(diào)整，建議遵守權重相加為1，并且每一個權重應不低于0.3。

對于搜索模型的搜索范圍，當N取2 000時，r應該在 (0,0.50σ]范 圍內(nèi)遍歷尋優(yōu)，步長為0.05σ。因為過高的r能使近似熵曲線整體趨勢放大，進而抹除了變異點處的變化趨勢，加大了耦合模型的檢測難度。L相 對于N不能過大或過小，認定L最大的偏差為2%，即最大取值為 2%N。在 [ 1,2%N]范 圍內(nèi)遍歷能整除N的數(shù)值是可行的，因此L分別取5、10、20、25和40。

2.2.3 搜索 (r,L)使 權重誤差 ξ最小

通過依次生成 (r,L)數(shù)對，在2 000樣本量下探索出最優(yōu)的參數(shù)對(r,L)，如表2和3所示。由表2和3可以看出，在L=5、10且r在 ( 0,0.50σ]取值時，耦合模型的檢測結果的精度比其余 (r,L)組合更高，且有聚集趨勢。

表2 N=2 000，L=5下各r取值下檢測到的結果分析Tab. 2 Analysis of the detected results under each r value when N=2 000 and L=5

進一步觀察表2可以看出,當L=5時，r∈{0.05σ,0.10σ,0.30σ,0.40σ,0.45σ}有較好的檢測結果，都僅離理想變異點t=1001一個樣本點之差，且權重誤差僅有0.15%。由表3可知，當L=10時，r∈{0.35σ,0.40σ,0.45σ,0.50σ}有較好的檢測結果，檢測到變異區(qū)間中包括了理想變異點，權重誤差為0.50%，比0.15%大。結合劉群群等[27]對非線性動力學指標近似熵在氣候突變檢測中模擬應用的研究，及馮文宏等[28]使用近似熵對涇河徑流量序列進行變異診斷，說明耦合模型對長序列變異診斷具有良好的效果。

表3 N=2 000，L=10下各r取值下檢測到的結果分析Tab. 3 Analysis of detected results under each r value when N=2 000 and L=10

因此，在N=2000,L=5及r∈{0.05σ,0.10σ,0.30σ,0.40σ,0.45σ}的參數(shù)設置下，耦合模型對于弱趨勢動力學結構變異檢測具有較好的性能。

2.3 弱趨勢變異情況下耦合模型對于不同樣本量下r 、L評價試驗

為了對模型的實用價值進行更深入的研究，接下來改變樣本量，進行之前的r、L搜索實驗，結果如表4所示。

表4 對IS2進行不同樣本量水平下的(r ,L)搜索實驗所得結果分析Tab. 4 Analysis of the results for IS2 with different sample sizes ( r,L)

由表4可知：當樣本量較高（如N=4 000～8 000）時，權重誤差均較小，耦合模型檢測到的變異區(qū)間中包括設置的理想變異點，且樣本量越大，檢測越精確；在樣本量N=60～3 000時，檢測的變異區(qū)間雖然無法包含真實變異點，但與真實變異點較接近，可以接受。

3 應用實例

3.1 窟野河溫家川徑流變異檢測

選取窟野河溫家川站1956—2019年的年徑流量序列作為研究對象，開展耦合模型的應用研究。根據(jù)表4給出的耦合模型參數(shù)取值結果，當樣本量為60左右時，取r=0.50σ,L=1。耦合模型的變異診斷結果如圖9所示。由圖9可知，近似熵序列被分成了3個不同均值的時段，代表3種不同的動力學結構的穩(wěn)定演變狀態(tài)，因此窟野河年徑流量序列在1983年和1996年發(fā)生了變異，且變異發(fā)生的概率接近100%。

圖9 使用耦合模型檢測窟野河溫家川1956—2019年年徑流量分析Fig. 9 Annual runoff analysis of 1956—2019 in Wenjiachuan，Kuye River by using the coupling model

3.2 變異成因分析

相關研究表明，由于全球氣候變暖和大規(guī)模的人類活動導致了水循環(huán)在時間和空間上的重新分配，使水循環(huán)要素如降水、蒸散發(fā)、徑流等發(fā)生顯著變化[1]。查閱窟野河溫家川1956—2019年的相關人類活動和氣候變化文獻[29–32]，可以得出以下分析結果：

1）1956—1978年間，降水是影響徑流的最主要因素。20世紀70年代末，窟野河流域大規(guī)模實行水土保持措施，包括修筑堤壩、植樹造林、梯田建設等，這些人類活動對徑流量的減少起到直接作用。20世紀80年代初期，流域內(nèi)多處中大規(guī)模的煤礦資源開采開工，加速了水土流失的進程[29–31]。自1979年以來，人類活動對徑流的影響的占比越來越大，且呈增加的趨勢，除此之外，1979—1995年降雨量減少直接導致徑流量的減少，因此在人類活動和自然因素的影響下，徑流序列在1983年發(fā)生了變異，符合實際情況。

2）20世紀90年代以來，人類活動已經(jīng)成為影響徑流的主導因素，煤炭的開采量于1996年開始出現(xiàn)大幅度的增加[32]。這顯然對徑流有著新一輪的影響，因此，徑流序列在1996年發(fā)生了變異，也符合實際情況。

4 結 論

本文建立了一種基于BG算法和近似熵的水文變異診斷耦合模型。主要建模步驟如下：首先，以近似熵刻畫水文序列各點的動力學狀態(tài)，基于數(shù)據(jù)滑動技術充分得到隨著滑動窗口移動而變化的近似熵序列；其次，運用BG算法依次對各點進行一分為二的迭代計算、對比，得到每個點的檢驗統(tǒng)計值和置信概率；最后，根據(jù)置信概率結果判斷是否需要對近似熵序列進行分割，因此，確定變異點的數(shù)目和位置。通過大量模擬實驗，證明耦合模型具有如下潛力：

1）對強趨勢和弱趨勢情況均有較好的檢測表現(xiàn)。對動力學結構弱趨勢變異，在特定L情況下具有良好的性能。

2）在N=2000,L=5 及r∈{0.05σ,0.10σ,0.30σ,0.40σ,0.45σ}的參數(shù)設置下，耦合模型對于弱趨勢動力學結構變異檢測具有較好的性能，檢測誤差較小。

3）當在樣本量較高（如N=4 000～8 000）時，權重誤差均較小，耦合模型檢測到的變異區(qū)間中包括設置的理想變異點，且樣本量越大，檢測越精確；在樣本量N=60～3 000時，檢測的變異區(qū)間雖然無法包含真實變異點，但與真實變異點較接近，可以接受。

利用耦合模型對窟野河溫家川站1956—2019年的年徑流量序列進行變異診斷研究，結果表明窟野河年徑流序列在1983年和1996年發(fā)生了變異，與實際情況基本吻合。

本文提出的耦合模型對于大樣本具有良好的性能，但是樣本量較?。ㄐ∮诘扔?0）時，模型的效果不是很理想，這一情況在滑動移除近似熵方法中更為嚴重[33]。而中國很多流域缺乏足夠的水文觀測數(shù)據(jù)，缺資料地區(qū)的水文變異診斷技術仍需進一步探索研究。