基于L-THIA模型与3S技术的大亚湾陆域非点源总氮污染研究

2019-12-19蒋婧媛徐姗楠黄洪辉刘华雪

应用海洋学学报 2019年4期

蒋婧媛,徐姗楠,黄洪辉,刘华雪

(1.广东省环境科学研究院，广东广州 510045； 2.中国水产科学研究院南海水产研究所，广东广州 510300；3.生态环境部华南环境科学研究所，广东广州 510530)

大亚湾是我国亚热带海域的重要海湾之一，湾内生物资源丰富、生境多样，设立了大亚湾水产资源自然保护区和国家级海龟自然保护区。大亚湾规划区于1992年被批准为国家经济技术开发区后，沿岸社会经济进入了快速发展时期，但伴随的是入海污染物通量的持续增长，导致大亚湾海域生态环境问题日益突出，生物资源和生态服务功能衰退，赤潮频发[1-3]。大亚湾海域的生态环境问题与陆域污染密切相关，包括点源污染和非点源污染。近年来，大亚湾沿岸区域建成了多家污水处理厂，点源污染控制力度日益加强。但非点源污染，即溶解的和固体的污染物从非特定的地点、在降水冲刷作用下通过径流过程汇入受纳水体并引起的污染，由于其具有广泛性、随机性、不确定性和难监测等特点[4]，防治措施仍然较为薄弱。随着入海污染物总量控制制度的实施，非点源污染的定量化研究逐渐引起广大学者的广泛关注[4-5]。由于非点源污染治理难度大，进行来源贡献的定量识别以及对关键源区的空间识别，有针对性地进行非点源污染控制管理具有重要意义[6]。

非点源污染的特点使得依靠监测对其进行定量分析评估存在难度大、投入高、耗时长等缺点，目前一般应用模型来模拟和分析非点源污染。非点源污染模型经历了从经验模型到机理模型、从纯数值模拟程序转向大型可视化专业软件的发展历程[7-8]。经验模型主要通过建立土地利用或径流量与非点源污染负荷之间的相关关系来估算污染负荷量，属于“黑箱”模型，数据处理方法简单，但由于缺乏机理基础，往往精度较低。机理模型围绕非点源污染的发生和迁移转化等具体过程构建，影响较大、广泛使用的有AnnAGNPS、HSPF、SWAT、SWMM等，属于“白箱”模型，通过参数的合理率定和校准能实现较准确的模拟预测，但参数众多，对数据量、数据精度要求高。模型的选择必须综合考虑研究目的和数据的可获得性[9]。

大亚湾陆域非点源污染研究尚处于起步阶段，数据累积少，在现阶段难以满足应用机理模型的数据需求。由美国环保署与美国普渡大学联合开发的L-THIA模型(Long-term Hydrologic Impact and Non Point Source Pollutant Model)相对SWAT等复杂机理模型对数据的要求较少，只需输入土地利用、土壤类型和长时间序列降雨数据就能模拟计算年均径流量和非点源污染负荷，能有效地考虑流域下垫面和降雨特点，并且很好地与GIS软件进行了结合，便于模型构建管理和模拟结果的可视化，对资料缺乏地区有很好的应用价值，在我国资料缺乏区域的非点源污染负荷模拟研究中得到了较广泛的应用[10-14]。因此，本研究选用L-THIA模型，并借助3S(GIS、RS和GPS)技术在模型参数确定和遥感信息之间建立直接的关系，同时利用入海河流通量观测结果等对关键参数进行率定，最后估算得到大亚湾陆域非点源总氮(TN)污染负荷量，并进行汇水区内来源贡献的定量识别以及对关键源区的空间识别。本研究一方面可为大亚湾落实污染物总量控制制度、制定污染治理对策等提供重要信息和决策参考，同时，也可为海湾地区采用易于推广的模型方法和容易获得的数据资料对非点源污染负荷进行定量估算提供可借鉴的方法。

1 研究对象与方法

1.1 研究区概况

大亚湾是一个典型亚热带溺谷型海湾(图1)，三面环山，湾口潮向东南，东接稔平半岛，与红海湾相接，西连大鹏半岛，与大鹏湾相邻，内、外湾面积共约1 200 km2。

图1 大亚湾地理位置Fig.1 Geographic location of the Daya Bay catchment

1.2 模型介绍

L-THIA模型的核心是SCS(Soil Conservation Service)径流模型[12-13]，能反映不同土壤类型、不同土地利用方式及前期土壤水含量对降雨径流的影响。该模型算法是基于集水区的实际入渗量与实际径流量之比等于集水区该场降雨前的最大可能入渗量(或潜在入渗量)与最大可能径流量(或潜在径流量)之比的假定基础上建立的，并假定集水区该场降水的初损为该场降水前潜在入渗量的0.2倍，即

(1)

式(1)中：Q为径流深(mm)，P为降水深(mm)，S为最大可能入渗量(mm)。

由式(1)可以看出：集水区的径流量取决于降雨量与该场降雨前集水区的潜在入渗量，而潜在入渗量又与集水区的土壤质地、土地利用方式和降雨前的土壤湿润状况有关。SCS 模型通过一个经验性的、综合反映上述因素的参数来推求S，即

S=25 400/CN-254

(2)

式(2)中：CN为径流曲线数。可见CN值越大，越容易产生径流。

在径流量模拟的基础上，L-THIA模型通过对不同土地利用类型赋予相应的事件平均浓度(Event Mean Concentration, EMC)数据，从而得到每一种用地类型的非点源污染负荷总量，其公式如下：

L=CEMC·Qn·k

(3)

式(3)中：L为非点源污染负荷(t/a)，CEMC为事件平均浓度值(mg/dm3)，Qn为年径流总量(m3)，k为单位转换系数。

1.3 数据处理

1.3.1 汇水区划分大亚湾汇水区的划分运用了美国ESRI公司开发的Arc/INFO软件的ArcHydro工具，利用数字地面高程(Digital Elevation Model，DEM)进行流域划分。本研究使用的是基于“先进星载热发射和反辐射计(ASTER)”数据计算生成ASTER GDEM数据产品，分辨率为30 m。通过洼地填充、流向计算、汇流累积量计算和河网提取等步骤，最后划分得到大亚湾汇水区范围(图1)，涉及3个区(县)共11个镇(街)，如表1所示，惠州市占大亚湾汇水区范围的78.68%，深圳市占21.32%。

1.3.2 土地利用分类本研究采用的遥感影像数据为下载自地理空间数据云平台的2017年Landsat 8 OLI_TIRS遥感影像，数据分辨率为30 m。利用IDRISI地理信息系统与遥感综合软件对遥感影像进行解译以获得大亚湾汇水区的土地利用图。

进行标准假彩色合成后，根据2017年实地考察数据建立各种土地利用类型的解译标志，在此基础上通过屏幕数字化建立遥感影像土地利用分类训练区样本。最后采用自组织映射(Self-organizing Map, SOM)神经网络方法进行监督分类。根据大亚湾汇水区土地利用类型，结合各地类非点源污染产污特征，确定进行土地利用分类的主要类型包括：林地、草地、园地、耕地、建设用地、滩涂水域等6种类型。

与传统分类方法相比，神经网络分类方法具有更强的非线性映射能力，能处理分析复杂空间分布的遥感信息，一般可获得更高精度的分类结果，因此在遥感影像土地利用分类中得到广泛应用。IDRISI中利用SOM进行土地利用监督分类首先是进行非监督聚类分类：输入层在接受输入样本之后进行竞争学习，SOM算法首先根据特征映射确定在输出空间中的最佳匹配或获胜神经元，神经元的突触权值向量可视为神经元指向输入空间的指针；所有权值向量在输入向量空间互相分离，在每个获胜神经元附近形成一个聚类区，各自代表输入空间的一类模式(即某种土地利用类型)，并将输入空间的样本模式类有序地映射在输出层上[15-16]。上述工作完成后，通过与用户提供的各土地利用类型训练区数据对比，程序再将非监督分类结果一一进行土地利用类型标识。大亚湾汇水区各土地利用类型空间分布见图2，其面积及所占比例见表2。

表1 大亚湾汇水区内行政区划Tab.1 Administrative division in Daya Bay catchment

续表1

图2 大亚湾汇水区土地利用分类Fig.2 Land use classification in Daya Bay catchment

表2 大亚湾汇水区土地利用分类结果Tab.2 Classification of land use in Daya Bay catchment

1.3.3 土壤水文单元分类根据土壤的渗透特征将土壤分为A、B、C、D 共4类土壤水文单元[10-14]，水文单元定义为相似暴雨和植被覆盖条件下具有相似径流潜力的一组土壤(表3)。L-THIA模型的土壤数据输入需将大亚湾汇水区内的土壤按水文条件归类，并在GIS软件中生成相应的土壤水文图。本研究采用的土壤数据来源于联合国粮农组织和维也纳国际应用系统研究所所构建的世界土壤数据库(Harmonized World Soil Database Version 1.1， HWSD)，中国境内数据源为南京土壤所提供的第二次全国土地调查1∶100万土壤数据。基于HWSD中的分层土壤砂粒、粘粒占比及美国农部制土壤质地分类数据，利用美国华盛顿州立大学开发的土壤水分特性软件SPAW(Soil-Plant-Atmosphere-Water)进行土壤水文单元的判断，结果如图3所示。

1.3.4 降雨数据模型所需的多年逐日降雨数据采用中国气象科学数据共享服务网中国地面气候资料日值数据集(V3.0)中的惠州市站点2000—2017年的逐日降雨量观测数据，并根据模型需要对数据进行处理和格式转换。使用长时间序列降雨数据驱动模型可以了解流域不同暴雨级别的平均径流及污染物产出状况，而不是某一暴雨的径流和污染物输出，因此，对流域规划和管理具有很好的指导意义[12]。

表3 土壤水文单元分类标准Tab.3 Classification standard for soil hydrologic groups

图3 大亚湾汇水区土壤水文单元分类Fig.3 Classification of soil hydrologic groups in Daya Bay catchment

1.4 参数校正及有效性检验

由于大亚湾海域不是完整的封闭流域，无法用整个汇水区的实测流量和通量数据来校正模型参数和检验模拟结果的合理性。研究表明，在缺乏实测资料的情况下，通过比对不同方法和模型的计算结果能够帮助验证模型结果的有效性[17]。因此，本研究采用如下方法进行参数校正及有效性检验：

1.4.1CN值校正及径流量模拟结果验证根据《惠州市水资源综合规划》[18]，大亚湾区域综合径流系数为0.61，由此计算得到多年平均径流量为9.47亿m3。以SCS-CN模型的CN值查算表为基础，结合土地利用和土壤类型数据，参考国内学者建立的CN值表[11,19]，根据径流系数法计算结果对CN值进行校正(表4)，最终模型模拟得到多年平均地表径流量与径流系数法计算结果相差8.97%，说明模型径流量模拟结果是可靠的。

1.4.2 EMC值校正及负荷量模拟结果验证 L-THIA模型中内置的EMC参数值源自北美地区的监测数据，由于不同地区自然和社会经济条件的不同，模型自带的EMC值与我国有较大差异，应根据研究区域实际情况进行调整。

表4 大亚湾汇水区CN值Tab.4 CN value in Daya Bay catchment

为此，首先基于大亚湾入海河流污染物通量的研究成果[20]，在大亚湾沿岸的3个区(县)各选择一条流域内人类活动相对较少、点源污染以生活源为主、全年排放情况较为稳定的入海河流(大亚湾经济技术开发区的岩前河、惠东县的大埔屯河、龙岗区的南涌河)，采用水文估算法[21]得到3条河流非点源TN入海通量。然后，本研究参考研究区有关研究成果[22-23]获得初步EMC进行模型运算，通过对比上述河流入海通量与模型模拟结果对EMC进行调整。EMC值校正结果如表5所示，河流入海通量与模型模拟值比对结果如表6所示，参考Moriasi等(2007)[24]提出的模拟结果满意程度分级标准，相对偏差在可接受范围内。

《惠州市水资源保护规划》[25]中对大亚湾经济技术开发区非点源TN负荷的估算结果为589.90 t/a，与本研究的模拟结果668.67 t/a相对偏差为13.35%，同样在可接受范围之内。

综上，可认为本研究所构建模型的模拟结果基本合理，可用于研究分析大亚湾陆域非点源TN污染。

表5 大亚湾汇水区各种土地利用类型TN的EMC取值Tab.5 TN EMC value of various land use types in Daya Bay catchment

表6 模型有效性检验结果Tab.6 Result of model validity test

2 结果与讨论

2.1 来源贡献定量识别结果

2.1.1 土地利用类型贡献量运用L-THIA模型模拟得到大亚湾陆域非点源TN负荷量，结果如表7和图4所示。大亚湾陆域非点源TN负荷量为2 559.00 t/a，其中耕地贡献量最大为871.44 t/a，占总量的34.05%；其次为建设用地(798.30 t/a)，占31.20%；再次为园地(543.36 t/a)，占21.23%；耕地、建设用地和园地合计贡献了86.48%的非点源TN负荷量。

从单位面积TN流失强度来看，耕地最高(145.30 kg/hm2)，其次为园地(52.91 kg/ hm2)和建设用地(51.92 kg/ hm2)，林地和草地的单位面积TN流失强度远远小于其他用地类型，分别为7.12、7.04 kg/ hm2。综合贡献量和单位面积TN流失强度来看，耕地对大亚湾汇水区内非点源污染的影响最大，TN负荷量和单位面积TN流失强度均为最高，其次是建设用地和园地，而林地、草地对污染的截留效果较好，流失强度低。

表7 大亚湾汇水区各土地利用类型非点源TN贡献与单位面积流失强度Tab.7 Contribution and loss intensity for non-point source TN load of various land use types in Daya Bay catchment

图4 大亚湾汇水区单位面积TN流失强度分布Fig.4 Distribution of TN loss intensity in Daya Bay catchment

2.1.2 行政区贡献量利用GIS的分区统计功能，计算得到汇水区范围内11个镇(街)的贡献量如表8所示。从各行政区非点源TN负荷量来看：①惠州市贡献了总量的86.29%，其中惠东县稔山镇和平海镇分别贡献了28.76%和23.92%，白花镇和平山街道仅有很小一部分区域位于汇水区范围内，因此贡献量非常小，惠东县合计贡献54.53%，为贡献最大的区(县)级行政区；大亚湾经济技术开发区贡献量占比为26.13%，分别来自霞涌街道(11.57%)、澳头街道(9.84%)和西区街道(4.72%)；惠阳区仅有淡水街道位于大亚湾汇水区内，贡献量占比为5.64%。②深圳市龙岗区贡献了大亚湾13.71%的陆域非点源TN负荷量，其中大鹏街道占5.43%、南澳街道占5.40%、葵涌街道占2.88%。

由于各镇(街)在大亚湾汇水区内的面积存在差异，故进一步采用单位面积TN流失强度进行分析。大亚湾汇水区内的平均单位面积TN流失强度为29.13 kg/hm2，在11个镇(街)中，有5个镇(街)单位面积TN流失强度超过汇水区平均值，均位于惠州市，按从大到小排序分别为淡水街道、白花镇、稔山镇、霞涌街道和西区街道。淡水街道因为在大亚湾汇水区范围内的建设用地和耕地合计所占比例高于其他镇(街)，单位面积TN流失强度在11个镇(街)中最大，为44.66 kg/hm2。深圳市的南澳街道单位面积TN流失强度最小，为15.44 kg/hm2，这是因为该街道位于大亚湾汇水区内的区域约70%的土地为林地，水土保持和污染截留能力较好。

表8 大亚湾汇水区各镇(街)非点源TN贡献与单位面积流失强度Tab.8 Contribution and loss intensity for non-point source TN load of the towns in Daya Bay catchment

2.2 关键源区空间识别结果

采用单位面积TN流失强度进行非点源TN流失关键源区的识别，运用自然断点法将流失强度划分为低、较低、中等、较高、高等5个等级(表9)。由表9可知，高TN流失强度的区域在大亚湾陆域汇水区内占6.70%，但贡献了汇水区33.69%的非点源TN负荷量，这部分区域是大亚湾汇水区TN流失的关键源区，是治理非点源污染的优先控制区。

从表10可以看出，分别有43.71%和33.99%的关键源区位于惠东县的稔山镇和平海镇内，同时通过上文行政区贡献量的分析结果可知这两镇的TN贡献量居汇水区前两位，因此需注意加强两镇的非点源污染防治。

从图5可见，关键源区有较明显的沿河分布特征，例如在惠州稔山镇的白云河、竹园河，平海镇的油麻园水，霞涌街道的霞涌河，深圳大鹏街道的龙岐河等入海河流沿岸流失强度较高，这与沿河区域人类聚集程度高，城镇开发建设、农业种植等人类生产生活行为对土地的扰动大有关。因此，在制定非点源污染防治方案时，应优先考虑在入海河流沿岸人类活动密集区采取措施减少非点源污染入海通量。