章 晶，張新華，蔡育杰，金 濤，唐 瑋

（四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065）

隨著全球經(jīng)濟(jì)發(fā)展，面源污染已經(jīng)成為水污染的主要原因[1].其中污染入河量計(jì)算不僅有識(shí)別水體污染主要來(lái)源的作用，更能對(duì)污染源的特征進(jìn)行明確分析，是一項(xiàng)必要的基礎(chǔ)性工作[2]. 但面源污染所具有的隨機(jī)性、廣泛性和滯后性等特征使得在監(jiān)測(cè)與計(jì)算中的難度較大[3]，而空間插值因其能夠很好反映區(qū)域變量的空間真實(shí)變化等優(yōu)點(diǎn)，使得其在面源污染計(jì)算的應(yīng)用已經(jīng)越來(lái)越廣.目前在面源污染負(fù)荷計(jì)算應(yīng)用最廣泛的方法是1996 年Johnes[4]提出的將土地利用類型、畜禽養(yǎng)殖、農(nóng)村生活排放等不同面源類型綜合考慮的輸出系數(shù)法. 蔡明等[5](2004)通過(guò)考慮流域損失改進(jìn)了輸出系數(shù)模型，探究了面源污染入河量的情況；胡昱欣等[6](2015)利用克里金插值方法對(duì)東遼河30 個(gè)子流域出口的入河系數(shù)進(jìn)行空間分析，得到東遼河流域泥沙入河量總體較大，而面源入河量呈現(xiàn)東南小西北大的特征；胡富昶等[7](2019)采用反距離權(quán)重法對(duì)降雨影響因子進(jìn)行空間插值，分析射洪縣污染空間的分布狀況.

然而不同的插值方法因其原理不同和對(duì)樣本點(diǎn)的要求不一致，所反映的空間效果存在一定差異，這對(duì)面源污染入河量的計(jì)算存在很大影響.本文基于輸出系數(shù)法對(duì)研究區(qū)3 種污染源進(jìn)行負(fù)荷總量計(jì)算，分析該區(qū)2018 年TN 負(fù)荷總量空間分布特征；此外，針對(duì)反距離權(quán)重法受樣本點(diǎn)密度的影響較大的特征，對(duì)該方法進(jìn)行均勻增加樣本的修正手段，并進(jìn)一步對(duì)比普通克里金法和修正前后的反距離權(quán)重法在計(jì)算污染入河量的結(jié)果差異，利用渠江出口實(shí)際負(fù)荷量檢驗(yàn)插值結(jié)果的精確性，探究不同插值方法在計(jì)算面源污染物入河量的適用性.

1 研究區(qū)域概況

渠江發(fā)源于川陜交界處米倉(cāng)山系鐵船山，是長(zhǎng)江支流嘉陵江左岸最大支流. 自東北向西南經(jīng)過(guò)平昌縣、渠縣、廣安等縣境，在重慶市合川區(qū)匯入嘉陵江.由于處在亞熱帶濕潤(rùn)季風(fēng)區(qū)，受季風(fēng)影響，降水充足，多年平均降水量在1 000 ～1 600 mm 之間.徑流主要由降雨補(bǔ)給，5 ～10 月為汛期，11 月至次年4 月為枯季，徑流年內(nèi)分配極不均勻，汛期徑流量占年徑流量的86.1%，而枯期僅占13.9%[8].渠江流域作為川東北糧倉(cāng)，全流域內(nèi)耕地占比55.47%，特別是干流兩岸及支流州河流域，農(nóng)業(yè)生產(chǎn)水平較高，受面源污染影響嚴(yán)重.

研究區(qū)域?yàn)榍饔蛑饕婕暗降?9 個(gè)區(qū)縣，如圖1 所示.涉及總?cè)丝? 634.24 萬(wàn)，其中農(nóng)業(yè)人口有951.69 萬(wàn)人.本文以III 類水質(zhì)為標(biāo)準(zhǔn)，依據(jù)渠江出口在2015 ～2018 年的監(jiān)測(cè)斷面資料可知，全流域主要污染物為總氮TN，其中2016 ～2018 年渠江出口處TN 濃度全面超標(biāo)，2018 年渠江出口TN 的平均濃度為1.53 mg/L，根據(jù)地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)為V 類水質(zhì).

圖1 渠江地理位置圖Fig.1 The geographical position of Qujiang River

2 研究方法

2.1 輸出系數(shù)法

輸出系數(shù)法通過(guò)采用輸出系數(shù)估算各類農(nóng)業(yè)污染源(土地利用、畜禽養(yǎng)殖、農(nóng)村生活)排污量，計(jì)算得到各區(qū)縣污染物負(fù)荷總量.輸出系數(shù)模型為[9]:

式中:Sj為污染物j在區(qū)域的總負(fù)荷量，t/a；j為污染物類型；Eij為污染物j在區(qū)域第i種土地利用類型的輸出系數(shù)或牲畜、人口的輸出系數(shù)，t/(km2·a)；i為區(qū)域土地利用類型的種類或牲畜、人口數(shù)量；Ai為第i種土地利用類型的面積(km2)，或牲畜數(shù)量(頭)，或人口數(shù)量(人).

基于各區(qū)縣污染物負(fù)荷總量，采用入河系數(shù)計(jì)算各區(qū)縣面源污染物入水體負(fù)荷量. 入河量表達(dá)式為[10]:

式中:Lj為污染物j進(jìn)入河流水體的負(fù)荷量，t/a；λ為入河系數(shù).

2.1.1 輸出系數(shù)的確定

本文根據(jù)國(guó)內(nèi)相似區(qū)域已有的研究結(jié)果并結(jié)合研究區(qū)的實(shí)際情況[7，11-14]，確定出不同土地利用、畜禽以及農(nóng)業(yè)生活的輸出系數(shù)取值，如表1 所示.

表1 渠江流域面源污染物輸出系數(shù)表Table 1 Qujiang river basin non-point source pollutant output coefficient table

2.1.2 入河系數(shù)的確定

為了更好反映區(qū)域污染物的實(shí)際入河程度，通常以入河系數(shù)作為重要參數(shù)，利用空間插值方法得到全區(qū)域污染物的入河情況[15]. 本文利用兩種插值方法進(jìn)行空間分析，獲取全區(qū)域TN 的入河情況，其中采用的入河系數(shù)點(diǎn)來(lái)自對(duì)渠江流域構(gòu)建的SWAT 模型，根據(jù)TN 模擬結(jié)果，計(jì)算出渠江流域50 個(gè)子流域的入河系數(shù)，得到入河系數(shù)的范圍在0.005 ～0.99 之間.

2.2 普通克里金法(OK)

普通克里金法是隨著地質(zhì)統(tǒng)計(jì)學(xué)的發(fā)展而衍生出的一種無(wú)偏估計(jì)的空間插值方法.其原理是利用已知的插值樣本加權(quán)平均值估計(jì)其余待預(yù)測(cè)點(diǎn)值，該方法能使預(yù)測(cè)值等于實(shí)際值的數(shù)學(xué)期望值并保持方差最小，能保證局部誤差最小等優(yōu)點(diǎn)，但也會(huì)產(chǎn)生平滑效應(yīng)，即較小值易被夸大，較大值易被低估.其計(jì)算公式如下[16]:

式中:Z為預(yù)測(cè)值；λ表示克里金法權(quán)重系數(shù)；Z(Xi)表示實(shí)測(cè)點(diǎn)Xi處的值.

2.3 反距離權(quán)重法(IDW)

反距離權(quán)重插值是基于相近相似原理，其內(nèi)涵為離得越近的樣本之間有著相似的性質(zhì)；相反，樣本之間離得越遠(yuǎn)其相似性就越小.具體計(jì)算公式如下[16]:

式中:Z是預(yù)測(cè)值，n為總樣本數(shù)，Zi為第i個(gè)樣本點(diǎn)值，di為預(yù)測(cè)點(diǎn)與樣本點(diǎn)的分離距離，p為加權(quán)冪.

這種方法在計(jì)算中相對(duì)容易，但在入河系數(shù)的應(yīng)用上有一定局限性:1)與地質(zhì)統(tǒng)計(jì)學(xué)方法相比，IDW無(wú)法在無(wú)樣本處估計(jì)預(yù)測(cè)值方差；2)當(dāng)樣本點(diǎn)空間分布不均勻時(shí)，IDW 可能存在問(wèn)題；3)IDW 的計(jì)算精度受樣本點(diǎn)密度的影響較大[17].

針對(duì)以上問(wèn)題，本文在利用IDW 對(duì)入河系數(shù)進(jìn)行空間插值時(shí)，因樣本點(diǎn)分布不均勻，在樣點(diǎn)未覆蓋到的岳池縣、儀隴縣和大竹縣區(qū)域根據(jù)相近相似原則增加了樣本點(diǎn)，并且對(duì)污染入河量較大的渠縣區(qū)域的樣本點(diǎn)密度進(jìn)行了加密，并進(jìn)行了修正前后的對(duì)比.

3 結(jié)果與分析

3.1 面源污染物負(fù)荷總量空間分析

通過(guò)輸出系數(shù)法計(jì)算得到研究區(qū)2018 年面源污染負(fù)荷總量Sj(表2)，由表2 可知研究區(qū)產(chǎn)生的TN負(fù)荷量為67 512.46 t/a，土地利用產(chǎn)生的污染負(fù)荷量最大，其占總負(fù)荷量的84.4%；其中耕地對(duì)土地利用產(chǎn)生的污染負(fù)荷總量的貢獻(xiàn)為51 191.22 t/a，占整個(gè)土地利用產(chǎn)生的75.8%. 從行政單元來(lái)看，巴中市和達(dá)州市由于耕地面積較大，且農(nóng)村人口較多，故產(chǎn)生的面源污染負(fù)荷總量也較大.

表2 研究區(qū)面源污染負(fù)荷總量計(jì)算結(jié)果Table 2 Calculation results of total non-point source pollution load in the study area

通過(guò)Arcgis 軟件對(duì)TN 負(fù)荷量進(jìn)行空間分析，得到研究區(qū)2018 年各區(qū)縣TN 負(fù)荷量分布圖(圖2)和各區(qū)縣3 類污染源貢獻(xiàn)情況(圖3). 由圖2 看出，通江縣和宣漢縣產(chǎn)生的TN 負(fù)荷量最大，而通川區(qū)、開江縣、前鋒區(qū)和華鎣市產(chǎn)生的TN 負(fù)荷量較?。挥蓤D3 可以看出，對(duì)每個(gè)區(qū)縣來(lái)說(shuō)土地利用產(chǎn)生的污染負(fù)荷量比重均最大，農(nóng)業(yè)生活產(chǎn)生污染量次之，畜禽養(yǎng)殖污染產(chǎn)生量最小.可以推斷是因?yàn)檠芯繀^(qū)人口密度大且農(nóng)業(yè)人口占一半以上，農(nóng)民以農(nóng)業(yè)生產(chǎn)為主，加上耕地面積占比大，因而土地利用產(chǎn)生的污染負(fù)荷量偏大，另外由于農(nóng)村居民分布不集中，集水管網(wǎng)和污水處理設(shè)施的建設(shè)遠(yuǎn)遠(yuǎn)不夠，導(dǎo)致生活污染物以及生活廢水任意排放.

圖2 渠江各區(qū)縣TN 負(fù)荷總量空間分布Fig.2 Spatial distribution of total TN load in each district and county of the Qujiang River

圖3 渠江各區(qū)縣不同污染源占比情況Fig.3 The proportion of different pollution sources in each district and county of the Qujiang River

3.2 不同插值方法下TN 入河系數(shù)與入河量空間分布

1)不同插值方法對(duì)TN 入河系數(shù)的空間分布影響

由圖4 可以看出，兩種插值方法一定程度上均能反映研究區(qū)入河系數(shù)的空間分布特征，整體上大致呈東北部向西南逐漸增加后減少的趨勢(shì).研究區(qū)東北部位于渠江源頭處，植被覆蓋面積大，對(duì)污染物入河有一定程度的阻擋作用；又因上游巴河、州河在渠縣開始匯入渠江干流，污染物入河系數(shù)也逐漸增加至渠縣境內(nèi)達(dá)到最大；最后隨著地勢(shì)變緩，研究區(qū)西南部的污染物入河系數(shù)漸漸減小.

圖4 不同插值方法下的入河系數(shù)空間分布圖Fig.4 Spatial distribution of river entry coefficients under different interpolation methods

具體來(lái)看， IDW 法在修正前的插值結(jié)果整體小于其他兩個(gè)，且呈現(xiàn)出入河系數(shù)大的位點(diǎn)分布不均的狀態(tài)，這是由于插值點(diǎn)的不均勻性導(dǎo)致的. 而OK 法和修正后的IDW 法的插值結(jié)果趨勢(shì)上差距不大，但OK 法相比IDW 法存在平滑效應(yīng)，使得總體插值范圍在0.134 ～0.849 之間，與樣本點(diǎn)實(shí)際的入河系數(shù)范圍0.005 ～0.99 相比縮小了許多. 另外，IDW 法受樣本距離遠(yuǎn)近的影響，對(duì)上游南江縣、萬(wàn)源市、宣漢縣的入河系數(shù)插值結(jié)果偏小. 綜上，由于插值方法原理的不同和插值樣本的處理方式不一致，兩種插值結(jié)果在各個(gè)區(qū)縣的空間分布上存在差異.

2)不同插值方法對(duì)TN 入河量的空間分布影響

不同插值方法對(duì)TN 入河量的空間分布如圖5 所示，研究區(qū)的TN 入河量在空間上分布不均，但總體趨勢(shì)為中間高并依次向兩邊遞減. 其中通江縣、宣漢縣和平昌縣雖然位于研究區(qū)的渠江上游處，但因其較其他區(qū)縣的耕地面積大，導(dǎo)致入河量也大；渠縣的TN 入河量貢獻(xiàn)為研究區(qū)最大，這是因?yàn)榍h的入河系數(shù)大小為上游城市的兩倍左右；而廣安區(qū)雖然入河系數(shù)大，但其地域面積比渠縣小，導(dǎo)致TN 負(fù)荷總量小，入河量也減小.

圖5 渠江各區(qū)縣TN 入河量空間分布Fig.5 Qujiang River spatial distribution of TN intake in each district and county

總體來(lái)看，OK 法計(jì)算出的研究區(qū)TN 入河量在空間上更能反映各區(qū)縣污染的實(shí)際入河情況，而IDW法因在南江縣、儀隴縣境內(nèi)的插值點(diǎn)較少，使得這兩個(gè)區(qū)縣的入河量偏小. 其中IDW 法在修正前的計(jì)算結(jié)果整體上小于其余兩個(gè)計(jì)算結(jié)果，這是由于整個(gè)研究區(qū)的入河系數(shù)在空間插值后的結(jié)果偏小導(dǎo)致；而修正后的IDW 法與OK 法相比，入河量結(jié)果稍微偏大，但相差大為減小，滿足誤差要求.

3.3 結(jié)果驗(yàn)證

本文采用兩種插值方法對(duì)入河系數(shù)進(jìn)行空間分析，并對(duì)IDW 方法修正前后的結(jié)果也進(jìn)行了對(duì)比，結(jié)合輸出系數(shù)法計(jì)算出2018 年研究區(qū)TN 入河量，并選取渠江出口斷面的實(shí)際負(fù)荷量來(lái)對(duì)三種計(jì)算結(jié)果進(jìn)行驗(yàn)證.結(jié)果見表3，可以看出OK 法在對(duì)TN 入河量計(jì)算上有很好的效果，其相對(duì)誤差僅為0.73%；而修正后的IDW 法相對(duì)于修正前相對(duì)誤差降低了28.59%.由此表明IDW 法采用的相近相似原理在對(duì)TN 入河量計(jì)算上有一定的偏差，在使用時(shí)需要對(duì)插值點(diǎn)進(jìn)行修正，修正后的插值結(jié)果可以達(dá)到較好精度.

表3 結(jié)合不同插值方法的ECM 模擬結(jié)果對(duì)比Table 3 Comparison of ECM simulation results combining different interpolation methods

4 討論

從污染源角度分析，土地利用產(chǎn)生的污染負(fù)荷量達(dá)到總量的84.4%，這與研究區(qū)以農(nóng)業(yè)生產(chǎn)為主的情況相符合；從污染入河情況分析，研究區(qū)TN 平均入河系數(shù)為0.39，其中渠縣的污染流失程度最大，應(yīng)當(dāng)受到重視；從插值方法分析，不同插值方法計(jì)算出的污染入河系數(shù)空間分布略有不同，OK 法的插值結(jié)果不可避免存在平滑效應(yīng)，即夸大較小值，低估較大值，但OK 法可以保證局部誤差達(dá)到最小，其在面源污染入河計(jì)算結(jié)果上也相對(duì)最優(yōu)，而IDW 法的計(jì)算精度受插值點(diǎn)數(shù)量和密度影響較大，要得到良好的插值結(jié)果對(duì)插值點(diǎn)的空間分布情況要求會(huì)更高.

5 結(jié)論

通過(guò)輸出系數(shù)法對(duì)研究區(qū)2018 年TN 負(fù)荷總量情況進(jìn)行空間分析，并采用OK 法和修正前后的IDW法對(duì)入河系數(shù)進(jìn)行空間插值，探究了不同插值方法對(duì)TN 入河量的影響，得到結(jié)論如下:

1)在對(duì)研究區(qū)的TN 入河情況進(jìn)行空間分析時(shí)，兩種插值法均能反映研究區(qū)TN 入河情況. 采用OK法能夠更全面真實(shí)的反映各個(gè)區(qū)縣的污染入河情況，而IDW 法在使用時(shí)要注意對(duì)插值點(diǎn)的空間分布形態(tài)、均勻程度等進(jìn)行分析，對(duì)分布不均，密度不夠應(yīng)適當(dāng)修定.本文采用加密方法進(jìn)行適當(dāng)修正，修正后的IDW 法和OK 法在入河量計(jì)算上都較為精確.

2)從污染源角度來(lái)看，研究區(qū)農(nóng)業(yè)人口數(shù)量較大，耕地面積廣，在耕種過(guò)程中的各類化肥以及農(nóng)藥的使用造成耕地的TN 負(fù)荷量高，其中耕地產(chǎn)生的TN負(fù)荷量占總量的75.8%；從污染流失程度來(lái)看，渠縣作為污染流失最嚴(yán)重的地區(qū)，面臨的面源污染風(fēng)險(xiǎn)最大，是今后污染治理的首要對(duì)象.其次，要大力普及測(cè)土配方施肥，推廣病蟲害綠色防控技術(shù)，著力提升農(nóng)業(yè)科技水平，方能達(dá)到有效控制面源污染的目的.