李 昌 文,徐 照 明,甘 拯,游 中 琼,马 强
(长江勘测规划设计研究有限责任公司,湖北 武汉 430010)
湖泊不同水位条件下的面积、容积是一项最基本的水文特征值[1]和重要的生态变量[2-3],是表征湖泊受气候变化与人类活动影响程度的重要指示器[4-5]。目前,关于湖泊的水位-面积(容积)关系构建方法已取得了大量研究成果。研究表明,利用遥感可快速提取湖泊面积,然而遥感数据并非实时获取、云层厚度影响精度、卫星偏角易产生畸变,存在时效性不强、时间分辨率不高等缺点,在低水位下的模拟值误差较大[6-8];MODIS遥感影像反演湖泊面积的精度较好,可为云雾天气下水体面积遥感定量监测提供弥补方法[9-12];通过DEM与遥感结合,可解决遥感影像受时空分辨率和云覆盖影响不能反演关键水位下的面积问题[13];利用水动力模型构建的实测水位与湖面面积的相关关系整体误差较小,但易受“涨退水”过程的影响,即同一水位下存在不同的湖面面积[8];在湖泊水面面积提取的基础上,可结合湖区水位站同期水位观测数据,建立统一的水位-面积关系[14]或刻画不同时期的水位-面积绳套关系[1, 15];传统方法基于地形图等高线,建立水位-面积关系,采用体积公式计算水位-容积关系,山区由于等高距较大,无法满足工程精度要求,平原地区由于地势平坦,等高线绘制困难带来较大误差,而基于DEM可解决该问题[16]。
洞里萨湖地处柬埔寨西北部,湄公河三角洲金边河段北岸,是东南亚最大的天然淡水湖泊,流域面积8.6万km2,被称为柬埔寨的心脏,是人民的“生命之湖”。洞里萨湖是湄公河的一个过水型、吞吐型、季节型湖泊,水位涨落、面积变化受控于湖区支流及湄公河的多重影响,入湖支流来水变化和湄公河干流水位变化的不同组合,造成洞里萨湖水位与面积年内、年际变化较大[17-18]。因此,合理构建洞里萨湖的水位-面积(容积)关系对洞里萨湖区与湄公河三角洲的防汛抗旱、调蓄能力评估、水资源开发利用和生态环境保护具有重要指导意义。
Sogreah[19]和柬埔寨农业、森林和渔业部水文司Sopharith T[20]分别于1966年和1997年基于1∶250 000 Sogreah地形图建立了洞里萨湖的水位-面积(容积)关系,Sogreah湖区地形图于1963年由3组航拍照片绘制而成,平面精度为100~250 m,垂向上有5条等高线,由于没有建立地面控制点进行几何精校正,精度较差,很难与其他地形图进行比较。1998年,Teng Peng Seang[21]基于湄公河委员会秘书处提供的洞里萨湖区Certeza测绘图和卫星影像数据分别构建了洞里萨湖水位-面积(容积)关系。Certeza 测图为1∶100 000地形图,于1964年完成,水平精度为65 m,垂向精度为1 m,包括1~13 m共计13条等高线,精度较差。卫星影像由于不能识别洞里萨湖区的大片洪泛森林和湖周的大片农田,造成中水位时模拟面积偏小(湖泊水体受洪泛森林遮挡),高水位时模拟面积偏大(湖周农田被误判划入了湖区),由于不能准确反演湖水深度,提取的洞里萨湖容积小于实际值[22]。2001年,Jantunen T[23]基于Certeza测绘图和更新的水文图集UHA构建了DEM和不规则三角网TIN,并基于Arc Info构建了洞里萨湖的水位-面积(容积)关系曲线,柬埔寨水文图集包括FINNMAP 1999年完成的洞里萨湖1∶100 000地形图和洞里萨河1∶20 000地形图,水平精度为1~2 m,垂直精度为0.2 m,计算精度得到了进一步提高。近几年,湄委会和柬埔寨国家湄公河委员会[24-25]结合1998~1999年的旱季水文调查成果、1993年的洞里萨湖洪泛平原Certeza测绘图、2003年的SRTM等数据,构建了更高精度的DEM,在此基础上建立了洞里萨湖水位-面积(容积)函数关系。
总体来看,洞里萨湖的水位-面积(容积)关系已有研究成果质量随着地形资料精度的提高逐步提升,但仍存在以下不足:这些关系均基于洞里萨湖在1,2,3,4,5,6,7,8,9,10,11 m等高线的面积和容积关系构建,其他高程下的面积和容积只能通过插值计算得出,插值计算成果与实际值有所差别,而洞里萨湖是开敞湖泊,湖水位的年内、年际变化较大,水位波动尤其是枯水位的较小变化均会造成湖泊面积、容积的大幅改变[26],进而对洞里萨湖的结构和功能产生重要影响。因此迫切需要构建更为精细的洞里萨湖水位-面积(容积)关系来反映湖区地形的实际变化情况。同时,洞里萨湖的水位-面积(容积)关系尚未与湖区控制性水位站的实测水位进行匹配,响应关系不明晰,不符合防汛抗旱精准预报及水资源综合管理的实际工作需求。
1.1.1采用的地形数据
本次研究收集了洞里萨湖区的陆上地形图、水下地形图和DEM资料。
(1)陆上地形。收集到柬埔寨全国1∶50 000陆上地形图(纸质),该图于20世纪60~70年代测绘,基本等高距为10 m,部分平坦地区有5 m间曲线。
(2)水下地形。收集到柬埔寨水文图集(UHA),UHA包括洞里萨湖1∶100 000、洞里萨河1∶20 000水下地形图,水平精度为1~2 m,垂直精度为0.2 m,该图于1992~1993年测量,1998~1999年更新,无等深线只有水深点。水下地形资料有纸质和GIS矢量数据两种,其中纸质地图有属性信息,水深单位图面注明为分米;GIS地图没有属性信息,洞里萨湖水深单位为厘米、洞里萨河水深单位为分米。经检查、核对,纸质和GIS两种资料为同一数据源。因此,本次研究采用GIS数据,相关属性信息采用纸质图。
(3)DEM。收集到芬兰环境研究所(Finnish Environment Institute)基于洞里萨湖区1993年1∶100 000 Certeza测绘图、柬埔寨水文图集和2003年SRTM数据构建的洞里萨湖湖区DEM,间距为100 m,格式为xyz文本文件。芬兰环境研究所提取的DEM数据相对1∶50 000陆上地形图+洞里萨湖1∶100 000水下地形图+洞里萨河1∶20 000水下地形图信息更全,精度更高,近年来常被湄公河委员会和柬埔寨水利气象部用来构建洞里萨湖的水位-面积(容积)关系[24-25]。
1.1.2采用的坐标基准
陆上地形和水下地形的平面基准分别为INDIAN DATUM 1960 UTM Zone 48N和INDIAN DATUM 1975 UTM Zone 48N,高程基准分别为M.S.L.Hatien datum和当地基面LLW,为与柬埔寨境内水文站网的坐标基准保持一致,本次研究将洞里萨湖区地形和河流水系、行政区划等有关数据的平面基准统一转换为WGS 1984 UTM Zone 48N,高程基准统一转换为M.S.L.Hatien datum。
1.1.3量算范围
洞里萨湖出湖控制站为波雷格丹站,实测最高水位为10.54 m(2011年10月20日),最低水位为1.11 m(2010年6月8日),结合洞里萨湖区地形以及湄公河委员会和柬埔寨水利气象部对湖区范围的界定情况,本次研究将洞里萨湖区的量算范围界定为波雷格丹站以上1~11 m等高线之间的区域。
1.2.1量算思路
基于方法一(采用1∶50 000陆上地形图+洞里萨湖1∶100 000水下地形图+洞里萨河1∶20 000水下地形图)和方法二(采用芬兰环境研究所DEM)分别构建洞里萨湖湖的水位-面积(容积)曲线。
首先,对洞里萨湖区地形资料进行矢量化和接合等处理,在此基础上基于ArcGIS构建不规则三角网TIN,量算不同水位下的洞里萨湖面积和容积。洞里萨湖湖区水位-面积(容积)曲线构建的技术路线见图1。
图1 洞里萨湖湖区水位-面积(容积)曲线构建技术路线
1.2.2量算方法
不规则三角网(TIN)能随地形起伏变化的复杂性而改变采样点的密度并决定采样点的位置,因而它能够避免地形起伏平坦时的数据冗余,又能按地形特征点如山脊,山谷线,地形变化线等表示数字高程特征。为了提高地形表达精度,运用ArcGIS来构建不规则三角网对洞里萨湖区域进行高程模拟。提取已有地形图高程点,等高线,水系特征线等关键要素生成TIN。不同的几何类型可以提供不同的表面要素类型,主要包括以下几个方面。
(1)水深点及高程点-离散多点。离散多点是TIN中的主要输入要素,由它们来决定表面的总体形状。通过离散点构造数字高程模型效果如图2所示。
图2 通过离散点构造数字高程模型
(2)洞里萨湖边线及双线河-水系特征线(隔断线)。隔断线通常用于呈现自然要素(如山脊线或河流)或建筑要素(如道路)。隔断线有以下2种:硬隔断线与软隔断线。隔断线可以有高程信息,也可以没有高程信息。硬隔断线:硬隔断线用于表示表面坡度的不连续性。河流和道路断面可作为硬隔断线包括在TIN中。硬隔断线能够捕获表面的突变并能改进TIN的显示和分析质量。软隔断线:软隔断线是不会改变表面局部坡度的线状要素,例如表示研究区范围边界的线等。
(3)量算范围-裁剪多边形。裁剪多边形用于定义TIN表面的边界。位于裁剪多边形之外的输入数据将被从插值和分析操作(例如,等值线或体积计算)中排除。
(4)不同水位下的洞里萨湖面积和容积量算。利用已建立的不规则三角网(TIN),根据实际形态特征将水体微分成若干个棱柱体,通过对每个柱体的体积求和,即可求得整个湖泊的容积。在此次运算过程中运用ArcGIS中3D Analyst中Surface Volumne工具进行计算分析,获取指定参考平面以下的TIN数据集表面的面积和容积。洞里萨湖不规则三角网(TIN)如图3所示。
图3 洞里萨湖不规则三角网(TIN)
1.2.3基础数据处理
(1)地形图矢量化。因收集到的1∶50 000地形图均为纸质数据,需要进行矢量化。
精度控制要求:矢量化依据CH/T 1015.4-2007《中华人民共和国测绘行业标准》基础地理信息数字产品1∶10 000,1∶50 000生产技术规程。具体精度要求为:图纸扫描分辨率不低于400 dpi,图廓定向点点位误差小于0.1 mm,线状地物采集误差一般小于0.2 mm,点状地物采集误差小于0.1 mm。
纸质地图扫描:图纸扫描采用卡莱泰克GX+38C彩色大幅面扫描仪进行扫描,扫描分辨率为400 dpi。
地图校正:由于纸质地图存在变形,且扫描得到的地图JPG文件缺少坐标定位,需要对地图JPG图片进行纠正及配准。首先将扫描好的影像数据进行地理几何纠正,生产带标准坐标的DRG数据,具体坐标系统为Indian 1960 UTM Zone 48N。为了消除图像局部变形,选取图廓四周的图廓点以及所有的公里网格为控制点进行纠正,纠正精度≤0.2 mm。
地图矢量化:洞里萨湖矢量化涉及1∶50 000地形图共计65幅,接合图如图4所示。矢量化软件平台采用Geoway 3.6,矢量化范围设定为洞里萨湖及其周边0~20 m等高线之间的区域,按精度≤0.2 mm对0~20 m高程范围内高程点、等高线、湖泊、双线河进行矢量化,部分不闭合的5 m和15 m间曲线保持其状态,生成shp格式地形数据。
图4 洞里萨湖1∶50000地形图接合示意
(2)水下地形数据处理。对洞里萨湖1∶100 000的4幅水下地形图和洞里萨河1∶20 000水下21幅地形图进行接合,如图5所示。
图5 洞里萨湖水下地形接合示意
高程系统转换:采用GIS数据,利用纸质水下地形图中的转换公式(LLW=Ha Tien MSL + 1.20),将水深点转换为HaTien平均海平面的高程点。每幅图水深点转换为HaTien平均海平面的高程点加常数不同。
本次洞里萨湖水位容积关系量算,采用了以下2种途径:
(1)途径一,柬埔寨1∶50 000陆上地形图+洞里萨湖1∶100 000水下地形图+洞里萨河1∶20 000水下地形图;
(2)途径二,芬兰环境研究所DEM数据。
上述2种途径与2006年湄委会WUP-FIN、2011年Huon Rath[24-25]构建的洞里萨湖水位-面积(容积)关系成果比较如图6和图7所示,可以看出,4套成果的水位-容积关系均较为接近。途径一构建的洞里萨湖的水位-面积关系曲线不光滑,在水位为4,10 m时出现了陡变,分析其主要原因是高水位部分采用的数据为1∶50 000陆上地形图,该图基本等高距为10 m,精度相对较差。途径二,湄委会WUP-FIN和Huon Rath基于芬兰环境研究所DEM数据构建的洞里萨湖水位-面积关系曲线较为光滑,结果极为接近。根据地形资料分析,洞里萨湖区地势非常平坦,水面比降小,多年平均水面比降为0.0 062‰,最大水面比降为0.0 189‰,水位-面积关系曲线不应出现突变或陡变,因此,途径二量算的洞里萨湖水位—面积(容积)关系更为合理。从这个角度也反映出了2002年以前有关洞里萨湖水位-面积(容积)关系研究成果的缺陷性。因此,本次研究选取途径二的量算结果作为最终成果。
图6 洞里萨湖水位-面积关系曲线研究成果对比
图7 洞里萨湖水位-容积关系曲线研究成果对比
为反映湖区地形的实际变化情况,满足防汛抗旱精准预报的实际工作需求,将途径二量算的水位-面积(容积)关系成果高程等间距由1 m细化为0.1 m,精细化的研究成果如图8和图9所示。从图8,9可以看出:不同水位下的洞里萨湖面积、容积变化较缓,具有“低水似湖、高水湖相”的自然景观特点,这与鄱阳湖的“低水河相、高水湖相”、“枯水一线、洪水一片”差异明显。
图8 洞里萨湖水位-面积关系精细成果
图9 洞里萨湖水位-容积关系精细成果
洞里萨湖的湖盆平坦,平均比降为0.002 4‰,洞里萨湖主湖区有甘邦隆1个水位站,该站距离洞里萨河河口约172 km,对洞里萨湖水位的分析具有较强的代表性。洞里萨河为洞里萨湖与湄公河之间的连接河道,波雷格丹站为洞里萨湖的出湖水文控制站,距洞里萨河河口约32 km。
本次研究依据收集到的1999年和2000年倒灌期洞里萨湖入湖、出湖站点实测流量资料,分析其水量变化,并结合湖区控制站甘邦隆的实测水位,对洞里萨湖水位-面积(容积)关系进行复核。具体复核方法如下:依据支流入湖控制站的实测流量资料以及支流流域面积与出口控制站集水面积的2/3次方计算支流入湖径流;依据波雷格丹的实测流量资料计算湄公河倒灌入湖水量;根据有关研究成果[24, 27],按波雷格丹站实测流量的5%近似计算漫滩入湖水量,根据洞里萨湖区洪泛平原地形,漫滩时机考虑为水位高于8 m;依据洞里萨湖出、入湖水量计算成果,合成入湖洪水过程,分析湖容的变化过程;依据洞里萨湖水位-容积关系和湖容的变化过程计算洞里萨湖的水位变化过程(以下简称模拟水位);比较甘邦隆站的实测水位及模拟水位,分析洞里萨湖水位-面积(容积)关系的合理性。由于本次研究未收集到湖区气象站点的实测降雨、蒸发资料,因此暂不考虑湖面降雨量和蒸发量。
1999年和2000年湄公河倒灌期间洞里萨湖的实测与模拟水位过程如图10所示。由图10可以看出:实测与模拟的水位过程非常接近,1999年绝对误差最大为0.21 m,平均为0.06 m,相对误差最大为9.9%,平均为1.5%;2000年绝对误差最大为0.22 m,平均为0.03 m,相对误差最大为5.1%,平均为0.1%。由此可认为,本次研究构建的洞里萨湖水位-面积(容积)关系是合理的。
图10 倒灌期内洞里萨湖的实测与模拟水位过程
目前,洞里萨湖湖区防汛抗旱水文情报预报以洞里萨湖的甘邦隆水位站和洞里萨河的波雷格丹水文站为代表,故以甘邦隆站和波雷格丹的水位为代表,建立洞里萨湖水位-面积(容积)关系与实测湖水位的响应关系。
根据洞里萨湖水位-面积(容积)关系复核成果,洞里萨湖水位-面积(容积)关系中水位数据直接采取甘邦隆站实测水位是合理的。因此,本次重点分析洞里萨湖水位-面积(容积)关系与波雷格丹站水位的响应关系。
甘邦隆站位于洞里萨湖湖区内,其与波雷格丹站的距离约为140 km,受复杂河湖水情及地形的扰动,其与波雷格丹站水位呈绳套曲线,如图11所示。可以看出,45°线以上对应为湄公河向洞里萨湖倒灌期,即河水位高于湖水位,甘邦隆站与波雷格丹站水位关系点据局部相对散乱,总体相关度较好,相关系数为0.947 6;45°线以下对应为汛后洞里萨湖向湄公河补水期,即湖水位高于河水位,两站的水位关系点据较为密集,相关度较好,相关系数达到0.964 4。
图11 甘邦隆站与波雷格丹站的水位相关关系
根据水位-面积(容积)关系及甘邦隆站1999~2011年逐日水位,计算得到洞里萨湖的逐日面积和容积值。
点绘补水期间洞里萨湖面积、容积与波雷格丹站水位的相关关系点群图,见图12。可以看出,洞里萨湖面积与波雷格丹站水位、洞里萨湖容积与波雷格丹站水位的相关关系均较好,相关系数分别为0.988 7和0.990 7。
图12 补水期间波雷格丹站实测水位与洞里萨湖面积和容积的相关关系
考虑到湄公河倒灌入湖流量与河湖水位差密切相关,且随着河湖水位差的增加而增加[17],而倒灌水量的变化最终引起湖面面积和湖容的变化。因此,点绘倒灌期不同河湖水位差(洞里萨河河口金边港站与洞里萨湖甘邦隆站的水位差)下的洞里萨湖面积、容积与波雷格丹站水位的相关关系点群图,如图13和图14所示。可以看出,在相同河湖水位差条件下,洞里萨湖的面积和容积与洞里萨河水位的相关关系较好,相关系数达到0.98以上,且随着洞里萨河水位的升高而增加。
图13 倒灌期间波雷格丹站实测水位与洞里萨湖面积的相关关系
图14 倒灌期间波雷格丹站实测水位与洞里萨湖容积的相关关系
基于本次收集到的不同来源的洞里萨湖湖区地形资料,在进行矢量化和接合、高程系统转换等处理的基础上,基于ArcGIS构建了不规则三角网TIN,量算得到不同水位下的洞里萨湖面积和容积。综述了洞里萨湖水位-面积(容积)关系的已有研究成果,并与本次成果进行了对比分析,评估出了最为合理的洞里萨湖水位-面积-容积关系。依据水文学方法采用洞里萨湖入湖支流站点与洞里萨河出湖站点的实测流量资料,分析其水量变化,并结合湖区控制站甘邦隆的实测水位,对洞里萨湖水位-面积(及容积)曲线进行了复核分析。结果表明,本次研究构建的洞里萨湖水位-面积(容积)是合理的。
以洞里萨湖出湖控制站波雷格丹站为代表,建立了洞里萨湖水位-面积(容积)关系与实测湖水位的响应关系。结果表明:洞里萨湖面积、容积与波雷格丹站水位呈绳套曲线,绳套两侧分别对应汛期湄公河向洞里萨湖倒灌期和汛后洞里萨湖向湄公河补水期,结合河湖关系研究成果构建了不同时期、不同河湖水位差下的洞里萨湖面积、容积与波雷格丹站水位的相关关系。
本次研究基于实测地形资料构建的洞里萨湖水位-面积(容积)关系,从精度上达到了现有条件下的最高水平,避免了基于遥感技术反演的精度问题,刻画了倒灌、补水等不同时期的水位-面积(容积)绳套关系,剖析了不同时期及河湖水位差对面积和容积的影响,同时精细化的量算成果不仅满足了防汛抗旱精准预报的实际工作需求,还为下一步构建湄公河水文-水动力学-水质耦合模型奠定了基础。