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

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

1.環(huán)境基準(zhǔn)與風(fēng)險評估國家重點實驗室，中國環(huán)境科學(xué)研究院 2.西南交通大學(xué)

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

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

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

1 研究區(qū)與研究方法

1.1 研究區(qū)

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

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

1.2 控制單元劃分

通過GIS軟件進(jìn)行研究區(qū)控制單元劃分?；谠搮^(qū)域的數(shù)字高程模型(DEM)利用水文分析模塊提取子流域邊界，并根據(jù)河流水系的實際情況進(jìn)行修正；結(jié)合當(dāng)?shù)氐男姓^(qū)劃，以污染控制單元劃分原則為指導(dǎo)，在地圖上將水系分布、子流域邊界、控制斷面分布、行政邊界等指標(biāo)的空間數(shù)據(jù)進(jìn)行疊加，獲得控制單元，將研究區(qū)劃分為16個控制單元(圖1)。

1.3 水質(zhì)與污染源分區(qū)

1.3.1水質(zhì)數(shù)據(jù)來源

采用張家口市環(huán)境監(jiān)測站提供的9個監(jiān)測斷面2014年1月—2017年9月的水質(zhì)數(shù)據(jù)，除部分?jǐn)嗝嬗捎诙镜蜏睾铀鈨?、夏季高溫河水?dāng)嗔鞫鵁o法采樣外，其余斷面采樣頻次均為每月1次，最后得到左衛(wèi)23組、響水鋪及雞鳴驛各36組、八號橋37組、石匣里34組、溫泉屯35組、小渡口33組、老鴉莊29組、北泵房27組，共290組監(jiān)測數(shù)據(jù)。選取GB 3838—2002《地表水環(huán)境質(zhì)量標(biāo)準(zhǔn)》中14個指標(biāo)，包括溫度、pH、電導(dǎo)率(EC)、溶解氧(DO)、高錳酸鹽指數(shù)(CODMn)、化學(xué)需氧量(CODCr)、五日生化需氧量(BOD5)、總磷(TP)、氨氮(NH3-N)、鋅(Zn)、氟化物(F-)、硫化物、砷(As)和糞大腸菌群數(shù)，按照《水和廢水監(jiān)測分析方法》[23]進(jìn)行測定。為消除量綱影響，在進(jìn)行聚類分析前，對數(shù)據(jù)進(jìn)行Z-分?jǐn)?shù)預(yù)處理，以實現(xiàn)正態(tài)分布。

1.3.2面源污染物流失量計算

污染源數(shù)據(jù)來自張家口市生態(tài)環(huán)境局提供的2017年張家口市環(huán)境統(tǒng)計數(shù)據(jù)和《張家口經(jīng)濟年鑒》[21]，主要包括各行政區(qū)面積，工業(yè)企業(yè)廢水排放量，CODCr、TN、NH3-N排放量，農(nóng)業(yè)種植面積，畜禽養(yǎng)殖量，農(nóng)村常住人口數(shù)等?；谑占臄?shù)據(jù)，進(jìn)行面源污染中農(nóng)業(yè)種植、畜禽養(yǎng)殖和農(nóng)村生活污水TN流失量計算。計算公式如下：

I1=A×F×E1

(1)

I2=Q2×P2×E2

(2)

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

(3)

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

1.3.3水質(zhì)與污染源空間特征分析

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

1.4 受體模型與污染源解析

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

(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為自變量，各指標(biāo)的標(biāo)準(zhǔn)化值作因變量，利用多元線性回歸分析得到每個因子的貢獻(xiàn)率。計算公式如下：

(8)

(9)

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

數(shù)據(jù)預(yù)處理和多元統(tǒng)計過程采用Excel 2016和SPSS 24軟件，空間分布展示采用ArcGIS 10.2軟件。

2 結(jié)果與討論

2.1 河流水質(zhì)污染特征

永定河張家口段各監(jiān)測斷面2014—2017年水質(zhì)指標(biāo)統(tǒng)計見表1。由表1可知，洋河、桑干河、清水河均受到不同程度的污染，CODMn、CODCr、BOD5和TP、NH3-N、F-濃度均出現(xiàn)超過GB 3838—2002中Ⅲ類水質(zhì)標(biāo)準(zhǔn)的情況。其中，洋河F-、TP污染較為顯著，F(xiàn)-濃度平均值超過Ⅲ類水質(zhì)標(biāo)準(zhǔn)，TP濃度最大值超過Ⅳ類水質(zhì)標(biāo)準(zhǔn)；清水河BOD5超標(biāo)較為顯著，最大值超過Ⅳ類水質(zhì)標(biāo)準(zhǔn)?？傮w來看，洋河、清水河的水質(zhì)污染程度較桑干河更嚴(yán)重。邵志江等[30]曾對永定河上游區(qū)域2013—2016年的監(jiān)測數(shù)據(jù)進(jìn)行統(tǒng)計，發(fā)現(xiàn)主要污染指標(biāo)為TP和F-，且洋河污染較嚴(yán)重，清水河污染較輕，而桑干河水質(zhì)清潔。

表1 研究區(qū)2014—2017年水質(zhì)指標(biāo)統(tǒng)計Table 1 Statistical description of water quality indicators in the study area during 2014-2017

2.2 河流水質(zhì)空間特征

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

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

圖3 研究區(qū)水質(zhì)指標(biāo)空間差異Fig.3 Spatial variation of water quality indicators in the study area

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

圖4 研究區(qū)水質(zhì)污染分區(qū)及工業(yè)企業(yè)分布Fig.4 Water pollution classification zonation results and industry distribution in the study area

2.3 水質(zhì)污染因子識別

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

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

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

A組提取出6個因子，累計方差為67.11%，分別提取出高荷載變量Zn、DO(F1)，F(xiàn)-、糞大腸菌群數(shù)(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 研究區(qū)土地利用現(xiàn)狀Fig.5 Land use in the study area

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

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

2.4 污染源貢獻(xiàn)率分析

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

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

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

2.5 污染源空間特征分析

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

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

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

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

3 結(jié)論

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

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