基于受體模型和控制單元分區的流域污染源解析
——以永定河張家口段為例

趙含嫣，趙銳，孫源媛，鄭明霞，蘇婧*，傅雪梅，丁鴻羽

1.環境基準與風險評估國家重點實驗室，中國環境科學研究院 2.西南交通大學

流域水污染是多種污染源作用的結果，隨著城市化進程的不斷發展，人口急劇增加，工業廢水、農業灌溉以及畜禽養殖排放都可能造成流域水污染[1-2]。近年來，我國加大了對流域水環境的治理與管理，隨著許多流域環境監測系統的建立與完善，獲取了較為系統的流域水質監測數據。受地理環境、氣候因素、土地利用方式及人類活動的影響，流域水質在空間上往往呈現出異質性[3]。因此，迫切需要從監測數據中挖掘有用的信息，識別流域水質空間分布特征，判斷水體污染主要來源，幫助決策者建立有效合理的水環境管理方案[4-5]。

受體模型是通過對樣品中檢測到的具體污染物數據，建立污染物與來源因果對應關系的源解析方法。常用的受體模型有化學質量平衡(CMB)[6]、正定矩陣因子分析(PMF)[7]、絕對主成分-多元線性回歸(APCS-MLR)[8]及Unmix模型[9]等。其中APCS-MLR模型能夠從少量數據中挖掘污染源信息，識別污染源類型并計算其貢獻率，廣泛應用于河流[10-11]、湖泊[12-13]、城市河網[14]和地下水[15-16]的污染源解析中。然而對于某些集水面積較大的流域，各子流域的污染源構成往往顯示出空間差異性，受體模型難以對不同子流域的污染源進行辨析，目前鮮有針對不同水質分區的污染源解析的研究。流域水質目標管理是從流域到控制單元再到污染源的多層次體系，一些學者開展了松花江流域[17]、黃河甘肅流域[18]、大遼河[19]和贛江流域[20]的控制單元劃分研究。若能將污染源解析與控制單元分區相結合，兼顧流域和陸域環境的特點，有效識別重點管控行業與優控單元，將為流域污染治理提供更具體、更具操作性的建議。

永定河流域是海河流域最長的支流，流經河北省張家口市的6區和6縣，流域面積占張家口市總面積的48%[21]。張家口市作為2022年冬季奧林匹克運動會場地之一，其生態環境、水質保護受到全社會的關注。以永定河流域張家口段為研究對象，以控制單元作為流域水污染防治的基礎單元，采用聚類分析(CA)方法開展流域水質和陸域污染源空間特征分析并進行分區；利用因子分析(FA)和APCS-MLR模型提取各區域水質中主要污染因子并定量解析其貢獻率，結合地理信息系統(GIS)與污染源分布，確定重點管控行業與優控單元，以期為永定河流域水質改善和水環境分級分區管理提供科學建議。

1 研究區與研究方法

1.1 研究區

永定河流域張家口段(113°49′E～115°58′E，39°33′N～41°32′N)全長573 km，流域面積為1.8萬km2。流域年降水量約370 mm，境內主要支流有桑干河、洋河及清水河，其中桑干河和洋河是河北省張家口市及山西省大同市、朔州市和忻州市等流域上游地區工農業的主要水源，供給人口超過1 000萬[22]。張家口市地處北京、河北、山西、內蒙古4省(區、市)交界處，面向沿海，背靠內陸，是溝通中原與北疆，連接中西部資源產區與東部經濟帶的重要紐帶，也是河北省礦產資源較豐富市之一，主要產業有金屬礦、非金屬礦、電力、熱力生產和農副食品加工業。

張家口監測站在永定河流域張家口段共設有9個監測斷面(圖1)：洋河自西向東有左衛、響水鋪、雞鳴驛、八號橋4個監測斷面，其中八號橋位于洋河、桑干河交匯處；桑干河自西向東有小渡口、石匣里、溫泉屯3個監測斷面；清水河自北向南有北泵房、老鴉莊2個監測斷面。

1.2 控制單元劃分

通過GIS軟件進行研究區控制單元劃分。基于該區域的數字高程模型(DEM)利用水文分析模塊提取子流域邊界，并根據河流水系的實際情況進行修正；結合當地的行政區劃，以污染控制單元劃分原則為指導，在地圖上將水系分布、子流域邊界、控制斷面分布、行政邊界等指標的空間數據進行疊加，獲得控制單元，將研究區劃分為16個控制單元(圖1)。

1.3 水質與污染源分區

1.3.1水質數據來源

采用張家口市環境監測站提供的9個監測斷面2014年1月—2017年9月的水質數據，除部分斷面由于冬季低溫河水封凍、夏季高溫河水斷流而無法采樣外，其余斷面采樣頻次均為每月1次，最后得到左衛23組、響水鋪及雞鳴驛各36組、八號橋37組、石匣里34組、溫泉屯35組、小渡口33組、老鴉莊29組、北泵房27組，共290組監測數據。選取GB 3838—2002《地表水環境質量標準》中14個指標，包括溫度、pH、電導率(EC)、溶解氧(DO)、高錳酸鹽指數(CODMn)、化學需氧量(CODCr)、五日生化需氧量(BOD5)、總磷(TP)、氨氮(NH3-N)、鋅(Zn)、氟化物(F-)、硫化物、砷(As)和糞大腸菌群數，按照《水和廢水監測分析方法》[23]進行測定。為消除量綱影響，在進行聚類分析前，對數據進行Z-分數預處理，以實現正態分布。

1.3.2面源污染物流失量計算

污染源數據來自張家口市生態環境局提供的2017年張家口市環境統計數據和《張家口經濟年鑒》[21]，主要包括各行政區面積，工業企業廢水排放量，CODCr、TN、NH3-N排放量，農業種植面積，畜禽養殖量，農村常住人口數等。基于收集的數據，進行面源污染中農業種植、畜禽養殖和農村生活污水TN流失量計算。計算公式如下：

I1=A×F×E1

(1)

I2=Q2×P2×E2

(2)

I3=Q3×P3×365×E3×C3×10-9

(3)

式中：I1為農業種植的TN流失量，t；A為農業種植占用的耕地面積，hm2；F為農業種植單位耕地面積的肥料施用量，取207.15 kg/hm2(以TN計，全文同)；E1為農業種植耕地的TN地表徑流流失系數，取0.005 41[24]；I2為畜禽養殖的TN排放量，t；Q2為畜禽的飼養數量，換算為豬當量；P2為畜禽的飼養周期，取190 d；E2為畜禽的TN排污系數，取5.34 g/(頭·d)[25]；I3為生活污水排放量，t；Q3為農村常住人口數，人；P3為農村人均用水量，根據流域農村人口用水情況取70 L/(人·d)；E3為排污系數，取0.4；C3為生活污水中TN平均排放濃度，取33.05 mg/L[26-28]。

1.3.3水質與污染源空間特征分析

采用聚類分析方法開展流域水質和陸域污染源空間特征分析，聚類分析是根據對象距離遠近或相似性大小進行分類的多元統計方法[29]。基于監測斷面水質和污染源數據，先將每一個監測斷面或污染源看作一類，然后將相近程度最高的兩類合并組成一個新類，使同一類別中的樣品或變量之間的同質性盡可能高，而不同類別之間的差異性也盡可能大。根據聚類分析結果，按照控制單元對研究區進行分區，對水質和污染源空間分異特征進行探討。

1.4 受體模型與污染源解析

將因子分析用于定性確定監測斷面的污染源，通過最大方差旋轉后獲得最大因子，選取特征值大于1的因子。在獲得污染因子之后，通過多元線性回歸(MLR)分析計算污染因子的定量貢獻。將絕對因子得分作為自變量，將污染物的測量濃度作因變量，利用多元線性回歸分析得到每個污染源的貢獻率。具體計算公式如下:

(4)

A0(k,n)=S(k,m)×Z0(m,n)

(5)

Ax(k,n)=S(k,m)×Z(m,n)

(6)

APCS=Ax-A0

(7)

以APCS為自變量，各指標的標準化值作因變量，利用多元線性回歸分析得到每個因子的貢獻率。計算公式如下：

(8)

(9)

式中：C為每個樣本各污染物標準化濃度之和；b為回歸方程的常數項；tj為因子j(j=1,2,…,k)的污染貢獻率，%；rj為因子j對濃度之和的回歸系數。

數據預處理和多元統計過程采用Excel 2016和SPSS 24軟件，空間分布展示采用ArcGIS 10.2軟件。

2 結果與討論

2.1 河流水質污染特征

永定河張家口段各監測斷面2014—2017年水質指標統計見表1。由表1可知，洋河、桑干河、清水河均受到不同程度的污染，CODMn、CODCr、BOD5和TP、NH3-N、F-濃度均出現超過GB 3838—2002中Ⅲ類水質標準的情況。其中，洋河F-、TP污染較為顯著，F-濃度平均值超過Ⅲ類水質標準，TP濃度最大值超過Ⅳ類水質標準；清水河BOD5超標較為顯著，最大值超過Ⅳ類水質標準。總體來看，洋河、清水河的水質污染程度較桑干河更嚴重。邵志江等[30]曾對永定河上游區域2013—2016年的監測數據進行統計，發現主要污染指標為TP和F-，且洋河污染較嚴重，清水河污染較輕，而桑干河水質清潔。

表1 研究區2014—2017年水質指標統計Table 1 Statistical description of water quality indicators in the study area during 2014-2017

2.2 河流水質空間特征

采用平方歐式距離進行測量，選擇離差平方和(Ward)法將9個監測斷面采樣數據的平均值進行空間聚類分析，結果如圖2所示。由圖2可知，9個監測斷面在空間上可分為兩大類，其中響水鋪、雞鳴驛、八號橋、老鴉莊聚為一類(A組)，左衛、石匣里、溫泉屯、小渡口、北泵房聚為一類(B組)。A、B 2組監測斷面水質指標標準化值的箱線圖如圖3所示。由圖3可知，A、B 2組監測斷面水質污染程度差異較大，A組監測斷面電導率(EC)、營養鹽類指標(CODMn、BOD5、NH3-N、CODCr、TP)、重金屬類指標(Zn、As)、其他污染物(F-、硫化物)和微生物指標(糞大腸桿菌群數)均明顯高于B組監測斷面，且A組監測斷面CODMn、NH3-N、CODCr、TP、F-等指標較多存在超過Ⅲ類水質的情況。因此，A組監測斷面為污染較嚴重的區域，重點控制指標為NH3-N、TP和F-。

圖2 9個監測斷面基于Ward法的聚類分析譜系Fig.2 Pedigree of cluster analysis of 9 monitoring sites based on Ward method

圖3 研究區水質指標空間差異Fig.3 Spatial variation of water quality indicators in the study area

根據斷面分類結果，結合控制單元的劃分，取控制單元內部或下游斷面為控制單元所屬斷面，將研究區劃分為A、B 2個區，結果如圖4所示。由圖4可知，A區包含4、6、8～13號控制單元，B區包含1～3、5、14～16號控制單元。從位置上看，A區的老鴉莊位于張家口市橋東區、橋西區和經濟開發區交界處，響水鋪位于宣化區，雞鳴驛和八號橋位于懷來縣，是城鎮化水平較高和農業種植較為發達的地區；B區的左衛位于懷安縣，石匣里、小渡口位于陽原縣，溫泉屯位于涿鹿縣，北泵房位于崇禮區，屬縣域中心城鎮。從人口和面積來看，A區覆蓋的區縣人口數占張家口市總人口的43.1%，面積占張家口市總面積的17.4%；B區覆蓋的區縣人口數占張家口市總人口的31.4%，面積占張家口市總面積的48.6%[21]。根據張家口市2017年環境統計數據，A區包含了研究區80%的工業源(圖4)，因此，工業源與城市及農村生活污染可能是A區監測斷面水質較差的主要原因。徐華山等[31]在漳衛南運河、李義祿等[14]在蘇州古城區水質空間分布特征的研究中也發現，城中村、人口密度、工業廢水和生活污水排放是導致水質空間分異的主要原因。因此，建議在老鴉莊至八號橋斷面之間增加沿程水質監測點。

圖4 研究區水質污染分區及工業企業分布Fig.4 Water pollution classification zonation results and industry distribution in the study area

2.3 水質污染因子識別

將A、B兩組監測斷面的原始數據進行因子分析，分析各組水質變化的關鍵因子，識別主要污染源，結果如表2、表3所示。用Kaiser-Meyer-Olkin(KMO)和Bartlett球形檢驗方法對監測數據進行相關矩陣檢驗，A組、B組的KMO分別為0.563、0.524，Bartlett顯著性分別為0.000、0.000，滿足P<0.05置信區間，檢驗結果表明因子分析是有效的。

表2 A組因子分析旋轉成分矩陣Table 2 Group A factor analysis rotation component matrix

表3 B組因子分析旋轉成分矩陣Table 3 Group B factor analysis rotation component matrix

A組提取出6個因子，累計方差為67.11%，分別提取出高荷載變量Zn、DO(F1)，F-、糞大腸菌群數(F2)，BOD5、EC(F3)，NH3-N、CODCr(F4)，TP、pH(F5)和CODMn、As(F6)。B組提取出6個因子，累計方差為62.250%，分別提取出高荷載變量DO、As、糞大腸菌群(F1)，CODMn(F2)，溫度、EC(F3)，NH3-N(F4)，Zn(F5)和CODCr(F6)。

圖5 研究區土地利用現狀Fig.5 Land use in the study area

B組所對應的B區覆蓋范圍主要包括尚義縣、崇禮區、懷安縣、萬全區、陽原縣、蔚縣和涿鹿縣。由表3可知，B組F1的主要荷載為DO、As和糞大腸菌群數，DO濃度增加可能是夏季水生植物的光合作用增強導致，張家口市主要農作物為玉米和土豆，夏季6—7月正值耕種施肥期，因此F1表征農業種植施用的除草劑、糞肥所造成的污染[41-42]。F2的主要荷載為CODMn，B區除陽原縣外都是以旅游業為支柱型產業[32]，且B區城鎮用地較少而農村居民點分布較多，由此將F2表征為由農村生活及旅游產生的面源污染。F3的主要荷載為溫度和EC，表征季節變化導致的水體物理化學特性的改變[43]，為自然因素。F4的主要荷載為NH3-N，B區的尚義縣、懷安縣、萬全區及崇禮區的城市職能均以畜產品加工制造為主[32]，而區域內畜禽養殖場多未設配套糞便污水處理設施，污水直接排入農田灌溉，因此將F4表征為由畜禽糞便造成的面源污染。F5的主要荷載為Zn，B區內蔚縣是張家口市重要的煤炭基地，涿鹿和陽原縣以礦山開發、建材和輕工業為主[32]，結合文獻梳理結果，將F5表征為經雨水沖刷進入地表水體的礦區地表徑流[32,44]。F6的主要荷載為CODCr，結合圖4和文獻[32]，B區大部分區縣的工業以食品、煙酒制造等輕工業為主，各制造業得到迅速發展的同時，其排放廢水具有CODCr高、處理達標率低、水質變化幅度大等特點[45]，因此將F6表征為區域內工業點源的排放。

綜合分析A、B組的因子荷載矩陣可知，A組水質主要受到生活源、工業點源、農業種植污染的混合影響，主要污染行業為采礦、冶金和食品加工業。A區冶金類行業管控不當將導致土壤或流域內重金屬元素的富集，從而污染地表水和地下水，建議加強對重點企業排污口的監管，杜絕污水直排現象。B組水質主要受到農業種植、農村生活、畜禽養殖、旅游等面源和工業點源的影響，主要污染行業為采礦業和食品制造業。與A組不同的是，B組的采礦業污染主要源自露天礦區及礦山開發等面源污染，因此建議在加強畜禽養殖管理的同時，還應加強礦區內污水的收集與轉運設施的建設。

2.4 污染源貢獻率分析

為了更明確地了解每種污染源對水質的影響程度，根據因子分析的結果，由式(1)～式(6)計算每個因子的污染貢獻率，結果如表4所示。由表4可知，A組可決系數(R2)為0.713，在95%的置信水平下，A組中F2自然因素、F5化肥流失的顯著性分別為0.081和0.499，大于0.05，說明在該顯著性水平下不顯著，即F2、F5 2個自變量對因變量的影響不大；B組R2為0.853，顯著性整體接近零，說明所建立的回歸模型擬合較好，可以解釋原有自變量與因變量之間的關系。綜上，本次建立的回歸模型具有統計學意義。

從表4可知，A組污染因子F6所代表的農藥、除草劑等農業種植的有機污染貢獻率最大，為44%；其次是表征點源污染的F4和F1，分別23%和20%。B組污染因子F2所代表的農村生活及旅游產生的面源污染貢獻率占比最高，為30%；其次是F6所代表的食品、煙酒制造等行業的點源排放貢獻率，為19%；另外F1及F4等面源污染貢獻率也分別達到18%和17%。總體來說，A區受點源和面源的影響程度相當，這與趙建國等[46]提出的永定河懷來段化肥與農藥施加是流域面源污染主要來源的結論相一致。因此建議A區在加強排污口監管的同時，控制農業種植的施肥量和用藥量，對農村的排水渠、泄洪渠進行清污和綠化修復，尤其是對懷來縣葡萄種植基地的地表徑流進行凈化后排放，避免直接排入官廳水庫造成污染。B區受面源影響較大，清水河上游的崇禮區還存在部分生活污水直排入河現象，且自北京市成功申辦2022年冬季奧林匹克運動會之后，崇禮區旅游業發展較快，游客較多[30]，建議增加B區鄉鎮農村的污水處理設施，加強旅游景區的基礎設施建設，重視旅游垃圾、餐飲廢水的收集與處理。

表4 A、B兩組不同污染源對污染物總量的貢獻率Table 4 Contribution rate of different pollution sources of group A and B to total pollutant amount

2.5 污染源空間特征分析

為驗證水質分析的結果，探討污染源與水質空間分布的內在聯系，并確定優先管控的控制單元，收集了研究區內33家工業企業的信息，包括年廢水排放量，CODCr、TN、NH3-N年排放量等，將工業點源統計數據與面源計算數據進行聚類分析，并進行污染源空間差異性展示，結果如圖6、圖7所示。由圖6(a)可知，33家工業企業可分為2類，其中類別1有10家，類別2有23家；類別1位于A、B區各5家，類別2位于A區有20家、B區有3家(圖7)。聚類分析結果顯示，污染源分布具有一定的規律性，類別2中大部分工業企業位于A區，且沿洋河分布。經調查發現，類別2中的主要企業屬于食品制造業等輕工業，類別1中的主要企業屬于采礦業。據《2017年張家口環境質量報告書》[47]，2016年工業廢水排放量位于前三的行業分別是煤炭開采與洗選業、黑色金屬冶煉與壓延加工業及農副食品加工業，位于前三的行政區分別是蔚縣、宣化區和涿鹿縣。蔚縣是河北省重要的煤炭工業基地，煤田已探明儲量達15億t；位于宣化區的宣鋼集團是國家重點大型企業，其鋼材年產量達800萬t，是張家口市的主要鋼鐵產地；涿鹿縣以礦山開發和林果業加工為主。結合表4污染源貢獻率分析結果，將金屬冶煉與食品制造業作為A區的重點管控行業，將采礦業與食品制造業作為B區的重點管控行業，建議加強洗煤廢水循環利用、高耗水企業廢水深度處理與回用。

圖6 基于Ward法的污染源聚類分析譜系Fig.6 Cluster analysis pedigree of pollution sources based on Ward method

由圖6(b)及圖7可知，研究區不同控制單元面源污染物排放量聚類后，可分為面源類別Ⅰ區、面源類別Ⅱ區2類區域，該分區與水質分區大部分區域相吻合。面源類別Ⅰ區的負荷明顯高于面源類別Ⅱ區。張家口市的農村多傍水，根據張家口市相關資料，位于面源類別Ⅰ區的傍水村莊有143座，常住人口為152 376人；位于面源類別Ⅱ區的傍水村莊有186座，常住人口為143 963人，控制單元內的傍河村莊和人口數是造成面源污染排放差異的主要原因[48]。2類區域中，農村生活污染負荷均最高，這主要由于流域內農村生活污水未經污水廠處理，部分傍水農村甚至將污水直接排入河中或滲入地下導致的；其次是畜禽養殖污染負荷，這就要求畜禽養殖場應規范排水設施，減少使用水沖清糞，多使用墊草墊料清糞方式以減少污染物的流失。值得注意的是，2、3、5和14號控制單元處于水質較好的B區，但卻屬于面源排放較高的類別Ⅰ區，因此將其作為預防和優先防控的控制單元；對位于面源類別Ⅱ區而水質較差的10、12號控制單元，則應加強對區域內工業點源的監管。污染源解析與控制單元分區相結合的方法較單一的采用APCS-MLR受體模型更能反映水質的空間分異特征，提供更多污染源的信息，結合污染源空間特征分析，可提高源解析能力，是研究區內重點管控行業與優控單元識別行之有效的方法。

注：圖中黑色數字為控制單元編號；灰色數字為工業企業編號。圖7 研究區污染源空間差異性分析Fig.7 Analysis of spatial variation of pollution sources

3 結論

(1)結合DEM數據和水系分布，將研究區分為16個控制單元。按水質受污染程度分為污染較重的A區(洋河、清水河中下游區域)及污染較輕的B區(洋河、清水河上游及桑干河區域)。A區水質主要受工業點源(43%)和面源(44%)的混合影響，污染主要來源于冶金、采礦、農副食品加工、食品制造業和農業種植；B區水質主要受面源(76%)影響，污染來源于農村生活及旅游產生的面源污染。

(2)將污染源數據進行聚類分析，確定A區的重點管控行業為金屬冶煉和食品制造業，B區的重點管控行業為采礦業和食品制造業。將水質較好但面源排放量較高的覆蓋陽原縣的2、3號，覆蓋涿鹿縣北部和蔚縣北部的5號、覆蓋萬全區的14號控制單元作為面源污染優先控制單元。