APP下载

汕头港海域氮、磷营养盐环境容量及排放总量控制的研究

2012-06-06张静柯东胜方宏达孙省利

大连海洋大学学报 2012年3期
关键词:梅溪榕江沙河

张静,柯东胜,方宏达,孙省利

(1.广东海洋大学水产学院,广东湛江524025;2.国家海洋局南海环境监测中心,广东广州5103003.广东海洋大学海洋资源与环境监测中心,广东湛江524025)

汕头港位于广东省东部沿海 (116°38'~116°45'E,23°19'~23°21'N),属河口港,是中国主要港口之一,担负着粤东、闽西南、赣南地区对外贸易进出口货物的运输[1]。榕江接纳汕头市沿岸地区生活污水和工业废水,由西向东经牛田洋注入汕头港;韩江西溪叉道梅溪直接注入汕头港,新津河和外沙河经汕头港东侧流入汕头海区。此外,汕头市主要陆源入海排污口有8个,其中直接排入汕头港的有6个。监测结果显示,经河流入海及排污的主要污染物为石油、悬浮物、氨氮、磷酸盐和 COD[2-4]。目前,关于汕头港海域的研究较少,主要有水温、盐度和三氮营养盐的分布以及三维流场数值模拟和温排水热扩散三维数值模拟的研究[5-7],但未见对氮、磷营养盐环境容量的研究。为此,本研究中,作者根据2006—2007年度广东省“908”专项水体环境和专项物理海洋调查资料和各相关验潮站的验潮资料,以及2006年8月位于研究海域内的贝类养殖区1个站位和汕头全海域海洋环境质量监测5个站位的监测资料,采用数值模拟方法计算了汕头港海域氮、磷营养盐的环境容量,并通过模拟方法给出了排放总量的控制方案,为确定该海域氮、磷营养盐的总量控制与管理计划的实施提供参考。

1 模型描述

1.1 水动力模型

1.1.1 模型概述 采用Blumberg和Mellor开发的三维原始方程海洋环流模型 (princeton ocean model,POM)为水动力模型。鉴于汕头港沿岸具有广阔潮间带的特征,采用POM08版的干湿网格方法来模拟实际的潮汐涨落情况,同时在程序中添加了河流边界的输入,并考虑了河流流量输入对研究海域水位的贡献。此外,在ADVCT子程序中采用抑制数值耗散的一阶精度迎风格式以增加模式计算的稳定性。模型的基本特点、方程、边界条件及初始条件可参考POM使用手册[8]。

1.1.2 河流输入方式 将河流视为点源,将河流初始带来的流量和污染物平均添加到河流入海位置所在模拟区域的边界内的若干个水平网格面积上,并考虑河流输入对整个模拟海域水位的贡献[9-12]:

其中:Q为河流的径流量;A为河流点源所在网格的水平面积。

1.2 水质模型

在近岸海域中,水深一般较浅,污染物在水体中混合充分,没有明显的分层变化。容量研究往往考虑的是一个时期内污染物的平均状况,因此本研究中采用垂向平均的水平二维对流扩散方程来计算汕头港海域的污染物浓度。控制方程如下:

其中:Ci为污染物浓度,下标i依次代表无机氮和活性磷酸盐;u、v分别为x、y方向的速度分量(此处为上述三维流场计算的各层流速的垂向平均);Dx、Dy分别为x、y方向的涡动扩散系数,由Smagorinsky(1963)公式[8]给出,在实际计算中由POM模型实时地更新对流扩散模式中的垂向平均流速和涡动扩散系数;S0为源强,且

式中:Cs0为污染物的点源强,在数值计算中作为边界条件给出;Csi为非保守物质的衰减项;k为衰减速率,对于保守物质而言其值为0。

1.3 环境容量模型

1.3.1 环境容量的几个概念 环境容量又称为环境的承受力、环境的承载力、环境的忍耐力,在不造成环境不可承受的前提下,环境所能容纳某物质的能力[13]。根据这一定义,本研究中提出几个容量概念:1)天然环境容量,指假设污染物背景浓度为0(即假设水体完全不受污染)时,在特定的水质目标下所求出的各污染源污染物的允许排放量;2)背景环境容量,指当考虑背景浓度时在特定的水质目标下所求出的各污染源污染物的允许排放量;3)现状环境容量,指当考虑海域水体现状浓度时在特定的水质目标下所求出的各污染源污染物的允许排放量。

1.3.2 容量计算方法 目前容量计算方法主要有模型试算法、分担率法、分区达标控制法、地统计学和GIS方法[14]。本研究中采用分担率法,其主要步骤及计算公式参考文献 [15]。

2 模型设置与参数的选取

2.1 模拟区域及站位布设

研究区域地理概况及调查站位如图1所示。为了利用908专项调查中的调查资料,本研究中选取整个图1的区域作为研究的计算海域,两条虚线分别为计算域的东、南开边界。由于实际调查的站位和设计站位相差较大,因此,图中站位编号仍采用设计时的编号,但以不同的符号分别表示水文调查站位和水质调查站位。其中水质站位A8位于计算区域外,作为计算参考。

图1 汕头港地理位置及调查站位图Fig.1 Geographical position and station graph at Shantou Port

在容量研究中,根据2006年8月—2007年11月对汕头港海区4次水质调查的结果,取全年的平均值作为计算现状容量的背景浓度。鉴于此次专项调查布设站位都位于汕头港港湾以外,因此,同时取2006年8月位于研究海域内的贝类养殖区1个站位GD04和全海域海洋环境质量监测5个站位YQ001~YQ005的监测结果作为计算现状容量的背景浓度。计算中将氮、磷看作保守物质,从环境管理的角度讲,这是偏于安全的。

2.2 河流边界条件

本研究中主要考虑榕江、韩江西溪的叉道梅溪、新津河和外沙河4条河流的输入。入海污染物总量参考2006年和2007年《汕头市海洋环境质量公报》的数据。直接影响计算海区的市区排污口有6个,自西向东分别为新兴路、市委前、大华路、石炮台、特区码头和龙眼南路排污口。由于6处排污口相距较近,为了便于计算,将排污口合并,即将前4个排污口和后2个排污口各合并为一个污染源。因此,整个模拟中综合考虑4条河流和两个合并排污口的排污,作为河流边界输入。污染源及氮、磷污染物的年入海通量如表1所示。

表1 汕头港入海河流的流量及污染物的入海通量Tab.1 River flow and pollutants fluxes to Shantou Port

2.3 模式安排

汕头港动力场的计算采用POM08二维外模式,该模式在深圳湾的流场模拟中得到了很好的验证[16]。模拟计算的平面网格数为117×75,空间步长Δx=Δy=300 m,垂直方向取5个sigma层,各层厚度平均分布。以M2、S2、K1、O14个分潮和榕江、梅溪、新津河、外沙河共4条主要河流及两个合并排污口的年平均流量分别作为外界强迫源驱动程序。分潮的周期振荡分别从南开边界和东开边界加入,调和常数由妈屿岛、汕头验潮站的实测资料经调和分析得到。风速风向参考广东省东部海湾[1]的相关内容,汕头港夏季受西南季风影响,春、秋、冬季均以东北季风为主,平均风速为2.7 m/s。海面拖曳系数 CD通过 Large and Pond(1981)的方法计算得到:

经调试,外模时间步长取2 s,内模时间步长取60 s;底粗糙度z0=0.002;Horcon参数C=0.1;垂向背景黏度umol=5×10-5m2/s。计算的起始时间为2006年8月1日24:00时,连续积分32 d,其中8月13日 (A7、A8号站观测时间均为8月13日15:00至8月14日15:00)、8月19日 (A9号站观测时间为8月19日11:00至8月20日11:00)和8月21日 (A6号站观测时间为8月21日10:00至8月22日10:00)的模拟结果被用来和海流调查站位上的实测数据进行比较。

2.4 海域环境功能区域划分及水质控制点的选取

根据1999年7月印发的《广东省近岸海域环境功能区划》(粤府办[1999]68号)、《关于调整汕头市近岸海域环境功能区划有关问题的复函》(粤办函[2005]659号)及附《汕头市近岸海域环境功能区排污混合区划表》,涉及到研究海域的汕头港近岸海域环境功能区共6个,功能区涉及港口、排污、增殖、养殖和旅游,水质目标执行海水水质标准(GB3097-1997)[17]二类和三类海水水质标准。排污混合区共9个,结合混合区宽度、长度、水质保护目标和污染源所处位置,本研究中选取了12个水质控制点,控制点编号及需要满足的水质标准如表2所示。

表2 汕头港研究海域控制点及其水质标准Tab.2 Water quality control points in the coastal Shantou Port

3 结果与分析

3.1 流场的模拟

3.1.1 潮位验证 由图2可见:本模式所模拟的水位变化与实测值基本一致。其中,A8、A9号站位的实测值均不是很理想,总体趋势虽有,但监测中每小时的水位变化波动较大,因而导致了实测值和模拟值偏差较大。从地形上说,A6号站靠近莱芜岛,A8号站靠近表角,地形曲折变化较大,因此测量误差较大;而A9号站位于妈屿岛的背面,受岛屿阻隔的影响,外海潮波在传入后要发生折转变形,因此误差也较大;模拟结果最好的是A7号站,处在湾口较开阔处,基本反映了外海潮波传入后的情况。由图2还可见:汕头港海域的涨潮历时略大于落潮历时,观测期间最大潮差出现在A6号站,为2.4 m左右,其余站位潮差基本在1.6 m左右;在一日内出现两次高潮和两次低潮,潮高和潮时日内不等,可以判断汕头港海域的潮汐属不正规半日潮。根据潮汐系数:

对输出的一个月的潮位结果采用T-tide进行调和分析,得到汕头港海域的F=1.14,可进一步证明汕头港海域的潮汐属不正规半日潮。

3.1.2 潮位验证 本次调查获得的实测数据是垂向三层 (海面下1 m,0.6×水深,离海底1.5 m)的流速,为了增加二者的可比性,在数据处理时,对三层数据进行了矢量加权平均处理,相同的方法也用在对实测流向的处理上。

图2 A6~A9号站水深实测值和模拟值的比较Fig.2 The measured and simulation water depth comparisons at stations from A6 to A9

从图3可见,流速和流向的模拟值与实测值变化趋势基本吻合。其中,A6号站模拟的个别时刻流速大小偏差较大,但流向基本吻合;A7和A8号站模拟的个别时刻流速偏大,除个别时刻流向偏差较大外,其余各时刻均吻合良好;模拟最差的是A9号站,由于处于妈屿岛的背面,外海潮波传入时发生反射变形,而这个过程在模型中无法体现,因此流速偏差较大,从实测的流速也可以看出实测值所反映的变化趋势与其它站位相差较大,而流向的模拟结果吻合较好。总的来说,此模拟结果还是令人满意的,基本反映了汕头港海域的水动力情况。从模拟和实测的流向可以看出,汕头港海域潮流属典型的往复流性质,潮流受地形约束在A6号站附近海域基本呈西南-东北走向,而在其余站位海域基本呈西北-东南走向。从水位图和流速图的对比来看,汕头港落潮流速一般大于涨潮流速。

图3 A6~A9号站垂向平均流速和流向实测值与模拟值的比较Fig.3 Vertical average velocity comparisons between measured and simulation values at stations from A6 to A9

3.2 氮、磷营养盐环境容量的计算

3.2.1 响应系数场 响应系数场给出的是各污染源在单位源强条件下污染物质随潮变化的时变特征,反映的是各污染源所在区域的水体交换能力。为了获得相对稳定的响应系数场,本研究中仅以外海带入的氮污染物浓度作为背景浓度,查看每15 d后的浓度场 (每日平均)的变化情况。结果表明,整个汕头港海域水体在30~45 d后基本完成交换。为了保证计算的可靠性,响应系数场取第45天的模拟结果进行讨论。

在不考虑背景浓度和开边界浓度,仅以单位源强输入的条件下,汕头港海域的6个入海污染源中,榕江排污主要对牛田洋和汕头港西部海域的影响较大,由单位源强流入后至牛田洋海域,氮浓度下降了约90%。由于地形的影响,港内狭窄,水流交换不畅,且牛田洋又是增养殖区,因此,如果榕江排污过大的话将直接影响港内水环境质量。梅溪排污主要影响牛田洋东部至汕头港西部海域,除局部很小海域氮浓度较高外,影响程度基本在5%左右。新津河排污主要影响范围向西到达整个牛田洋养殖区,向东至汕头港口门外到达濠岛的东部海域,除入海排污处,氮浓度的较高值出现在汕头港西部海域的狭窄地形处,影响程度基本在10%左右。外沙河排污对全海域有较大影响,西至牛田洋养殖区,东至汕头港港外的大部分海域,就单位源强而言,影响基本在5%左右,只在外沙河入海的局部海域影响程度达10%以上。合并排污口一对牛田洋养殖区至汕头港西部海域有较大影响,就单位源强而言,这种影响比梅溪排污的影响更大,但影响程度最大不超过8%。合并排污口二从影响的强度和范围来说都较合并排污口一轻,影响程度最大不超过4%。

3.2.2 污染责任分担率 图4、图5分别给出了在海域背景浓度和开边界浓度为0、各污染源实际源强输入的条件下,模型运行45 d,汕头港海域氮、磷污染责任分担率的计算结果。

图4 汕头港海域各污染源无机氮污染责任分担率场Fig.4 Source-sharing rate contours of inorganic nitrogen of five pollutant sources in the coastal Shantou Port

图5 汕头港海域各污染源活性磷酸盐污染责任分担率场Fig.5 Source-sharing rate contours of active phosphate of five pollutant sources in the coastal Shantou Port

由图4可见,榕江排污在牛田洋西部海域氮污染责任分担率达到40%以上,东部海域氮污染责任分担率较小;梅溪排污在牛田洋、汕头港海域以及湾口的南部海域,氮污染责任分担率在20%以上,其余海域不足10%;新津河的污染责任分担情况与梅溪较相似,只是对全海域的影响,即分担率要更大一些;外沙河氮污染对全海域的影响最大,对湾口海域的氮污染责任分担率达到了60%以上,尤其是湾口北部海域,氮污染责任分担率高达80%以上,即使对港内的牛田洋海域,氮污染责任分担率也达到20%,表明外沙河排污严重影响了整个汕头港海域氮的含量;对于两个合并的排污口,排污基本只影响排污口附近海域很小的范围,对整个汕头港海域最大的氮污染责任分担率分别不足20%和10%,基本不影响海水的使用功能。

由图5可见,虽然污染源强不同,但各污染源对整个汕头港海域的磷含量的影响范围和趋势基本与氮相同。所不同的是对于全海域磷的责任分担,两个合并排污口对磷的分担率值较氮的分担率值要高一些,分别达到34%和15%,表明这两个排污口磷输入的强度比氮大。分担率场的计算结果表明,汕头港海域主要的污染来源还是陆源河流排污,其中榕江对港内海域,外沙河对整个汕头港海域,以及梅溪、新津河对牛田洋以东、汕头港口以西海域水质,均有着严重的影响,因此,在氮、磷营养盐的总量控制中,应重点关注榕江、外沙河、梅溪和新津河各污染物的入海通量。

3.2.3 各污染源天然容量、背景容量和现状容量(允许排放量) 由于污染责任分担是在开边界浓度和背景浓度为0的基础上算得的,因此当某污染源靠近开边界的时候,就会减轻该污染源对整个海域造成的水质影响,尤其当开边界 (即外海)某污染物浓度较高时,这种影响会更大。由表3可知:在汕头港海域的6个入海污染源中,除外沙河无机氮尚有很小的现状容量外,其余各污染源均无排放容量,在实际操作中,6个污染源的无机氮均应采取严格的削减措施;6个污染源的活性磷酸盐均有现状容量,从数值上看,除新津河活性磷酸盐的现状容量占其天然容量的比例较低,约为50%外,其余各污染源的现状容量占其天然容量的比例都较高,再次表明各污染源附近海域活性磷酸盐的现状浓度较低,只要采取一定的控制措施,即可保证海域活性磷酸盐的含量,满足水质标准。

4 氮、磷总量的分配

4.1 各污染源削减量和削减率的计算结果与分析

根据计算得到的汕头港海域6个入海污染源的现状容量 (允许排放量)和汕头港海域6个入海污染源的污染物入海通量,计算各污染源的削减量和削减率 (表3)。从表3可看出,在汕头港海域的6个入海污染源中,对于无机氮,只有外沙河有少许的允许排放量,但其值远小于现有排放量,基本需要100%削减;而其余5个污染源均需100%削减。对于活性磷酸盐,各污染源均尚有容量,除梅溪和外沙河需大幅削减,最低削减率达到83%以上外,其余各污染源也都需削减,削减率在23%~59%。计算结果表明:整个汕头港海域氮污染严重,现有监测质量已不满足规定的水质保护目标;整个汕头港海域磷污染相对较轻,采取一定的减排措施可以达到规定的水质保护目标。

表3 汕头港海域各污染源的容量及污染物的削减量和削减率Tab.3 Capacity of six pollutant sources,pollutants elimination amount and rates of five pollutant sources in the coastal Shantou Port

根据计算结果,汕头港海域的主要减排污染物为无机氮。因此,在设计减排方案之前需要先计算当各污染源均需100%减排时,在现状监测的背景浓度和开边界浓度条件下,各水质控制点达到水质标准所需的时间。值得一提的是,由于100%减排后相当于洁净的淡水流入,因此各污染源附近的水质控制点很可能首先达标而远离污染源的海域尚有未达标的区域。为此,根据功能区划,在水质超标区补充12个水质控制点,以增加计算的可靠性。

4.2 无机氮减排方案的设计

对于无机氮,依次设计8个减排方案,待浓度场稳定后结果如下:

方案一 (达标):各污染源均100%减排 (现有开边界浓度),待浓度场稳定后全海域水质可以达到相应的水质保护目标。

方案二(不达标):各污染源均90%减排(现有开边界浓度),当进一步降低减排率后,全海域无机氮浓度在一年后仍未达到各自的水质保护目标。

方案三 (不达标):各污染源均90%减排 (南开边界浓度满足一类水质标准),由图4-(b)可见,南开边界带来的无机氮海水浓度较高,对计算域内无机氮的浓度产生了一定的影响,因此设计南开边界无机氮浓度满足一类水质标准 (即为0.20 mg/L)运行后,全海域无机氮浓度在一年以后仍未达到各自的水质保护目标。可见,除了南开边界附近海域浓度稍有不同外,南开边界的改变并没有使汕头港内的无机氮浓度有明显变化。

方案四 (不达标):榕江90%减排,其余各污染源100%减排,除牛田洋海域达不到二类养殖功能外,其余各海域水质均可满足相应的水质保护目标。表明榕江排污即使大幅减排 (减排达90%),仍然达不到既定的水质目标,而只有100%减排后可以达到水质目标。

方案五 (达标):榕江100%减排,其余各污染源90%减排,这种减排方案能够使海域水质满足其相应的水质功能。

方案六 (不达标):榕江100%减排,其余各污染源80%减排,在汕头港内部达不到相应的水质保护目标。

方案七 (达标):榕江100%减排,外沙河70%减排,其余各污染源90%减排,这种减排方案能够使海域水质满足其相应的水质功能。

方案八(不达标):榕江100%减排,外沙河60%减排,其余各污染源90%减排,在达濠岛近岸海域水质超标,达不到相应的水质保护目标。

4.3 活性磷酸盐减排方案的设计

由于汕头港海域各个污染源活性磷酸盐均尚有容量,因此从实际操作和计算简便的角度出发,在容量计算的结果上,以榕江、梅溪和外沙河、其余污染源共3个分组设计减排方案。依次设计6个减排方案,待浓度场稳定后结果如下:

方案一 (达标):榕江20%减排,梅溪和外沙河均30%减排,其余污染源执行容量计算结果的减排率。可以看出,榕江带来的活性磷酸盐在到达牛田洋海域后浓度即降到很低的水平,牛田洋海域的活性磷酸盐浓度主要受其它污染源的影响。

方案二 (不达标):榕江20%减排,梅溪和外沙河均20%减排,其余污染源执行容量计算结果的减排率,全海域活性磷酸盐浓度不能达到相应的水质保护目标。影响海域水质的主要污染源是梅溪,其次是外沙河。

方案三 (达标):榕江10%减排,梅溪和外沙河均30%减排,其余污染源执行容量计算结果的减排率,这种减排方案能够使海域水质满足其相应的水质功能。

方案四 (不达标):榕江10%减排,梅溪和外沙河均20%减排,其余污染源执行容量计算结果的减排率,全海域活性磷酸盐浓度不能达到相应的水质保护目标。

方案五 (达标):榕江10%减排,梅溪和外沙河均30%减排,其余污染源10%减排,这种减排方案能够使海域水质满足其相应的水质功能。

方案六 (不达标):榕江10%减排,梅溪和外沙河均30%减排,其余污染源不减排,全海域活性磷酸盐浓度不能达到相应的水质保护目标。表明靠近港内的新津河、港内的两个合并排污口对海域水质达标还是有一定的影响。

4.4 氮、磷减排结果的讨论

从环境容量的计算结果可以看出,整个汕头港海域无机氮超标严重。通过设计的8个减排方案可以看出,榕江100%减排可以达标,但90%减排时不能达标;外沙河70%减排可以达标,但60%减排时不能达标;其余污染源影响较小,均设计为90%减排。因此,在实际执行时,榕江至少要减排95%,外沙河减排65%,其余污染源减排90%或适当降低减排率,才能使整个海域无机氮浓度基本满足相应的水质保护目标。

整个汕头港海域活性磷酸盐含量较低,只有局部海域 (主要是港内)超标,达不到相应的水质保护目标。从设计的6个减排方案中可以看出,榕江10%减排可以达标,但不减排不能达标;梅溪、外沙河30%减排可以达标,但20%减排时不能达标;其余污染源10%减排可以达标,但不减排不能达标。因此,在实际执行时,榕江至少要减排5%,梅溪、外沙河减排25%,其余污染源减排5%,才能使整个汕头港海域活性磷酸盐浓度基本满足相应的水质保护目标。

[1]中国海湾志编纂委员会.中国海湾志:第九分册:广东省东部海湾[M].北京:海洋出版社,1998.

[2]汕头市海洋与渔业局.2006年汕头市海洋环境质量公报[EB/OL].http://sdinfo.coi.gov.cn/hygb/dfhygb/2006/shantou

[3]汕头市海洋与渔业局.2007年汕头市海洋环境质量公报[EB/OL].http://www.ydtz.com/news/shownews.asp?id=32551

[4]汕头市海洋与渔业局.2008年汕头市海洋环境质量公报[EB/OL].http://www.soa.gov.cn/hyjww/hygb/yhsshyhjzlgb/lbn/webinfo/2009/05/1238639581189782.htm

[5]陈浩如,彭云辉,李少芬.汕头港沿岸水体温、盐度及三氮营养盐的分布特征[J].南海研究与开发,1990(3):14-20.

[6]黄平.汕头港水域潮流运动的三维数值模拟[J].海洋环境科学,1995,14(4):9-15.

[7]黄平.汕头港水域温排水热扩散的三维数值模拟[J].海洋环境科学,1996,15(1):59-65.

[8]Mellor G L.Users guide for a three-dimensional,primitive equation,numerical ocean model[M].New Jersey:Princeton University,2004.

[9]Kourafalou V H,Oey L Y,Wang J D,et al.The fate of river discharge on the continental shelf:1.modeling the river plume and inner shelf coastal current[J].Journal of Geophysical Research,1996,101(2):3415-3434.

[10]Kourafalou V H,Thomas N L,Sun J,et al.The fate of river discharge on the continental shelf:2.transport of coastal low salinity waters under realistical wind and tidal forcing[J].Journal of Geophysical Research,1996,101(2):3435-3455.

[11]Kourafalou V H.River plume development in semi-enclosed Mediterranean regions:North Adriatic Sea and Northwestern Aegean Sea[J].Journal of Marine Systems,2001,30:181-205.

[12]刘浩.渤海生态系统关键物理生物过程的数值研究[D].青岛:中国科学院研究生院,2005.

[13]GESAMP.Environmental Capacity:An Approach to Marine Pollution Prevention[J].Unepreg Seas Rep Stud,1986,80:62.

[14]朱静,王靖飞,田在峰,等.海洋环境容量研究进展及计算方法概述[J].水科学与工程技术,2009(4):8-11.

[15]张存智,韩康,张砚峰,等.大连湾污染排放总量控制研究:海湾纳污能力计算模型[J].海洋环境科学,1998,17(3):1-5.

[16]Zhang J,Lin J G,Sun X L,et al.Application of an Improved Wetting and Drying Scheme in POM[C]//2009 International Conference on Energy and Environment Technology,ICEET 2009,Guilin,China,2009:427-432.

[17]中华人民共和国国家标准.海水水质标准GB3097-1997[S].北京:中国标准出版社,1997.

猜你喜欢

梅溪榕江沙河
刘静设计作品(二)
云雾飘渺景如画
沙河板鸭营销策划方案
白沙河
白沙河
婷婷聊古诗
捏只可爱的哈士奇
榕江“萨玛节”
快乐童年——榕江侗族儿童
藏在深山中的蜡花
———榕江苗族蜡染