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

（1.泰安市水文局，山東 泰安271018；2.山東農(nóng)業(yè)大學(xué)，山東 泰安271018）

0 引 言

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

1 材料與方法

1.1 研究區(qū)概況

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

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

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

7 個水質(zhì)站點位于泰安市水功能區(qū)劃中的工業(yè)用水區(qū)，其目標(biāo)水質(zhì)為Ⅳ類水。

1.2 研究方法

1.2.1 季節(jié)性Kendall 檢驗法

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

當(dāng)有n年p月的水質(zhì)數(shù)據(jù)，則X表示為：

式中：x11，…，xnp為每月水質(zhì)質(zhì)量濃度數(shù)據(jù)，p≤12。

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

式中：1≤k

當(dāng)n≥10，S服從正態(tài)分布，且標(biāo)準(zhǔn)差Z滿足，

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

1.2.2 Copula 聯(lián)合分布

Coupla 函數(shù)是定義域為[0，1]均勻分布的多維聯(lián)合分布函數(shù)，可以利用它將多個隨機(jī)變量的邊緣分布聯(lián)結(jié)起來構(gòu)造聯(lián)合分布函數(shù)[13]，

式中：C為Copula 函數(shù)；θ為Copula 參數(shù)；F1，F(xiàn)2，…，F(xiàn)n為各隨機(jī)變量的邊際分布。

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

表1 3 種Copula 函數(shù)表達(dá)式Table 1 Three kinds of Copula function expressions

采用下式計算聯(lián)合分布函數(shù)的經(jīng)驗頻率Pei。

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

為選擇最優(yōu)的Copula 函數(shù)作為聯(lián)合分布函數(shù)，采用離差平方和最小準(zhǔn)則（OLS）、AIC信息準(zhǔn)則來評價Copula 函數(shù)的擬合優(yōu)度[20]。

OLS的計算式為：

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

AIC信息準(zhǔn)則計算式為：

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

2 結(jié)果與分析

2.1 水質(zhì)變化趨勢

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

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

采用季節(jié)性Kendall 檢驗方法對大汶河流域泰安段的逐月水質(zhì)數(shù)據(jù)進(jìn)行趨勢分析。

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

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

2.2 水質(zhì)水量變化關(guān)系

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

圖3 不同時期水質(zhì)變化關(guān)系Fig.3 The relationship of water quality change in different periods

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

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

圖4 不同時期水量與水質(zhì)變化關(guān)系Fig.4 The relationship between water quantity and water quality in different periods

2.3 水量水質(zhì)聯(lián)合分布

河流的水量與水質(zhì)是相互聯(lián)系、相互影響的，為反映水量變化影響下的水質(zhì)變化特征，有必要進(jìn)行水量水質(zhì)的聯(lián)合分析，采用Copula 函數(shù)構(gòu)建大汶河流域泰安段的水量水質(zhì)的二維及三維聯(lián)合分布，描述二者之間的相關(guān)性結(jié)構(gòu)。

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

表3 Copula 函數(shù)參數(shù)估計Table 3 Parameter estimation of Copula functions

對常用的聯(lián)合分布函數(shù)Clayton Copula 函數(shù)、Gumbel Copula 函數(shù)、Frank Copula 函數(shù)進(jìn)行擬合優(yōu)度檢驗，選擇OLS 與AIC 值最小的Copula 函數(shù)作為構(gòu)造水量水質(zhì)聯(lián)合分布的函數(shù)。由表4 和表5 可以確定，水量與水質(zhì)的二維聯(lián)合分布采用Frank Copula 函數(shù)，三維聯(lián)合分布采用Clayton Copula 函數(shù)。

表4 二維Copula 函數(shù)的擬合優(yōu)度檢驗Table 4 The goodness-of-fit test for two-dimensional Copula function

表5 三維Copula 函數(shù)的擬合優(yōu)度檢驗Table 6 The goodness-of-fit test for three-dimensional Copula function

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

圖5 二維Frank Copula 函數(shù)的聯(lián)合分布概率Fig.5 Joint distribution probability of two-dimensional Frank Copula function

圖6 三維Clayton 函數(shù)的聯(lián)合分布概率Fig.6 Joint distribution probability of three-dimensional Clayton function

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

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

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

采用Clayton Copula 函數(shù)構(gòu)建水量與水質(zhì)三維聯(lián)合分布圖，結(jié)果如圖9、圖10 所示。

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

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

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

3 討論

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

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

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

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

4 結(jié)論

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

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

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

4）通過擬合優(yōu)度檢驗，選擇Frank Copula 函數(shù)構(gòu)建大汶口站的水量與水質(zhì)的二維聯(lián)合分布，COD與NH3-N 的概率分布規(guī)律基本一致。選擇Clayton Copula 函數(shù)構(gòu)建水量與水質(zhì)的三維聯(lián)合分布，當(dāng)聯(lián)合概率值大于0.6 時，COD 與NH3-N 的概率分布規(guī)律發(fā)生變化。