大汶河流域泰安段水質變化及其與水量響應關系研究

侯恩光，張升第，王如巖，張春霞，楊會清，段建筑，高振勇，譚秀翠*

（1.泰安市水文局，山東 泰安271018；2.山東農業大學，山東 泰安271018）

0 引 言

【研究意義】受自然因素與人類活動的共同影響，地表河流水質發生了變化[1]，河流水質的時空變化特征研究已經成為水資源管理領域的重要工作[2]。大汶河作為山東省泰安市的母親河，其水量豐富，但年際和年內分配不均，近些年來，受沿岸工農業快速發展的影響，部分河段水質受到明顯影響，危及流域水生態系統。因此，開展大汶河流域泰安段的水質變化及其與流量的響應關系研究勢在必行。【研究進展】水質與水量是水資源的雙重屬性，二者之間有著緊密的關系。一般來說，河流水質在1年內的變化都是由不同時期水量大小不同的差異所引起的[3]。因此，為分析不同時期水量差異所引起的水質變化，常采用季節性Kendall 檢驗法分析水質的變化趨勢[4-5]，但該方法并未將水量與水質建立關系。采用數值模型研究水量水質關系，需要輸入水動力與水質的各項參數[6-7]，其直接決定模型的計算精度。為了更真實地反映流域的水質狀況，Copula 函數被應用于構建水量與水質的聯合分布，分析二者的變化關系。2003年意大利學者De Michele 等[8]將Copula 函數應用于水文水資源領域，之后該方法得到廣泛應用。張娜等[9]采用Gumbel-Hougaard Copula 函數構建了年最大日雨量與年最大7日雨量之間的聯合分布，并基于聯合分布推求設計暴雨值。杜懿等[10]研究顯示二元t-Copula 函數最能體現郁江南寧水文站洪水峰量聯合分布的特性規律。黃鋒等[11]采用Frank Copula 函數建立年最大洪水發生時間和洪峰流量兩變量的聯合分布，有效挖掘了入庫洪水信息。廖顯薇等[12]利用Copula 函數構建二維干旱變量的聯合分布。由于Copula 函數能夠靈活地構造邊緣分布為任意分布的水文變量聯合分布，張翔等[13]將Copula 函數引入到了水量水質的研究中，構建二者的聯合分布函數，并對淮河流域蚌埠閘的水量水質聯合概率分布函數進行計算。【切入點】目前對于大汶河流域水質的研究主要以各類水質指標為基礎，開展水質評價、納污能力計算[14]、水生態健康評價[15]等，而對于水量影響下的水質變化規律考慮不夠充分。【擬解決的關鍵問題】因此，以大汶河流域泰安段為研究對象，在分析水質時空變化規律及演變趨勢的基礎上，采用Copula函數構建二維水量（流量）與水質、三維水量（降水、流量）與水質的聯合分布，計算聯合概率值，明晰水質與水量的響應關系，研究成果對大汶河流域水資源管理和保護方案的合理制定、實現水資源高效利用具有重要意義。

1 材料與方法

1.1 研究區概況

大汶河是黃河下游最大的一條支流，發源于濟南市鋼城區，自東向西流經濟南、泰安2 地，大汶河流域面積8 945 km2，年平均降水量694.5 mm，流域年平均氣溫12.9 ℃[16]，流域分布概況見圖1。根據2018年《泰安市水資源公報》，大汶河區水資源總量15.48億m3，占泰安市水資源總量的82%，是地區發展的重要水源。大汶河水質主要受沿岸工礦企業、居民生活區及農業面源污染影響，2018年大汶河水系評價總河長307 km，其中，Ⅲ類水占4.6%，Ⅳ類水占67.7%，劣Ⅴ類水占27.7%。《泰安市2018年度區域水功能區水質監測報告》顯示，2018年泰安市主要入河排污口污水入河量為2.03 億m3/a，其中化學需氧量（COD）占27%，氨氮（NH3-N）占4%，這2 個水質指標為泰安市水功能區限制納污控制指標。

圖1 大汶河流域分布Fig.1 Distribution of the Dawen river basin

為研究大汶河流域泰安段水質變化規律及其與水量變化的響應關系，本文采用2010—2019年大汶河流域泰安段干流及支流分布的7 個水質站點（漸汶河、北望、樓德、大汶口、南李村、障城、戴村壩）的2 個主要污染指標化學需氧量（COD）和氨氮（NH3-N）及1 個國家級水文站（大汶口）的逐月實測資料（降水量、流量）進行分析。

7 個水質站點位于泰安市水功能區劃中的工業用水區，其目標水質為Ⅳ類水。

1.2 研究方法

1.2.1 季節性Kendall 檢驗法

月尺度的水質數據具有很強的季節性特征，因此，采用季節性Kendall 檢驗法對大汶河流域泰安段的水質進行趨勢分析。該方法可以實現對多年同一月份的水質質量濃度進行比對，能有效地避免流量周期性變化帶來的影響[17-18]，進而揭示水質隨時間的變化規律及發展趨勢。

當有n年p月的水質數據，則X表示為：

式中：x11，…，xnp為每月水質質量濃度數據，p≤12。

對各年第i月的水質質量濃度值進行比較，并計算Si，計算式為：

式中：1≤k

當n≥10，S服從正態分布，且標準差Z滿足，

對于已知的趨勢檢驗顯著性水平α，取值0.1 和0.01，當|Z|≤Zα/2，則滿足0 假設條件。

1.2.2 Copula 聯合分布

Coupla 函數是定義域為[0，1]均勻分布的多維聯合分布函數，可以利用它將多個隨機變量的邊緣分布聯結起來構造聯合分布函數[13]，

式中：C為Copula 函數；θ為Copula 參數；F1，F2，…，Fn為各隨機變量的邊際分布。

常用的Copula函數有多種類型，其中Archimedean型在水文領域最為常用，常用的函數有Clayton Copula函數、Gumbel Copula 函數、Frank Copula 函數3 種，其函數表達式見表1[19]。

表1 3 種Copula 函數表達式Table 1 Three kinds of Copula function expressions

采用下式計算聯合分布函數的經驗頻率Pei。

式中：mi為觀測值中滿足X≤xi，Y≤yi（二維聯合分布）或滿足X≤xi，Y≤yi，Z≤zi（三維聯合分布）的聯合觀測值的個數；n為總的數據對個數。

為選擇最優的Copula 函數作為聯合分布函數，采用離差平方和最小準則（OLS）、AIC信息準則來評價Copula 函數的擬合優度[20]。

OLS的計算式為：

式中：Pei、Pi分別為聯合經驗頻率和聯合理論頻率；i為樣本序號；n為樣本容量。OLS值越小，Copula函數的擬合效果越好。

AIC信息準則計算式為：

式中：m為模型參數的個數。AIC值越小，Copula 函數的擬合效果越好。

2 結果與分析

2.1 水質變化趨勢

2010—2019年，大汶河流域泰安段COD 與NH3-N 質量濃度均呈下降趨勢（圖2），分別由2010年的34.31、1.98 mg/L 降低到2019年的23.57、0.57 mg/L，水質呈變好趨勢。在空間上，由上游漸汶河站到下游戴村壩站，COD 與NH3-N 的質量濃度呈下降趨勢，但質量濃度最高的斷面出現在樓德站，其位于大汶河的支流柴汶河上，其COD 與NH3-N 的質量濃度分別為31.24、1.56 mg/L，均超過Ⅳ類水質標準，受其影響，位于其下游的干流水質站大汶口站的COD 與NH3-N 質量濃度也較高。通過比較可以看出，樓德站為水污染防治的重點治理斷面。

圖2 2010—2019年大汶河流域泰安段水質時空變化Fig.2 Temporal and spatial variation of water quality in Taian Section of Dawen River Basin from 2010 to 2019

采用季節性Kendall 檢驗方法對大汶河流域泰安段的逐月水質數據進行趨勢分析。

表2 為水質變化趨勢分析結果，7 個站點2 個指標的14 次評價結果中，10 個評價結果呈高度顯著下降趨勢，1 個呈顯著下降趨勢，表明整體水質呈改善趨勢。雖然，漸汶河站的NH3-N 質量濃度呈顯著上升趨勢，但其平均質量濃度較低，為0.62 mg/L，可達到Ⅲ類水質標準，未呈現污染狀態。有2 個評價結果為變化趨勢不明顯，其中漸汶河站COD 質量濃度滿足目標水質要求，但樓德站的NH3-N 質量濃度相對較高，在部分月份達到劣Ⅴ類水質，占比19%，是水污染防治的重點治理指標。

表2 水質變化趨勢分析Table 2 Results of water quality change trend analysis

2.2 水質水量變化關系

河流流量具有周期性變化，河流水質組分質量濃度大多受流量周期性變化的影響[21]，而呈現出不同的變化特征。因此，本文分析大汶河流域泰安段7 個水質站點汛期（6—9月）與非汛期（其余月份）的水質變化規律。由圖3 可知，大汶河流域泰安段的7 個水質站點，汛期的水質優于非汛期，汛期COD 和NH3-N 的質量濃度均達到Ⅳ類水質標準，其中南李村、障城、戴村壩3 個斷面的NH3-N 的質量濃度低于0.5 mg/L，達到II 類水質標準。在非汛期，水質質量濃度相對較高，其中樓德站的非汛期的COD 與NH3-N的質量濃度均不滿足水功能區目標水質要求。

圖3 不同時期水質變化關系Fig.3 The relationship of water quality change in different periods

水質和水量是河流不可分割的兩方面屬性，質以量為載體，量的多少直接影響其環境承載力的大小[22]。因此，本文以大汶口站逐月實測降水量（P）、流量（Q）、化學需氧量（COD）、氨氮（NH3-N）為水量水質的分析指標，進行相關分析。由非汛期到汛期，大汶口站月平均降水量由22.26 mm 增大到128.04 mm，月平均流量由8.86 m3/s 增大到52.02 m3/s，COD 與NH3-N 的質量濃度分別下降6%、37%，Ⅳ類水達標月份占比由73%提高到83%。

由大汶口站水量與水質的相關分析圖可以看出，見圖4，隨著P、Q的增大，COD 與NH3-N 質量濃度呈減小趨勢，質量濃度峰值分別為45、7.66 mg/L，出現在2012年5月與2011年2月，其對應的降水量與流量均低于月平均值。在非汛期，數據點比較集中，靠近Z軸沿縱向分布，在汛期，數據點比較分散，在空間上比較接近于XY平面。綜上可知，大汶河泰安段水量與水質之間存在較為密切的聯系。

圖4 不同時期水量與水質變化關系Fig.4 The relationship between water quantity and water quality in different periods

2.3 水量水質聯合分布

河流的水量與水質是相互聯系、相互影響的，為反映水量變化影響下的水質變化特征，有必要進行水量水質的聯合分析，采用Copula 函數構建大汶河流域泰安段的水量水質的二維及三維聯合分布，描述二者之間的相關性結構。

Copula 函數的參數估計一般采用相關性指標法、適線法、極大似然法[23]，采用相關性指標法確定二維Copula 函數參數，采用極大似然法確定三維Copula函數參數，參數取值見表3。

表3 Copula 函數參數估計Table 3 Parameter estimation of Copula functions

對常用的聯合分布函數Clayton Copula 函數、Gumbel Copula 函數、Frank Copula 函數進行擬合優度檢驗，選擇OLS 與AIC 值最小的Copula 函數作為構造水量水質聯合分布的函數。由表4 和表5 可以確定，水量與水質的二維聯合分布采用Frank Copula 函數，三維聯合分布采用Clayton Copula 函數。

表4 二維Copula 函數的擬合優度檢驗Table 4 The goodness-of-fit test for two-dimensional Copula function

表5 三維Copula 函數的擬合優度檢驗Table 6 The goodness-of-fit test for three-dimensional Copula function

大汶口站水量與水質的二維、三維經驗聯合分布概率值與Copula 函數理論聯合分布概率值相關關系，見圖5、圖6。由圖5、圖6 可知，對于選定的二維Frank Copula 函數和三維Clayton Copula 函數，理論聯合分布與經驗聯合分布的擬合效果較好，R2均在0.85 以上，說明選取的二維、三維Copula 函數是合理的，可用來分析水量水質聯合分布問題。

圖5 二維Frank Copula 函數的聯合分布概率Fig.5 Joint distribution probability of two-dimensional Frank Copula function

圖6 三維Clayton 函數的聯合分布概率Fig.6 Joint distribution probability of three-dimensional Clayton function

采用Frank Copula函數構建水量與水質二維聯合分布，繪制概率分布圖及等值線圖，如圖7、圖8 所示。流量與COD、NH3-N 的二維聯合分布規律比較類似，即當流量、水質指標質量濃度增大時，二者的聯合概率值越大；當聯合概率一定時，COD、NH3-N質量濃度隨著流量的增大而降低，并趨于穩定；當流量一定，在任意聯合概率分布下，COD 質量濃度高于NH3-N 質量濃度；等值線的疏密，表明了水質質量濃度隨流量變化速率的大小，聯合概率低于0.7 時，等值線較為密集，水質質量濃度變幅較小。根據等值線圖，可以判定任意流量情況下，COD、NH3-N 質量濃度的發生概率，其分析結果對于水污染防治、合理開展水環境規劃，實現水資源的可持續利用具有重要的實踐意義。

圖7 流量與COD 的二維聯合分布Fig.7 Two-dimensional joint distribution of flow and COD

圖8 流量與NH3-N 的二維聯合分布Fig.8 Two-dimensional joint distribution of flow and NH3-N

采用Clayton Copula 函數構建水量與水質三維聯合分布圖，結果如圖9、圖10 所示。

圖9 P、Q 與COD 的三維聯合分布概率Fig.9 Three-dimensional joint distribution probability diagram of precipitation,flow and COD

圖10 P、Q 與NH3-N 的三維聯合分布概率Fig.10 Three-dimensional joint distribution probability diagram of precipitation,flow and NH3-N

可以看出，COD 與NH3-N 的三維聯合概率分布規律并不一致。當聯合概率值小于0.6 時，隨著P、Q、水質指標質量濃度的增大，三者聯合概率值增大，其在XY 軸的投影形態比較類似，且主要沿著流量增大的方向變化。當聯合概率值大于0.6 時，COD 與NH3-N 的三維聯合概率的演變趨勢發生變化，其對應的流量基本在200 m3/s 左右，但對應的降水區間有較大差異。

3 討論

郝守寧等[24]研究顯示，水質隨河流方向有改善的趨勢，大汶河流域泰安段的水質也呈現了相同的變化規律，但不同河流水質的沿程變化特征會受到沿岸城市分布、人類活動的劇烈程度影響，陳善榮等[25]對長江干流59 個水質監測斷面進行水質變化特征研究，結果顯示上游水質好于中下游。

蔡帥[26]研究表明，受來水量影響，大連市蔡房身大橋斷面水質季節變化較為明顯，豐水期（8月）水質優于平水期（10月）和枯水期（4月）水質。這與大汶河流域泰安段7 個水質站點的汛期與非汛期水質變化規律一致。馮衛等[2]研究表明，受人為因素影響，河流沿岸點源的污水排放在豐水期、平水期和枯水期存在較大差異，會導致水質變化的時間差異。

綜上可以看出，河流水質的時空變化規律受自然因素、人為因素綜合影響，在不同的流域會呈現不同的變化特征。

本文分別采用Frank Copula、Clayton Copula 函數構建了大汶口站水量水質的二維、三維聯合分布，而張翔等[13]對淮河流域蚌埠閘的水量水質的聯合分布研究采用的是Gumbel-Hougard Copula 函數，主要是因為Copula 函數本身不具有水文物理基礎，不能從物理意義上確定各隨機變量的理論聯合分布模型，通常通過擬合優度檢驗，選擇與實測資料配合最優的Copula 函數[27]。劉章君等[27]綜述了Copula 函數在水文領域的研究進展，證實Copula 函數是一種靈活構造多變量聯合分布和處理多變量問題的有效工具，具有良好的適用性和不可替代的優越性。對于水量水質的研究，不僅用于構建二者的聯合分布模型，在水量水質聯合風險的綜合評價[28]中也體現出較大的潛力。本文對于水量水質關系的分析主要采用大汶河干流大汶口水文站的實測資料，受資料限值，缺少對其上游、下游及支流的水量水質關系分析，需要在未來的研究中做進一步分析討論。在構建水量與水質二維、三維聯合分布的基礎上，應繼續開展應用研究，開展水量水質風險分析，為水污染防治、水量水質綜合監測管理提供參考。

4 結論

1）在時間上與空間上，大汶河流域泰安段COD與NH3-N 的質量濃度均呈下降趨勢，其中樓德站為水污染防治的重點治理斷面。

2）季節性Kendall 檢驗方法結果表明，71%的評價結果為高度顯著下降趨勢。

3）大汶口站實測水量水質相關分析表明，隨著降水量與流量的增大，COD 與NH3-N 質量濃度下降。非汛期到汛期，COD 與NH3-N 的質量濃度分別下降6%、37%，Ⅳ類水達標月份占比提高10%。

4）通過擬合優度檢驗，選擇Frank Copula 函數構建大汶口站的水量與水質的二維聯合分布，COD與NH3-N 的概率分布規律基本一致。選擇Clayton Copula 函數構建水量與水質的三維聯合分布，當聯合概率值大于0.6 時，COD 與NH3-N 的概率分布規律發生變化。