不同插值方法在面源污染入河量計算中的適用性研究

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

（四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室，四川 成都 610065）

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

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

1 研究區(qū)域概況

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

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

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

2 研究方法

2.1 輸出系數(shù)法

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

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

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

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

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

本文根據(jù)國內(nèi)相似區(qū)域已有的研究結(jié)果并結(jié)合研究區(qū)的實際情況[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ù)，利用空間插值方法得到全區(qū)域污染物的入河情況[15]. 本文利用兩種插值方法進行空間分析，獲取全區(qū)域TN 的入河情況，其中采用的入河系數(shù)點來自對渠江流域構(gòu)建的SWAT 模型，根據(jù)TN 模擬結(jié)果，計算出渠江流域50 個子流域的入河系數(shù)，得到入河系數(shù)的范圍在0.005 ～0.99 之間.

2.2 普通克里金法(OK)

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

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

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

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

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

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

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

3 結(jié)果與分析

3.1 面源污染物負荷總量空間分析

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

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

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

圖2 渠江各區(qū)縣TN 負荷總量空間分布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)不同插值方法對TN 入河系數(shù)的空間分布影響

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

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

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

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

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

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

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

3.3 結(jié)果驗證

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

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

4 討論

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

5 結(jié)論

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

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

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