(武漢大學(xué) 水資源與水電工程國家重點(diǎn)實(shí)驗(yàn)室， 湖北 武漢 430072)

1 研究背景

對河流物質(zhì)通量的準(zhǔn)確估算，可以掌握河流各種物質(zhì)輸移的時(shí)空分布，為河流水環(huán)境治理工程提供決策依據(jù)。

河流物質(zhì)瞬時(shí)通量是瞬時(shí)流量和瞬時(shí)濃度的乘積，而時(shí)段通量則是由兩者乘積關(guān)于時(shí)間的積分[1-3]。精確計(jì)算河流物質(zhì)時(shí)段通量，理論上需要對該河段控制斷面的流量和物質(zhì)濃度做長時(shí)間連續(xù)觀測。而在實(shí)際監(jiān)測工作中，往往無法做到這一點(diǎn)。我國大中型河流水質(zhì)測站的觀測周期一般為15 d或1個(gè)月[4-5]。以這些觀測資料為基礎(chǔ)的不同年通量估算方法的準(zhǔn)確性如何一直受到廣泛關(guān)注[6-9]。

B.W.Webb[10]采用英國河流多年流量和水質(zhì)觀測數(shù)據(jù)分析了6種時(shí)段通量估算方法的表現(xiàn)和適用性。A.Coynel[11]分析了采樣頻率和采樣方法對各時(shí)段通量估算方法精度的影響。P.J.Johnes[12]從估算方法、采樣頻率、流域人口密度等方面分析了磷通量估算的不確定性。國內(nèi)學(xué)者中，富國[13-14]分析比較了5種通量估算方法的差異，并結(jié)合長江和黃河的實(shí)測數(shù)據(jù)，對5種方法做了誤差分析。郝晨林[15]對5種方法應(yīng)用于不同污染物時(shí)段通量計(jì)算時(shí)的誤差大小進(jìn)行了分析。

本文試圖以觀測較密的懸移質(zhì)泥沙數(shù)據(jù)為基礎(chǔ)，通過用不同時(shí)間間隔選取數(shù)據(jù)，對6種通量估算方法得到的結(jié)果進(jìn)行誤差分析，為水質(zhì)通量的估算提供參考。

2 數(shù)據(jù)獲取及方法

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

選取1980年長江干流32個(gè)水文測站和烏江16個(gè)水文測站的逐日流量和逐日含沙量數(shù)據(jù)，作為本次研究的基本數(shù)據(jù)[16]。

2.2 基準(zhǔn)通量的計(jì)算

根據(jù)實(shí)測每日流量和含沙量數(shù)據(jù)，基準(zhǔn)年通量計(jì)算式如式(1)所示[17]。

(1)

式中，F(xiàn)ref表示基準(zhǔn)年通量,t/a；Qi為實(shí)測日平均流量,m3/s；Ci為實(shí)測日平均含沙量,kg/m3；i為時(shí)間序列號。由于年通量的準(zhǔn)確值較難獲取，因此采用準(zhǔn)確度較高的基準(zhǔn)年通量Fref作為參照值，以評估不同方法估算結(jié)果的精度。

2.3 通量估算方法

在物質(zhì)濃度的實(shí)際測量過程中，很難做到連續(xù)監(jiān)測[18-19]。一般的觀測周期為15 d或1個(gè)月。表1給出了時(shí)段通量估算的6種常用方法。

表1 時(shí)段通量估算方法[1]Tab.1 Load estimation Methods

注：表中K表示時(shí)間轉(zhuǎn)換系數(shù)(取決于估算時(shí)間段的長短)。

2.4 評價(jià)指標(biāo)

對通量估算方法的評價(jià)分為兩方面：無偏性和穩(wěn)定性。

無偏性通過式(2)評估：

(2)

式中，α表示無偏性評價(jià)指標(biāo)，F(xiàn)mean表示通量計(jì)算均值。α表示方法本身的系統(tǒng)誤差，α越接近0，則系統(tǒng)誤差越小；α>0，則表明估算均值大于真實(shí)值；反之，則表明估算均值小于真實(shí)值。

穩(wěn)定性通過式(3)評估：

(3)

式中，β表示穩(wěn)定性評價(jià)指標(biāo)，δ表示標(biāo)準(zhǔn)差。β反映隨機(jī)計(jì)算結(jié)果相對均值的離散程度，表征算法的穩(wěn)定性。β越小，則算法越穩(wěn)定；反之，則算法穩(wěn)定性越差。

2.5 隨機(jī)模擬次數(shù)確定

將6種通量估算方法分別應(yīng)用到長江流域及烏江流域48個(gè)測站，采用7，15，30 d共3種采樣周期，分析各方法的適用性。應(yīng)用Monte Carlo方法[20]對不同采樣方案做隨機(jī)模擬，得到計(jì)算均值、標(biāo)準(zhǔn)差等統(tǒng)計(jì)學(xué)數(shù)據(jù)。為確定最小隨機(jī)計(jì)算次數(shù)，針對不同采樣方案做大量模擬計(jì)算。計(jì)算結(jié)果表明，不同采樣方案下(7，15，30 d)，當(dāng)隨機(jī)模擬次數(shù)達(dá)到600次以后，均值和標(biāo)準(zhǔn)差均趨于穩(wěn)定。因此，本次研究針對每個(gè)計(jì)算方案的隨機(jī)模擬次數(shù)均限定為600次。

3 計(jì)算結(jié)果分析

3.1 采樣頻率對通量估算結(jié)果的影響

本次研究涵蓋長江干流沿程36個(gè)測站以及烏江干流沿程12個(gè)測站。根據(jù)測站的控制流域面積以及空間分布，選取石鼓(21.42萬km2)、宜昌(100.55萬km2)、大通(170.54萬km2)、武隆(8.30萬km2)4個(gè)典型測站，分析采樣頻率(周期)對通量估算結(jié)果的影響。

圖1給出了4個(gè)典型測站在應(yīng)用不同通量估算方法時(shí)，α值與采樣周期之間的關(guān)系。由圖1可知，在宜昌站和大通站，采樣周期增大對6種估算方法的值幾乎沒有影響。在石鼓站，當(dāng)采樣周期超過15 d，C、D、E方法計(jì)算所得的α值略有微波動，而F方法計(jì)算的α值顯著增大。在武隆站，當(dāng)采樣周期超過5 d，C、D、E方法α隨著采樣周期增大而產(chǎn)生波動，且整體呈減小趨勢，減小幅度不大；而F方法所得的α值會隨著采樣周期繼續(xù)增大，急劇增大。

由圖1可知，不考慮采樣周期變化對值的影響，不同估算方法自身的值存在差異。在采樣周期較小的情況下，A、B方法α普遍小于0，而C、D、E、F方法則十分接近0。這表明，將A、B方法應(yīng)用到泥沙通量估算時(shí)，算法自身存在系統(tǒng)誤差，估算值較真實(shí)值偏小。在不同的測站，估算值較真實(shí)值偏小的百分比有較大差異。例如，大通站的僅為-0.25，而在武隆站值接近-0.50，也即估算值僅為真實(shí)值的一半。

圖2給出了4個(gè)典型測站在應(yīng)用不同通量估算方法時(shí)，β值與采樣周期之間的關(guān)系。由圖2知，A、B、C、D、E方法的β與采樣周期之間總體呈正相關(guān)關(guān)系，且接近線性相關(guān)。不同測站β的增長速率存在差異，例如采樣周期為30 d時(shí)，大通站β取值接近0.12，而宜昌站β取值則接近0.30。

當(dāng)測站控制流域面積較大時(shí)(如大通站和宜昌站)，F(xiàn)方法計(jì)算的β值與采樣周期接近線性正相關(guān)關(guān)系。而當(dāng)測站控制流域面積較小(如石鼓站和武隆站)時(shí)，F(xiàn)方法計(jì)算的B與采樣周期之間的關(guān)系較復(fù)雜，可分為兩段。當(dāng)采樣周期小于某一臨界值時(shí)，β與采樣周期接近線性正相關(guān)關(guān)系；而當(dāng)采樣周期大于這一臨界值時(shí)，β急劇增大，且β與采樣周期得相關(guān)性變?nèi)酰現(xiàn)方法基本失效。不同測站的采樣周期臨界值差異較大，如石鼓站的臨界采樣周期為21 d，而武隆站則為9 d。

圖1 α與采樣周期關(guān)系Fig.1 Relationship between α and the sampling period

圖2 β與采樣周期關(guān)系Fig.2 Relationship between β and the sampling period

在同一測站，不同估算方法值的大小關(guān)系不甚明確。以宜昌站為例，當(dāng)采樣周期為10 d時(shí)，βC>βE>βD>βA>βB>βF；而當(dāng)采樣周期為30 d時(shí)，βF>βC>βA>βE>βD>βB。βF對采樣周期的變化較為敏感，即當(dāng)采樣周期較小時(shí)，βF最小，而當(dāng)采樣周期增大到一定程度以后，βF反而最大。

3.2 流域面積對通量估算結(jié)果的影響

由前文分析知，幾種算法的表現(xiàn)與測站控制流域面積相關(guān)。當(dāng)測站控制流域面積較小時(shí)，算法的系統(tǒng)誤差會出現(xiàn)不同程度的增大，算法穩(wěn)定性也會變差。為了進(jìn)一步分析，選取長江干流沿程以及烏江沿程主要測站作為研究對象，研究測站控制流域面積對通量估算方法的具體影響。選取測站的控制流域面積最小不到0.1萬km2(如烏江大溪河鳴玉站)，最大則達(dá)到了100萬km2(如大通站)以上，基本可反映天然河流通量估算的所有情形。

圖3表示6種不同通量估算方法的α與流域面積之間的關(guān)系。由圖3知，A、B、C、D、E算法所得的α與測站控制流域面積呈正相關(guān)關(guān)系，而F算法算得的α值則呈現(xiàn)相反的規(guī)律。總體而言，計(jì)算斷面控制流域面積越大，算法系統(tǒng)誤差越小。而由α與流域面積的相關(guān)性知，當(dāng)流域面積較小時(shí)，A、B、C、D、E方法的通量估算值較基準(zhǔn)值偏小，且流域面積越小，偏小程度越大；而F方法的通量估算值較基準(zhǔn)值偏大，并且當(dāng)流域面積極小時(shí)，估算值遠(yuǎn)遠(yuǎn)大于基準(zhǔn)值，算法基本失效。

由圖3(f)可知，當(dāng)計(jì)算斷面的控制流域面積小于某一臨界值(約在1萬～10萬km2之間)時(shí)，α急劇變大，甚至超出圖中縱坐標(biāo)范圍，此時(shí)估算均值是基準(zhǔn)值的數(shù)倍，計(jì)算結(jié)果嚴(yán)重失真。流域面積臨界值與數(shù)據(jù)序列的采樣周期也有關(guān)系：采樣周期越大，這一臨界值越大。以上分析表明，F(xiàn)方法并不適用于較小流域的通量估算，尤其當(dāng)數(shù)據(jù)序列的時(shí)間跨度較大時(shí)，其適用范圍更窄。

圖4 β與流域面積關(guān)系Fig.4 Relationship between β and the drainage area

由圖4可知，β與流域面積之間呈負(fù)相關(guān)關(guān)系，即測站控制流域面積越大，算法穩(wěn)定性越好。

由圖4(a)～(e)可知，A、B、C、D、E方法規(guī)律比較接近，β與流域面積之間呈近似線性關(guān)系，且采樣周期越大斜率越小。當(dāng)流域面積較小時(shí)，β可能大于1。此時(shí)，單次估算值可能是均值的2倍以上，即使算法本身沒有系統(tǒng)誤差，估算結(jié)果也有可能偏離實(shí)際值較遠(yuǎn)。由圖4還可知，采樣周期較小時(shí)，可在一定程度上減小值。因此，在估算控制流域面積較小的斷面通量時(shí)，應(yīng)當(dāng)具備采樣頻率較高的數(shù)據(jù)序列，否則，計(jì)算結(jié)果會出現(xiàn)較大波動，可能與實(shí)測值相差數(shù)倍。

由圖4(f)知，當(dāng)流域面積較大時(shí)，F(xiàn)方法所得規(guī)律與A、B、C、D、E方法類似。當(dāng)流域面積小于某一臨界值(約在1萬～10萬km2之間)時(shí)，β遠(yuǎn)大于1，結(jié)果嚴(yán)重失真。此時(shí)，通過縮短采樣周期并不能顯著改善計(jì)算結(jié)果，F(xiàn)方法完全不適用。

綜合以上分析，6種算法的估算結(jié)果與測站控制面積之間均存在顯著相關(guān)關(guān)系。測站控制流域面積越大，算法表現(xiàn)越好。而當(dāng)測站控制流域面積小于臨界值時(shí)，算法的無偏性和穩(wěn)定性都會變差，計(jì)算結(jié)果出現(xiàn)不同程度的失真。

4 結(jié) 論

本文從無偏性和穩(wěn)定性兩方面討論了6種通量估算方法的表現(xiàn)及適用范圍，得出以下結(jié)論。

(1) A、B兩種方法存在顯著的系統(tǒng)誤差，估算值較真實(shí)值偏小。這一誤差并不會隨著采樣周期縮短而消失，且當(dāng)采樣周期變長時(shí)，系統(tǒng)誤差會增大。當(dāng)采樣頻率足夠高，且測站流域面積較大時(shí)，C、D、E、F方法的系統(tǒng)誤差幾乎為0。

(2) 6種通量估算方法的表現(xiàn)與采樣周期(頻率)、測站控制流域面積存在顯著相關(guān)關(guān)系。采樣周期變長、控制流域面積較小都會導(dǎo)致算法的無偏性和穩(wěn)定性變差。

(3) 當(dāng)采樣周期小于某一臨界值或當(dāng)測站控制流域面積大于某一臨界值時(shí)，F(xiàn)算法的表現(xiàn)急劇變差，通量估算值與真實(shí)值可能相差數(shù)十倍以上，算法基本失效。在本研究中，不同測站采樣周期的臨界值差別較大，而流域面積臨界值約在1萬～10萬km2之間。

以上結(jié)論均是基于長江干流和烏江沿程主要測站1980年的逐日流量及逐日含沙量數(shù)據(jù)得出，對評估水質(zhì)方面有關(guān)通量計(jì)算方法的精度有一定的參考意義。幾種方法應(yīng)用于其他污染物(如氮、磷等)具體的表現(xiàn)及適用范圍如何，有待進(jìn)一步研究。