小流域坝系监测的内容决定着是小流域坝系监测的方法。小流域坝系监测方法从目前的研究现状可以看出,监测方法又可以分解为常规方法和现代新技术方法。对于常规方法我们在就简单的罗列出其方法,而现代新技术较为详细的研究。现代新技术在本项目中主要为GIS、GPS、RS的3S技术研究。
(1)工程基本情况监测
根据调查统计以及实地查看,对流域淤地坝的分布、数量、面积、配置比例等基本情况进行监测。
(2)坝系安全监测
在项目实施期间对流域22座淤地坝的坝体、坝基、坝区、泄水设施以及坝体有无滑坡、裂缝、渗流等情况进行巡回监测。
(3)拦水淤积监测
选取流域内有代表性的6座淤地坝,其中骨干坝、中型坝、小型坝各2座,采用常规和GPS两种方法监测,对洪水期间淤地坝的平均淤积厚度及淤积量进行监测。
(1)水文监测
用王茂沟流域水文测站,实时观测把口站次洪水的输沙情况,并根据水位观测几率和采集水样,计算流域次洪水的洪水流量、含沙量、输沙量、侵蚀量等。
(2)降雨量监测
流域内布设自动雨量观测器一台,监测流域汛期降雨量。
(3)气象监测
气象数据采用韭园沟示范区试验场自动气象站监测数据,主要监测项目包括降水、气温、湿度、风速、风向五要素,蒸发采用口径为20cm蒸发皿人工观测。
(1)植被监测
植被监测采取选取典型样区,采用常规方法对流域各地类植被进行监测,每年监测两次。
(2)土壤水分监测
土壤监测主要对坡地、梯田、坝地,采用电子水分仪,取土深0~50cm,时间为每年6月1日~10月1日,每月15日、30(31)日取样两次,监测各地类土壤水分含量。
(3)土壤养分监测
土壤养分监测每年进行一次,并对各地类土壤养分含量进行数据测定。
(1)土地利用监测
利用1/10000地形图实地勾绘,对小板块用皮尺和测绳进行逐块测量,并在室内用求积仪量算。
(2)社会、经济监测
采用调查统计方法对流域人口、劳力、文化程度、农业、牧业、产业结构等进行监测。
(1)基本情况监测
流域内共选取典型农户6户,好、中、差各两户。监测人员直接进行调查方法,对家庭的人口、劳力、文化层次结构等进行监测。
(2)土地结构监测
根据典型农户自报和监测人员直接进行调查两种方式进行。对农户所经营各类土地面积、各类作物种植面积、投劳等情况进行监测。
(3)粮食产量监测
对典型农户各类作物的产量进行监测,根据典型农户自报和监测人员直接进行调查获得。
(4)投入产出监测
对农户所经营各类土地,种子、有机肥、化肥、劳力、农膜等监测,数据根据典型农户自报和监测人员直接进行调查获得。
(5)农户经济结构监测
根据典型农户自报和监测人员直接进行调查获得的数据。对典型农户的种植业、牧业、林果业等收入情况监测。
(6)农户支出监测
根据典型农户自报和监测人员直接进行调查获得。对典型农户的医疗保健、教育、生活、饮食等消费支出进行监测。
选择具有代表性的坝地、坡地、梯田、水地对其投入、产出、净效益等情况进行典型地块监测。数据根据监测人员直接进行调查获得。
小流域坝系建设作为黄土丘陵沟壑区水土保持综合治理中重要的措施,它不仅有效地防止了沟岸扩张、沟底下切和沟头延伸,对整个流域水沙运行起到了拦蓄和调洪的作用。同时在有效地增加基本农田的基础上,可以促进坡地的退耕还林还草,进而实现流域生态环境建设的可持续发展,所以坝系监测作为小流域水土保持或坝系建设综合效益评价的指标体系,不仅要监测其坝系自身建设安全稳定、蓄水拦沙效益和坝地增产效益等,还要包括坝系所处的小流域整体生态建设过程中的土壤性质、植被建设、小气候环境变化和沟道水文泥沙变化等生态指标以及其所引起的小流域社会、经济环境的动态效益。从而将坝系建设作为整个小流域生态建设的切入点和突破口进行系统监测和效益评价,全面、系统地监控和准确、及时掌握坝系建设的运行情况和小流域整体生态建设效果。
通过小流域坝系拦沙蓄水监测方法的研究,认识和掌握淤地坝对流域水沙的拦截、调节和蓄存机理,以及流域水沙在坝系中的演进过程,揭示坝系相对稳定的规律;量化小流域坝系水沙来源指标,为淤地坝工程规划、设计、施工、运行管理、河道水资源的优化配置与合理利用,以及坝地防洪保收技术提供科学的依据。
小流域坝系拦沙蓄水监测主要包括淤地坝拦沙和蓄水监测两个方面。拦沙蓄水监测主要是监测次暴雨后淤地坝的拦沙和蓄水情况。监测指标是淤地坝拦沙量、蓄水量。根据淤地坝是否设有溢洪道分为两种类型:不设溢洪道淤地坝的拦沙量、蓄水量监测和设溢洪道淤地坝的拦沙量、蓄水量监测。
(1)不设溢洪道淤地坝的拦沙量、蓄水量监测
1)监测站点的布设:根据洪水泥沙在淤地坝内的淤积规律(坝前淤积较坝尾淤积薄,坝前蓄水较坝尾蓄水深),从坝前到淤积末端,以控制淤积体平面变化为原则,按相邻间距小于淤积总长度1/6~1/10布设若干个断面,测量断面间的间距。在每一断面布设若干个测点(能控制淤积断面起伏),并在各高程监测站点布设带有水位标尺的水泥桩,测量各监测站点的高程和测点间的水平距离。汛前将带有水位标尺的水泥桩布设在监测站点上,记录下各个监测站点的原始高程数据。
2)监测数据的获取与计算:次暴雨后,我们到选定的监测淤地坝所在地进行实地监测,利用各监测站点水泥桩上的水位标尺的读数,读取次暴雨后监测淤地坝的各个监测站点上的水面高程 Z j ,然后采用平均淤积高程法,计算出所监测淤地坝的总库容。
各断面的平均高程和水面的平均高程用下式计算:
式中 Z ——坝区淤积面的平均高程(m);
Z i ——第 i 断面的平均淤积高程(m);
Z j ——第 i 断面第 j 测点的水面高程(m);
Δ L i ——相邻断面的间距(m);
L =∑Δ L i ——坝前到淤积末端的长度(m);
B j ——同断面相邻测点间的水平距离;
B i =∑ B j ——第 i 断面淤积面的宽度(m)。
由原始库容曲线查得与高程 Z 相应的库容,即为次暴雨后监测淤地坝的总库容 W 总 。
若干天后,监测淤地坝内的蓄水排干,利用各监测站点水泥桩上的水位标尺读数,读取次暴雨后监测淤地坝的各个监测站点上的淤积高程 Z j ,然后采用平均淤积高程法,计算出所监测淤地坝的淤积库容。
各断面的平均高程和水面的平均高程计算公式同式(4-1)、(4-2)。
由原始库容曲线查得与高程 Z 相应的库容,即为次暴雨后监测淤地坝的淤积库容 W 淤 ,即为淤积体积,淤积体积再根据泥沙容重换算成重量,即为监测淤地坝的拦沙量。
然后,由监测淤地坝的总库容减去淤积库容后,剩余部分就是监测淤地坝的蓄水库容 W 蓄 ,即: W 总 - W 淤 = W 蓄 ,计算出的蓄水库容就是监测淤地坝的蓄水量。
(2)设溢洪道淤地坝的拦沙量、蓄水量监测
在淤地坝淤满前,设溢洪道淤地坝的拦沙量、蓄水量监测与不设溢洪道淤地坝的拦沙量、蓄水量监测的方法相同,淤地坝淤满后,在监测淤地坝的溢洪道上布设一个测流断面,当淤地坝溢洪道排水时进行流量、含沙量监测。
1)测流断面的布设:溢洪道纵断面多为矩型,布设测流断面时,在溢洪道便于观测的一侧将直立式水尺垂直安装地基面上;也可用红漆画在观测位置。并在设置水尺的上下游有长约10m以上的平直段,且不受下游回水的影响,床质均一,无影响水流的杂草、淤泥和碎石等杂物。
2)水位观测及沙样提取:在小流域坝系水文泥沙监测中,实时观测洪水期水位的涨落变化过程显得尤为重要。通过观测的水位与流量建立的关系,可以推求洪水流量。次暴雨洪水观测水位时,根据洪水涨落情况增加测次,以能测得完整的水位变化过程为原则。次暴雨洪水的沙样提取次数以能控制住含沙量的变化为原则,在洪水起涨、峰腰、峰顶、落平和水位转折变化点均提取沙样,每次较大洪峰过程一般应不少于7~10次。水沙样的处理经常采用置换法。
3)次暴雨洪水总量及泥沙总量的计算
利用观测到的水位及实验方法率定的水位与流量之间的关系式,计算出瞬时流量。同样将取回的沙样,经过处理,求出单位含沙量。然后利用加权平均流量法计算出次暴雨洪水总量和泥沙总量(详细内容见水文泥沙监测方法)。
淤地坝坝系安全稳定监测,根据其工程组成和生产运行主要包括两个部分,即坝体及泄水建筑物监测和坝系安全运行监测。
(1)坝体及泄水建筑物监测
1)监测指标:由于该流域坝系建设已经较为完善,而且都是20世纪90年代以前建设的坝体和泄水建筑物,所以坝体及泄水建筑物监测的主要内容包括坝体滑坡、位移、冲刷情况、渗流、沉陷和裂缝等涉及坝体和泄水建筑物安全的主要内容。
2)监测方法及数据资料获取:坝体及泄水建筑物安全监测的主要方法采用巡视检查法;对流域卡口站骨干坝要布设监测点进行定位监测。巡视检查直接在现场进行勘测、调查统计和现场拍摄影像进行监测,巡视频次应在每年汛前、汛后按照规定的巡查项目进行一次;定位监测要在坝体左、中和右岸进行埋设监测桩,进行GPS定位,在汛期或年度开展标尺数据记载,巡查其位移、沉陷和滑坡情况。在流域出现特大洪水、感应性地震或持续高水位等比较严重的破坏性自然灾害因素和出现其他危情时,应对规定巡视的项目内容进行逐项连续监测,详细记载时间、部位、险情,并绘制出草图等适时现场情况。
3)监测资料成果分析:对常规无安全性的监测资料进行整理、归档,建立流域坝系监测资料数据库,数据性资料进行整理、输入和图形、表述性处理。对病危性坝体和泄水建筑物要及时对危害情况上报有关上级部门,提供快速决策的依据。
(2)坝系安全运行监测
1)监测指标:主要是对整个流域坝系内正在运行的生产或拦洪工程,逐个进行病坝和险坝的具体坝名、数量、位置,毁坏情况包括水毁、人为毁坏程度、位置、尺寸以及毁坏原因(设计标准、工程质量等)。
2)监测方法:在监测流域按计划每年在汛后和汛期较大暴雨洪水发生后,对全流域分别采取现场勘察测量、调查统计核实坝系安全确定的病坝和险坝的具体坝名、数量、位置,毁坏情况包括水毁、人为毁坏程度、位置、尺寸等,对毁坏性工程现场判断确认毁坏的主要原因(设计标准、工程质量等)。
3)监测资料成果分析:在建立小流域坝系监测资料库的基础上,每年对次暴雨或年度监测资料进行一次整理分析,对工程毁坏性数量和具体尺寸程度等定量指标,要进行数理统计和处理,编制成果图、表等,计算其小流域坝系的安全比。对毁坏性工程位置、名称、人为或水毁等定性指标,要进行质量评定,确定流域坝系安全程度。通过多年坝系动态监测,可以掌握和了解流域坝系动态运行情况,确保流域坝系监测成果的科学性和系统性,为流域坝系工程建设提供技术支持和科学决策依据。
小流域坝系水文泥沙常规监测,监测的主要内容是降雨量、流量和含沙量,以求出小流域坝系次暴雨洪水总量和来沙总量。采用的方法是建立雨量站进行降雨观测和在小流域坝系出口修建卡口站进行水文泥沙项目的观测。卡口站可监测小流域坝系每次暴雨洪水上游来水来沙和出口断面水沙情况,监测典型淤地坝坝系的水沙演进过程。通过卡口站对小流域坝系水沙的监测,可计算出小流域坝系的来水来沙量,为小流域坝系效益分析评价提供科学依据。水文泥沙观测的主要项目有降雨量、水位、流量、含沙量。
监测指标的确定为小流域坝系水文泥沙监测的指标有降雨量(mm)、流量(m 3 /s)、含沙量(kg/m 3 )。
(1)降雨量监测方法
1)观测场地选择与雨量点密度:收集降雨受周围环境影响很大,因而观测场地应避开强风区,建在四周开阔、平坦,并保证降水在倾斜下降时,四周物体不致影响降水落入雨量器内。四周障碍物与仪器的距离不得小于障碍物顶部与仪器口高差的2倍。
2)小流域雨量点密度:小流域雨量点的多少受流域面积、形状和地形变化大小的制约,也随降雨观测服务的目的而变。以小流域坝系监测水文泥沙为服务目的的降雨观测,一般采用1km 2 面积,至少应有1个雨量点。
3)仪器的安装:观测雨量要用雨量器、自记雨量计等。雨量器安装在有人驻守的雨量点上。雨量器安装高度为:器口至地面高0.7m,且保持器口水平,三角架深入地面以下要牢固。自记雨量计常用的为虹吸式的,它可记录降水过程及雨量变化,需要观测人员经常检查、换纸、加墨水等工作,因而常安装在径流站。它的安装高度为0.7m或1.2m。为保持稳定、牢靠,应在仪器底部埋设基桩或砖块,并将三条线拉紧埋实,注意保证器口水平。
4)降雨量观测及其结果计算:降雨量观测多采用雨量器用2时段观测,即每天早8时和晚20时观测,这是为满足研究的需要。若遇少雨或无雨天也可在早8时观测一次。自记雨量计观测在每日8时整进行,当降雨量多或暴雨时要经常巡视,以排除故障,防止漏记降雨过程。
小流域雨量观测结果的资料整理与计算是建立在所设点的雨量资料校正准确的基础上。整理的内容主要是求流域的平均次降雨量。流域平均雨量推求的方法一般有两种:算术平均值法、面积加权法。
算术平均值法:在小流域布置雨量点密度足够,且地形单一、分布均匀时,可用算术平均值法推求小流域平均雨量。公式为:
式中 ——流域平均次雨量(mm);
P i ——流域所设点雨量的校正次雨量(mm);
n ——流域布设的点雨量数。
面积加权法:若流域雨量点较稀、分布不均,或地形起伏大,则用面积加权法计算平均雨量。各雨量点面积的确定多帮泰森多边形法。此法是将布设的雨量点依据紧邻关系连成多个三角形,从周边雨量点起作三角形每个边的中垂线,再延长连接流域中部雨量点周围各线的中垂线,形成围绕雨量点的多边形(含流域边界线),这些多边形所包围的面积即为被包围雨量点的面积。用各雨量点的面积求出占总面积的份数,即面积权数 f i ,用下式计算流域平均雨量:
式中 ——流域平均次雨量(mm);
P i ——各点雨量(mm);
f i ——各点雨量的面积权数。
f i = F i / F ( F i 为该雨量点覆盖面积, F 为流域总面积)。
(2)小流域坝系来水来沙量的监测方法
小流域坝系次暴雨洪水的来水来沙变化快、变化大,其流量、含沙量的监测一般采用测流建筑物量测和断面量测。在水土保持水文泥沙的常规监测中,经常使用的量水建筑物有巴塞尔量水槽、薄壁量水堰(按出口形状分为三角形、矩形、梯形等)、三角形量水槽等。本次课题研究我们以王茂沟小流域坝系为示范基地。王茂沟小流域坝系水文泥沙常规监测方法,采用在王茂沟小流域坝系出口处建立卡口站,进行小流域坝系水文泥沙监测,具体分为以下几个步骤:
1)站址选择及量水建筑物:选址的好坏,直接影响到测流的精度。我们选址时是按以下标准进行的,即水流流动顺畅,无弯道和宽窄变化的沟道,且床质均一,无突兀大石头和凹陷坑穴,边岸杂草不影响水流。王茂沟小流域坝系水沙监测布设的量水建筑物为三角形量水堰和矩形宽顶堰两种,测流断面设在王茂沟小流域出口处1 # 骨干坝的溢洪道上。当洪水流量小于0.43m 3 /s时,在三角形量水堰量水建筑物上进行观测;当洪水大于0.43m 3 /s时,在矩形宽顶堰量水建筑物上进行观测,量水堰底宽3m。
2)水位泥沙观测方法:王茂沟小流域水沙监测过程中水位观测的频次,根据当时次暴雨洪水水位变化的急缓情况,以能测出完整的水位变化过程,满足流量计算为原则。一般采用的标准是“快一分、慢一厘”,即洪水水位涨或落变化快时采用一分钟记一次水位,洪水水位涨或落变化慢时采用水位变化1厘米记一次水位。次暴雨洪水的水样提取频次以能控制住含沙量的变化为原则,在洪水起涨、峰腰、峰顶、落平和水位转折变化点均提取水样,每次较大洪峰过程一般应不少于7~10次,采样一般用小桶或者储存在其他采集器皿。
水样的处理经常采用的方法是置换法。主要工作内容有:量沙样容积;测定比重瓶盛满浑水的总重量;测定浑水的温度;计算泥沙重量。比重瓶一般采用容积为100cm 3 的,注入水样时,不能太急,应使浑水沿壁徐徐流下,当浑水注入后,可用手指轻击比重瓶的四周,以助气泡外逸;然后再注入少量清水,使水面到达一定的刻度时,加上瓶塞,用干毛巾擦干瓶外水分。称瓶加浑水重,并用水温度计迅速测定其温度。泥沙的重量用下式计算:
式中 W s ——泥沙重量(g);
W ws ——瓶加浑水重(g);
W w ——同温度下瓶加清水重(克),由曲线上查得。
3)次暴雨洪水总量及泥沙总量的计算:次暴雨洪水总量计算首先根据观测到的水位及水位与流量的关系式求出瞬时流量。王茂沟小流域坝系的水沙监测选用的量水建筑物有两种:三角堰和矩形宽顶堰。根据水位计算流量的公式分别是:
三角堰流量公式为:
式中 Q ——瞬时流量(m 3 /s);
H ——水位高度(m);
u ——常数,为0.6;
θ ——三角堰顶角,为120 °;
g ——常数,为9.8。
矩形宽顶堰流量公式为:
式中 Q ——瞬时流量(m 3 /s);
H ——水位高度(m);
m ——常数,0.785;
b ——底宽(m);
g ——常数,为9.8。
根据上式求出的瞬时流量,用面积包围法计算出次暴雨洪水总量。即用相邻两个瞬时流量相加,除以2得出时段平均流量,再乘以相邻两瞬时流量间的历时,得出该时段的洪水流量。以此类推,求出整个次暴雨洪水过程相邻两个瞬时流量在某一时段内的洪水流量,然后将这些求得的时段洪水流量全部累加起来,所得的结果就是此次的次暴雨洪水的总量。
用流量加权法计算出次暴雨洪水泥沙总量。首先根据以上观测计算出的泥沙重量,求出单位含沙量;一次洪水过程一般处理水样7~10个,也就是说能求出7~10年单位含沙量,然后用直线插补法,计算出每一个水位流量对应的单位含沙量。整个次暴雨洪水过程所有时段内的洪水流量乘以单位含沙量得出的每个时段内的产沙量,然后将这些求得的时段产沙量全部累加起来,所得的结果就是此次的次暴雨洪水的泥沙总量。
生态监测是指坝系建设所引起的小流域生态环境土壤、植被状况和小气候变化的动态监测。土壤环境的监测主要对不同地类影响土壤性质的水分含量、孔隙度、干容重等物理指标和不同地类土壤有机质(%)、碱解氮(mg/kg)、速效磷(mg/kg)、速效钾(mg/kg)、pH值等化学指标的监测;植被状况监测是对流域内主要分布,并能反映小流域林草植被建设的动态变化和建设效果的人工乔木林、经济林、灌木林、牧草和草灌混交的5种地类测定其郁闭度和盖度;小气候监测也是对坝系建设流域的主要气象要素降水、气温、空气湿度、蒸发和风力进行动态监测,降水主要是监测其汛期(植物有效生育期)逐日降水量和发生暴雨时的时段降水量,气温和空气湿度主要监测其汛期(植物有效生育期)逐日温度、湿度,监测时段最高、最低气温和湿度等,蒸发主要监测其汛期(植物有效生育期)逐日蒸发量,风力主要是对风速和风向进行定时监测。
(1)生态、气象要素监测方法确定
1)土壤物理性质监测方法:在选定土壤监测样区内,将地表上部2~3cm带有植物残体的表层土刮去,采用环刀法向下取10~20cm的土样在实验室测定其土壤孔隙度、干容重等,土壤含水率采用时域反射电子水分仪,在汛期(植物有效生育期)内每月15、30日进行测定,根据林草根系分布确定水分监测深度即农牧地类型监测20、30、50、100cm四层,林地类型监测20、30、50、100、150、200cm六层水分含量。
2)土壤化学性质监测方法:在选定的每个监测样区内,将表层2~3cm的土壤及其植物残体刮去,取0~25cm深的土壤,在样地上采集点用梅花分布采集法取各点土样混合约1kg,在室内风干后再测定。有机质采用重铬酸钾容重法,碱解氮采用扩散法,速效磷采用01SEN法,速效钾采用原子吸收法,pH值采用比色法。
3)植被状况监测方法:选择流域内分布较广、代表性较强的人工草地、灌木林、乔木林、草灌混交林和经济林五大类型,各种类型选择3个标准样地,考虑到植物生长的开始和结束,于每年的6月和9月各监测1次。林木郁闭度采用树冠投影法,灌木盖度采用线段法,草地盖度采用针刺法。
4)小气候环境变化监测方法:考虑到监测流域的建设条件,在相同类型区建立小气候地面自动监测站,监测站为一周监测一次,定时回传降水、温度、空气湿度和风力等实时动态变化数据,然后对回传的数据进行所需成果数据的编报。软件支持系统包括实时监测系统和地面气象编报两个系统。
(2)生态效益监测样区布设
根据监测所确定的监测项目土壤性质、植被状况等,监测的样地类型有坡地、梯田和坝地,土地利用方式有乔木林、灌木林和农地等10种(详见表4-1)。在小流域中选择有代表性的人工乔木林地、灌木林地、草地、果园分别布设监测小区,监测其土壤理化性质、郁闭度和盖度。其中林木、牧草样方的位置应在坡面的上、中、下部位分别选定,乔木林样方面积20m×20m,灌木林5m×5m,草地2m×2m。
表4-1 土壤理化性质观测小区基本情况
(3)监测数据成果分析
1)土壤性质监测成果分析:土壤物理性质根据监测数据进行相同类型地类的动态比较分析,不同地类进行数据的类比分析。化学养分分析,根据监测要素中的各个参评项目将土壤养分评价划分为5个等级,每个等级所对应的得分按等比数列划分指数。把参评的各个数据按照参评项目划分的级别分别列出,然后用参评因素中各个参评项目的得分与各参评因素权重的乘积之和,除以参评项目的最高得分与各参评因素权重的乘积之和,得到该土壤某参评因素的得分值,其值越大土壤肥力就越高,反之则越小。土壤养分评价项目与评价指数见表4-2。
表4-2 土壤化学养分参评项目与评价指数
2)林草植被监测数据分析:林草植被覆盖情况变化,按照覆盖度等级划分的方法,即Braun-Blanguet的覆盖度等级划分方法,将植被覆盖度划分为5级,详见表4-3。其中同一类型植被覆盖度变化评价,通过对同一类型植被覆盖度在6月、9月两次监测,计算植物年内、年际增加幅度的大小次序,作为林草植被建设好坏指标评价次序;不同类型植被覆盖度变化评价,分析其在年内、年际增加幅度等级次序,作为指标评价依据。
表4-3 植被覆盖度等级
3)小气候环境变化监测数据分析:温度、蒸发监测数据分析。采用气象要素数据分析技术规范中的差值法分析,按照监测流域各年气候要素月、季的平均值,计算分析监测流域小气候环境的变化;应用气候学理论分析监测区小气候变化的原因,进而评价坝系和相关联的生态建设对小气候的影响程度。也可应用对比法分析小气候的变化,用绥德气象站的资料与监测流域的监测数据进行各要素的变化比较。由于监测流域气候资料系列短,故需将短系列的资料延长,处理方法采用称为膨化处理,经膨化处理的数据资料均需进行距平处理,其计算公式为:
标准化处理公式:
式中 X j ——处理后的新序列;
K ——汛期(6~9月)4个月;
X j ( K )——逐年各月气候要素值组成的原序列;
——多年各月气候要素的平均值;
S n - i ( K )——月气候要素的标准差。
对比分析时,采用回归分析法求各参照点(即自变量)与基本点(即应变量)对监测年以前绥德县气象站数据进行回归计算新序列的最佳方程。然后将本次监测数据参照点的新序列值代入(4-8)式或(4-9)式,求出监测流域的气候要素的值,将基本点气候要素天然值减去实测值,其差值即为坝系建设流域对当地小气候影响所引起的变化部分,最终评价其变化原因。
降水量、湿度、风速等监测数据分析。该类监测数据采用比值法分析,假定大范围气候变化所引起的气候要素在项目区建设前后的比值不变,就可以利用参照点前后气候要素的平均比值,推算出监测年份相应流域对气候要素影响的增、减量,然后分析评价其变化原因。
黄土丘陵沟壑区坝系建设做水土保持生态环境建设的基础。生态环境且是该区人们生存的基本条件,所以从在该区建设沟道坝系,已决定了其具有广泛的社会属性,而坝系建设监测流域的经济社会效益则是这种属性的反映。沟道坝系建设的社会效益的含义即坝系建设后,对社会环境系统和所在流域经济系统的影响及产生的宏观经济变化和社会效应。也就是坝系建设离不开为实现社会发展目标和全面实现国民经济发展目标所作的贡献与影响程度。所以一个完整坝系建设小流域,其对流域不同社会群体和个体经济利益影响、社会进步要素提高的总和即为经济效益和社会效益。同时效益监测与评价密不可分,后者必须以前者为基础。小流域坝系建设区域主要是广大农村,农村应是实施监测的主体,而广大农户是项目建设的直接参与者,同时又是受益人——基本主体单元。因此,项目社会效益监测应以农户监测为基础,并附以宏观的社会经济调查,其主要反映农民经济收入的增加和项目实施后流域水土保持生态建设的投入产出两个方面。而流域社会效益监测和经济效益监测的侧重点有所不同,社会效益监测的内容更为广泛一些,在经济效益监测的基础上,即农业增产、农民增收的基础上,促进小流域农村教育、卫生、交通和文化科技水平的进步程度高低来反映。
(1)经济效益监测
1)监测布设:以坝系建设的监测流域为单元,选择典型农户进行定点调查为主,并辅助布设典型地块监测。典型农户应在监测流域内按农户经济水平分好、中、差3个档次选取、布设典型农户6户,各个层次按2户进行经济效益监测。好、中、差的标准按年人均经济纯收入、年人均占有粮食、恩格尔系数等3个指标确定,各指标选取标准见表4-4。为了获取典型农户在项目实施后获得的增产效益,在典型农户经营坝地、梯田等各类土地的基础上,选择典型地块进行监测。典型地块为坝地、梯田、农坡地,监测中将3种农地进行对比,同时对乔、灌、草、经济林及养殖户的羊、猪、兔等进行专项监测,共布设监测点6处。
表4-4 典型农户选取指标标准
注:恩格尔系数是指食品指出占生活消费总支出的比例。
2)监测方法:典型农户监测,对典型农户和监测人员都要事先进行技术培训,使之熟练地掌握调查内容与调查方法,根据年内家庭经营、收支等情况进行连续监测;典型地块监测,对坝地、梯田、农坡地以及果园、林、草地典型地块的粮食和果、林、草产品(果品、枝柴、草料)产量均应单打单收,求得其单位面积产量,并了解其增产的具体原因。
3)监测指标:典型农户监测指标。Ⅰ家庭人口、劳力、文化程度,承包经营的农地(梯田、坝地、水地、河滩地、坡地)、林草地,农业、非农业用劳(含妇女用劳,农业、非农业剩余劳动力),粮食、林、牧、果和运输、加工、编制、劳务输出等收入,人均经济纯收入、上缴国家税收、集体提留,家庭生活消费、生产性资产、生活性资产。Ⅱ粮食作物、经济作物种植面积及其产量。Ⅲ农作物生产投入及产出情况。Ⅳ果业生产情况。Ⅴ畜牧业生产及投入产出情况。Ⅵ农业机械、役畜及产品情况。Ⅶ住房、电视、电冰箱、摩托车、自行车等情况。Ⅷ生活、教育、医疗卫生等消费。
典型地块监测指标。Ⅰ典型地块农作物生产情况。Ⅱ典型地块农作物投入产出情况。Ⅲ典型地块人工乔木林、灌木林监测:地类、面积、造林时间、树种、整地方式、造林方式、成活率、地径、树高、树冠、枯落物厚度、郁闭度。Ⅳ典型地块人工种草监测:地类、面积、种植时间、草种、播种量、密度、收割次数、投入、总产量、总产值。Ⅴ典型地块果园和经济林监测:地类、面积、品种、种植时间、投入、成活率、密度、总产量、总产值。
4)监测数据成果分析:在各类监测数据归类整理总结的基础上,着重分析和反映监测流域在监测年内投入、产出动态净现值和内部回收率,监测的投入成本主要考虑措施成本和运行成本,措施成本主要是监测流域各类措施的投资数值,运行成本是典型农户和典型地块监测数据与市场调查结果综合分析确定。产出主要包括农业效益(粮食、经济作物、瓜果蔬菜等主副产品产量和产值)、林果业效益(木材、枝条、果品和种籽等)、种草效益(饲草和草籽价值)。最后根据以上动态数据和利率计算分析净现值(NPV)和内部收益率。
(2)社会效益监测
1)监测布设:社会效益监测可分为农户监测和监测流域宏观社会经济状况监测。农户监测布设又分为典型农户监测和样本农户监测。
典型农户监测。典型农户监测的目的是选取典型,以便深入、细致地了解监测点上的生产经营运行情况。在监测流域按经济水平分好、中、差3个档次选取、布设要监测的农户。好、中、差的标准可按年人均经济纯收入、年人均占有粮食、恩格尔系数等3个指标确定。为了解典型农户在实施各项水土保持措施后所获得的增产效益,在典型农户经营的土地上布设典型地块进行监测。典型地块中对水地、坝地、梯田、坡地进行对比监测,对乔、灌、草、经济林及养殖户中羊、猪、兔等进行专项或典型地块监测,示范区共布设了6个典型农户监测点。
样本农户监测。为了解监测流域面上生产经营运行情况,保证监测成果具有一定的代表性,在监测流域内按随机抽样的方法,抽3家农户进行样本监测。
项目区宏观社会经济状况监测。根据监测流域的实际情况以村为单位进行,监测的内容有:土地、人口、劳力;土地利用结构;耕地面积变化;农村产业结构变化;社会进步情况。
2)监测方法:典型农户定点监测。对被选取的典型农户随着监测年度的开展进行长期、连续、定点监测。监测布设完毕后,对所选的典型农户进行监测要求培训并定期指导,按设计的固定格式连续进行记录,定期汇总记录成果。
一次性调查。将监测流域样本农户调查与社会经济调查结合起来,在项目实施中每年调查一次,调查的内容需事先设计好,并采取问卷笔录式逐户调查。
3)监测数据成果分析:在典型农户、样本农户和监测流域宏观社会经济状况监测布设和数据监测的基础上,主要分析监测流域农民生活水平的提高程度,即人均粮食增加幅度、人均现金收入增加程度;土地利用结构的变化,即基本农田、压缩坡耕地和林草地的增加程度;农村生产结构变化,即在总产值中农、林、牧、副和劳力输出等各业产值变化情况及因素分析;劳动生产率的提高程度以及促进农村社会化进步的文化、教育、卫生、科技和交通等事业的投入程度等。
GPS技术在水土保持监测中,宏观方面可建立GPS控制网,在控制网的基础上,进行控制点测量,为卫星遥感技术的定向提供加密点,也可以于宏观区域和重点区域水土保持信息的采集、提取;在微观方面,可利用GPS技术监测水土保持措施图斑内容、几何量算、措施质量等,如梯田、淤地坝、乔木林、灌木林、人工种草等措施的面积监测。同时可以进行水土流失实时监测,如沟头前进、沟底下切、沟岸扩张的速度,甚至可以监测典型样点水土流失量等。本次坝系监测主要应用于流域高程网点控制测量和淤地坝面积以及淤地坝淤积监测。
(1)测量仪器
本次投入9600型GPS仪器一套(三台套),其静态相对定位精度为(5mm+1ppm),动态水平精度优于1m,高程精度为±(10mm+2ppm)。
采用软件是SOUTH南方测绘仪器公司的9600型GPS测量系统GPSADJ专用基线处理与CADR14平差软件。
普通水准仪、DJ6-1光学经纬仪各一台。
(2)利用资料及座标系
1)王茂沟流域1∶10000的地形图。
2)王茂沟流域内国家三角点东山圪塔Ⅲ级控制点、流域附近魏家焉国家三角点魏家焉Ⅱ级控制点。
3)坐标系统采用WGS—84坐标系进行无约束平差,然后进行1954北京坐标系,高程采用1956黄海高程。
(3)控制网的建立
1)选点埋石:GPS控制网的平面观测于2004年3月起,我们对王茂庄沟及其附近的国家三角点进行勘查,共勘测5个点,最后确定流域内东山圪塔一个国家三级控制点,流域外魏家焉一个国家控制二级制点,两控制点相距6.7 km,两标石完好无损。为了便于王茂沟流域淤地坝监测,我们在王茂沟黑兔焉、1 # 、2 # 坝增设了三个控制点,高舍沟增设一个点。
根据监测设计要求,流域内共埋设D级GPS点标石3座,流域附近埋设1座。标石规格为20×20×100cm。GPS点位选在土质坚实,交通便利,方便使用,便于长期保存的地方。
2)GPS野外观测测量:GPS控制网的平面观测于2004年4月10日至12日完成。GPS淤地坝现状碎部测量于2004年4月23日完成。
GPS技术规定:为了能够使观测数据可靠精确,以免造成人力、物力、资金的不必要浪费,依照国家GPS的有关规范进行设计实施。(详见表4-5)。
表4-5 基本技术规定
观测方案按设计书的要求进行,利用魏家焉和王茂沟东山圪塔国家三角点与4个GPS点联合构网进行观测,网中最长边6.2km,最短边0.7km,平均边长3.3km。作业基本技术指标见表4-6。
表4-6 作业模式及基本技术指标
(4)基线解算
1)基线解算采用南方测绘软件GPSADJ4.0基线处理软件。对野外采集的原始数据*. STH观测文件转换为标准的RINEX观测文件,此文件可在不同的路径下任意选择。读入观测文件后进行向量解算及网平差。
2)基线解算时参考1992年国家GPS测量规范,按以下要求进行解算:基线边长相对中误差要求如表4-7所示。
星历采用广播星历; 采用“码和相位”解算;
电离层采用计算模型; 对流层采用标准模型;
卫星高度角采用15 °; 采样频率5s。
3)闭合环检验:所有基线组成闭合环,闭合环总数4个,同步环总数4个,闭合环最大节点数3个。各坐标分量闭合差的限差按下式规定:
其中: a =10, b =5。闭合环相对闭合差0-1ppm3个,1-2ppm1个,优于国家规范的等级限差(表4-8)。
表4-8 同步闭合环表
4)平差计算及精度:GPS基线向量网首先进行了WGS-84坐标系三维自由网平差。在平差计算时,利用魏家焉( W 001)和王茂沟东山圪塔(D101)国家1954北京坐标系3度分带三角点为起算点进行二维约束平差计算。中央子午线经度采用Ⅲ度。计算结果见环闭合差报告(表4-9、表4-10):
WGS84-坐标系下经自由网平差的平差结果;
三维自由网平差;
单位权中误差:0.015852(m)。
表4-9 基线及其改正表
表4-10 平差后WGS84坐标和点位精度表
当前坐标系统:国家1954北京坐标系(表4-11、表4-12):
椭球长半径:6378245.000000 椭球扁率:1/298.30
M0:1.000000 H=:0.000(投影高)
B0:0.000000 L0=:0.000000(中央子午线)111.000000
N0:0.000000(北向加) E0=:500000.000(东向加)
采用网配合法进行转换;单位权中误差0.000844(m)
表4-11 54北京坐标系点位精度
表4-12 平差后坐标和点位精度表
参数拟合高程109.276507(表4-12~表4-14)。
拟合公式:
内符合精度中误差±206.972(mm),表4-14。
表4-13 拟合后高程残差
表4-14 54北京坐标系拟合高程表
表4-15 坐标及拟合高程
从以上基线解算精度及网平差精度来看:GPS平面测量最弱边相对中误差为0.84ppm,平差计算精度远远优于设计书和相应国家GPS测量规范的要求。说明本次GPS网布网合理,观测安排科学,数据处理和平差严密。
至此通过我们实地选点埋石,野外观测,内业基线解算,在王茂沟小流域黑兔焉、王茂沟1 # 、2 # 坝建立三角控制网。见GPS测量控制网点布设图4-1。
图4-1 王茂沟流域GPS测量控制网图
(5)控制测量结果分析
可以看出,WGS84-坐标系下三维自由网平差,基线解算最差边相对误差1∶745851,精度达1.3×10 -6 ,国家北京54坐标系,采用网配合法进行转换,基线解算最差边相对误差1∶369511,精度达2.7×10 -6 ,基线解算完全符合《GPS测量基本技术要求规范》。高程精度通过已知点正常高与大地高拟合,内符合精度中误差±206.972mm,拟合精度为±146.351mm。为了了解GPS所测同一点高差对淤积厚度的影响,我们对黑兔焉(WW01、WW03点)控制点进行不同时期两次测量,从测量结果看,WGS84-坐标系下两次测量高程差4mm,国家北京54坐标系下,两次测量高程差3mm,完全可以满足淤地坝淤积监测要求。
(6)小流域GPS沟道工程动态监测
采用两台GPS接收机同时作业观测,其中一台作为基准站,一台作流动站,流动站与基准站相距小于15km。流动站在起始点静止初始化10分钟,然后从起点开始,流动在运动过程中,按预先设定的时间间隔自动观测,自动记录数据。在移动过程中保持同步观测4颗以上的卫星,并且保持连续跟踪,在待测点观测一个以上的历元(默认值5秒),实现了“STOP AND GO”的测量过程,定位精度(相对基准点)可达1cm。观测工作结束后,将储存在采集器的数据文件传输到计算机中进行后处理。处理后直接输出坐标成果并显示所有的航迹线。数据经过处理可方便的进入CAD进行图形编辑,数据成果可导入MapInfor等GIS系统。主要应用于勘界测量和面积测量。
为了摸清该流域目前的水土流失情况和淤地坝的淤积情况等,本次对流域22座淤地坝的布局、数量、淤地面积、淤积分布、淤积量、坝地利用以及各淤地坝之间的关系等进行监测。结合流域坝系特征,将流域分为上、中、下三个断面,在各断面选取有代表性的王茂沟1 # 、2 # 骨干坝,黄柏沟、马地咀中型淤地坝和何家峁、王踏沟小型淤地坝共6座进行监测。
1)GPS差分动态监测:动态相对定位模式主要应用于地籍勘测和面积测量,我们利用GPS差分动态相对定位模式的这一特点,通过7天时间对王茂沟小流域的22座淤地坝汛前数量、淤积面积进行监测。监测结果见表4-16、表4-17。
表4-16 小流域沟道工程总体建设统计表
表4-17 小流域坝系GPS差分监测结果表
续表
2)高程网点布设:根据洪水泥沙在淤地坝内的淤积规律(坝前淤积较坝尾淤积薄),对流域内王茂沟1 # 、2 # 骨干坝,黄柏沟、马地咀中型淤地坝和何家峁、王塔沟小型淤地坝共6座淤地坝进行GPS汛前、次暴雨(汛期)结束后的高程监测网点数据收集。我们将各监测淤地坝分为上中下三个断面,在各断面布设高程控制监测网点,并在各高程监测点布设带有水位标尺的水泥桩(便于验证GPS监测数据),然后利用流域已建立的三角控制网点,测量上述6座监测坝各个淤积断面的控制网点高程。黄柏沟淤地坝高程监测网点布设见图4-2。
图4-2 黄柏沟淤地坝高程网点布设示意图
利用GPS对淤地坝汛前以及次暴雨(汛期)结束后的面积监测和高程测量,计算出汛前以及次暴雨(汛期)结束后的平均淤地面积和平均淤积厚度,然后计算出次暴雨(汛期)结束后的坝地淤积量。
3)测量结果分析:GPS差分动态相对定位模式,定位精度(相对基准点)可达1cm。根据我们对100m 2 的地块进行实际量算与GPS差分动态相对定位模式多次测量结果标明,GPS差分动态测量面积误差在±0.45%左右,完全可以满足淤地坝的面积测量。在实际测量过程中,其几何精度因子PDOP一般都在2与6之间,定位精度较高。但在部分沟道测量中,其定位精度有的超过分米级,特别是后差分动态测量时,在沟道狭窄的地方,几何精度因子PDOP有的超过了7,甚至达到10,定位精度较差,当然这还不能满足水土流失动态监测对精度方面的要求。造成这种精度差异除了仪器本身的原因外,还有一个重要的原因,即与黄土丘陵沟壑区地面开阔度对GPS接收机接收卫星信号强弱有关。但通过用其他测量仪器(全站仪、经纬仪、水准仪)进行修测和补测。
数字地形模型是地形表面形态属性信息的数字表达,是带有空间位置特征和地形属性特征的数字描述。数字地形模型中地形属性为高程时称为数字高程模型(Digital Elevation Model,简称DEM)。高程是地理空间中的第三维坐标。由于传统的地理信息系统的数据结构都是二维的,数字高程模型的建立是一个必要的补充。DEM通常用地表规则网格单元构成的高程矩阵表示,广义的DEM还包括等高线、三角网等所有表达地面高程的数字表示。在地理信息系统中,DEM是建立DTM的基础数据,其他的地形要素可由DEM直接或间接导出,称为“派生数据”,如坡度、坡向、粗糙度等,并进行通视分析、流域结构生成等应用分析。因此,DEM将会在水土保持领域中得到广泛使用。
本项目由于涉及到对地形地貌、淤地坝淤积、土壤侵蚀等研究,所以完成了王茂沟小流域乃至整个韭园沟流域的DEM。同时本项目研究中,由于使用高分辨率QuickBird卫星数据,需要与DEM数据叠加分析,定量研究淤地坝坝系监测中应用,将该研究区域的DEM分辨率定为1m。并把本次应用的所有数据都统一在高斯-克吕格投影,1954年北京坐标系。
DEM的数据采集方法通常有以下几种:现有地图数字化、地面测量、空间传感器、数字摄影测量方法。本项目研究采用的是现有地图数字化。其工作流程如图4-3所示。
图4-3 DEM建立的技术路线
(1)小流域数字高程模型的建立
1)数据资料准备:
a .数字化底图:1∶1万地形图,共9幅,等高距为10m。
b .图幅控制点坐标:图幅控制点坐标是用来进行图幅定向的,它能够确定地图的地理位置和比例大小。
c .确定地图的分层与分幅:GIS是以图层的方式管理地图的,将点、线、面等地理实体按其性质的不同分别归入不同的图层进行分层管理是GIS管理空间数据的基本形式,本研究主要分为:等高线层(terlk)、线状水系层(wtlpt)、面状水系层(wtlnt)。
d .设计代码:代码的设计非常重要,它是计算机存贮、检索、识别的基础,使之能够满足各种应用分析需求。
e .在计算机上建立自己的工作目录,将地图影像放入,文件格式为TIF。
2)DEM制作步骤与方法:
a .地图预处理:
变形纠正 扫描过程中由于图形倾斜,造成扫描后的地图产生变形,因此在数字化前,需对地图进行纠正。如果存在几何变形,可利用ERDAS的几何校正模块进行纠正;如果是扫描过程造成的图像倾斜,则可在PHOTODHOP下进行纠正。首先,利用[标尺]沿着图幅边缘画一条横线,然后选择[图像]菜单中的[旋转画布]下的[任意角度],在弹出的对话框中会自动计算出旋转角度,点击[好]即可。
二值化 不少数字化软件接受的数据为二值化的数据,因为二值化后的数据量小,很大程度地提高了数字化的速度。二值化也是在PHOTODHOP下处理。选择[图像]菜单中[调整]下的[阈值],移动小三角直到图像满意为止。
处理完将图像保存为TIF文件格式。
b .数字化采集:本次数字化采集应用GEOSCAN软件。
调入底图 打开GEOSCAN,在[调图]中选择[调入栅格图像],在弹出的对话框的[文件类型]中选择TIF,找到存放地图的位置,将图打开。
地图定向 在[地图]中选择[图形定向],在弹出的对话框中选择[齐次方程定向(至少四点)],点击[OK]。将鼠标移到左上方点在角点上,在弹出的放大图上精确定位,并输入坐标,点击[接受量测],其他控制点按顺时针方向以此类推,并回到第一个点,把第一个点再做一遍,选择[结束量测]。
创建图层 在[设置]中选择[图层控制],在这里分别创建所需图层,并以不同颜色区分开来,先设置等高线层terlk为当前图层。
数字化跟踪 用GEOSCAN的数字化工具,可以选择自动跟踪和手动跟踪。
属性赋值 点击要赋值的对象,在弹出对话框中输入高程值及代码。
成果输出 数字化完成之后首先要存盘,选择[调图]中的[保存矢量图形]。然后将成果输出,选择[调图]中的[输出外部格式,并选择[输出到AutoCAD],在弹出的对话框中选择默认值,此时,在工作目录中就会多一个与地图名称相同的dxf文件。
c .数据后处理:数字化完的数据都不可避免存在着错误和误差,属性数据在输入时,也难免会存在错误,因此对图形数据和属性数据进行检查、编辑和处理,是保证数据正确可用的必要条件。本次研究的数据处理是在ArcInfo下进行的。
首先将dxf文件转入到ArcInfo中,形成coverage文件,然后在对各层数据创建拓扑关系,进而消除其图层内的悬挂与相交现象,检查其属性正确与否,修改完所有错误之后,即可存盘。由于对点线作过移动、删减与增添,原先建立的拓扑关系遭到破坏,必须重新建立拓扑关系。
d . DEM的建立:经过数字化采集与编辑处理,为数字高程模型的建立做好了准备,即可创建不规则三角网TIN,再将TIN内插成Grid,即可生成DEM。
(2)利用数字高程模型提取流域特征数字信息
由于数字高程模型DEM可以用于描述地面起伏状况,进而用于提取各种地形参数,其中有:坡度、坡向、坡长、等高线、沟道长、沟壑密度等。本项目应用建立的DEM,利用Arcview及Arcgis系统下完成了各种流域特征数据的提取。
1)坡向的提取:坡向定义为坡面法线在水平面上的投影与正北方向的夹角。在ArcGIS和ArcView两个软件中Aspect表示每个栅格与它相邻的栅格之间沿坡面向下最陡的方向。在输出的坡向数据中,坡向值有如下规定:坡向以度为单位,按逆时针方向从0°(正北方向)到360°(绕圆一圈后的正北方向)来度量。坡向图中每个栅格单元的值表明此栅格单元所在坡的朝向。水平的坡没有方向,被赋值为-1。
a .利用ArcView来制作坡向图:坡向可在数字高程模型DEM或TIN数据的基础上提取,在DEM基础上提取坡向的步骤如下:
ⅰ 在视图目录表中添加DEM并激活它。
ⅱ 从Surface菜单中选择Derive Aspect命令。
ⅲ 显示并激活生成的坡向主题Aspect of DEM(如图4-4)。
图4-4 提取坡向
在DEM或TIN的面主题中坡度为0°(平地)的栅格在输出的坡向主题中被赋值为-1,如果围绕中心栅格的任何相邻栅格是No Data数据,它们将被赋予中心栅格的值,然后计算坡向。在坡向主题的图例中表示了八种主要方向,例如:东[67.5-112.5 °],东南[112.5-157.5]。
b .利用ArcGIS软件来生成坡向图:在水保工作中,坡向有时比高程更为有用,尤其是在造林工作中,我们可以在选择树种的时候,阴坡由于其水分条件较好,适宜栽植的树种就比较多,而阳坡则适宜栽植的树种比较少。对于坝系建设来说,阳坡取土场要比阴坡取场早解冻10天左右,这给淤地坝施工带来很大的方便。
对于Grid或TIN来说,其坡向可以动态地进行快速计算,可以生成坡向度。具体步骤如下:
ⅰ 点击3D Analyst菜单,将鼠标指向Surface Analysis,选择Aspect命令(图4-5)。
ⅱ 单击Input surface下拉列表的箭头,选择用来计算生成坡向栅格图。这里选择wmgdem。
图4-5 Aspect命令菜单
图4-6 指定名称图
ⅲ 指定输出栅格图像的文件名。选择px(图4-6)。
2)坡长的提取:坡长的提取,通常在ArcView软件中不是直接进行求取的,而是采用先求负地形,再通过其的水文分析功能,求出负地形的水流方向、水流长度等,而后得到正地形的山脊线、坡长。
a .激活DEM主题层。
b .在Analysis菜单下使用Map Calculator命令,公式为[[DEM-10000]×(-1)],提取DEM层的负地形,把DEM层的负地形记为A。
c .激活A层,调用Hydro菜单下的FlowDirection命令,则生成负地形的流向层,记为B。
d .激活B层,调用Hydro菜单下的FlowAccumulation命令,则生成负地形的水流累计量层Flow Accumulation,记为C,编辑C层的图例,使其值分为三个范围:0-固定值、固定值-最大值和NO DATA值,这三个值的颜色分别设为无填充色、任一彩色、黑色。这个固定值的确切值,可以通过C层和Hill shade层(通过)的共同显示来选定,选定固定值的标准是彩色值能比较好的反映山脊线,如图4-7。
图4-7 选定的固定值的图形显示
e .调用Analysis菜单下的Map Query命令,查询的表达式为:FlowAccumula tion≥“固定值”,得到新的主题 层Map Query 1,记为D,D是对C的二值化。
f .激活D层,点击Theme菜单下的Table命令或快捷按钮打开D主题的表,选择值为1的数据,再调用Analysis菜单下的Find Distance命令,则得到新的主题层,记为E。E上的每一个栅格的值是距最近的山脊线之间的垂直方向上的栅格数。
g .调用Analysis菜单的Calculator命令,公式为[E×栅格单元的尺寸],则得到的新的主题F的值为距最近的山脊线的垂直距离。(见图4-8)。
图4-8 坡长
因为水流的方向不是严格的和山脊线成90°,(如图4-9),大多数的水流方向只是接近90°,实际的坡长应是沿水流方向的长度,所以求得的主题F的值只是一种坡长的近似值,这种方法求坡长只是一种求取坡长的快速的、近似的方法。
图4-9 水流的实际方向
3)等高线的提取:利用DEM还可以生成等高线,我们使用ArcGis和Arc View的Contours功能可以生成一个新的线主题,每条线表示了具有相同高度、数量或者浓度的连续位置集合。生成的等值线经过平滑处理,真实地再现了表面等值线,这些等值线可以是等高线、等温线、等降水量线等。
等高线是地面上高程相同的各点连成的闭合曲线。根据等高线图形,可以判读地貌形态特征,量算各点的高程、坡向和坡度。
a . ArcView生成等高线:
ⅰ 在视图目录表中添加DEM并激活它。
ⅱ 从Surface菜单中选择Create Contours命令。
ⅲ 在出现的Contours Parameters对话框(如图4-10)中输入等高距Contour interval和基础等高线的值Base Contours。
图4-10 等高线制作对话框
ⅳ 生成等高线主题Contours of DEM(如图4-11)。
图4-11 生成等高线
生成单根等高线:则采用工具条上的等值线工具 可以生成单根的等值线,其步骤如下:
ⅰ 激活DEM,单激 按钮。
ⅱ 在DEM中点击选定的位置,视图上即描绘出一条通过所选点位代表所选点值大小的等高线,其高程值显示在窗口的左下角的状态栏中(如图4-12)。
图4-12 提取单根等高线
b . ArcGIS生成等高线
ⅰ 点击3D Analyst菜单,将鼠标指向Surface Analysis,选择Contour命令。
ⅱ 点击Input Surface下拉列表的箭头,选择用来生成等高线的DEM。在Contour interval文本框中输入等高距,一般为10m。在Bese contour文本框中输入等高线的基准高程,一般为最低高程或0。在Z factor文本框中指定高程变换系数。一般为1。保存等高线文件。单击OK(图4-13)。
图4-13 用ArcGIS生成的20m等高线图
(3)使用ArcGIS进行坐标系统定义
1)对于Coverage的定义:在Coverage Properties对话框的Projection(投影)选项卡,包含Coverage坐标系统及其投影参数信息。如果Coverage还没有定义坐标系统,可以通过单击投影选项卡中的Define按钮,为其定义一个坐标系统,定义投影向导会帮助我们完成投影坐标工作。在定义坐标系统时,既可以交互式地定义坐标系统的各参数,也可以将另外一个Coverage、Grid或TIN中的相关信息复制过来。如果Coverage已经定义了坐标系统,还可以在定义投影向导的帮助下改变现有的坐标系统。
a .交互式定义Coverage坐标系统:首先打开ArcCatalog模块,选择ArcCata log目录树中需要定义的Coverage文件,然后再用交互式定义其坐标系统。
ⅰ 在ArcCatalog目录树中选择需要定义的Coverage文件。
ⅱ 单击主菜单中的File菜单下拉菜单中的Properties命令。或在文件名上点击右键下拉菜单中的Properties命令。打开Coverage Properties对话框(见图4-14)。
图4-14 Coverage Properties对话框
ⅲ 单击Define按钮,打开定义投影向导选择方式对话框(图4-15)。
图4-15 定义投影向导选择方式
ⅳ 选择如图,进入交互方式定义投影坐标系统。单击Next按钮,打开定义投影向导确定投影类型对话框(图4-16)。
图4-16 根据定义向导确定投影类型
Ⅴ 在Propertions列表中选择投影类型,单击Next。输入需要的投影参数:坐标单位、投影分带、中央经线等(如图4-17)。
图4-17 确定投影坐标参数
ⅵ 单击Next按钮,确定坐标基准对话框。ArcGIS平台的所有预定义供用户选用的基准面(Datum)中,均没有我们国家的基准面定义。假如精度要求不高,可利用前苏联的Pulkovo1942基准面,代替北京54坐标系。椭球体(Spheroid)选用Krasovsky(图4-18)。
图4-18 确定投影基准
ⅶ 单击Next按钮,查看投影坐标系统,如果满意单击Finish按钮,返回Coverage Properties对话框,按确定按钮,完成坐标系统的定义。
b .用已知坐标系统来匹配未知坐标系统
首先打开ArcCatalog模块,选择ArcCatalog目录树中需要定义的Coverage文件,然后再来匹配方式定义其坐标系统。
ⅰ 在ArcCatalog目录树中选择需要定义的Coverage文件。
ⅱ 单击主菜单中的File菜单下拉菜单中的Properties命令。或在文件名上点击右键下拉菜单中的Properties命令。打开Coverage Properties对话框。
ⅲ 单击Define按钮,打开定义投影向导选择方式对话框。
ⅳ 选择第二选项,进入匹配方式定义投影坐标系统。单击Next按钮,打开定义投影向导确定已知投影系统的Coverage对话框(图4-19)。
图4-19 确定已知坐标系统的Coverage
Ⅴ 在Dataset文本框中单击浏览按钮,确定使用已知投影系统的Coverage。如图选择的china文件。查看出现的坐标系统参数。如果满意单击Finish按钮,返回Coverage Properties对话框,按确定按钮,完成坐标系统的定义。
2)对于Shapefile定义:首先打开ArcCatalog模块,选择ArcCatalog目录树中需要定义的Shapefile文件,然后再来定义其投影坐标系统。
ⅰ 在ArcCatalog目录树中选择需要定义的Shapefile文件。
ⅱ 单击主菜单中的File菜单下拉菜单中的Properties命令。或在文件名上点击右键下拉菜单中的Properties命令。打开Shapefil Properties对话框(见图4-20)。
图4-20 打开Shapefil Properties对话框
ⅲ 单击Fields选项卡,单击Fields Name列表中的Shape,打开Fields Proper ties列表中相应的字段特性列表,单击特性列表中Spatial Reference后的省略号(…)按钮,打开Spatial Reference属性对话框(图4-21)。
图4-21 Spatial Reference属性对话框
ⅳ 单击New按钮,选择Projected,打开定义投影坐标系统对话框(图4-22)。在Projection选项卡中单击Name文本框下拉箭头,选择一种投影类型,可以输入该投影的相关参数。在Geographic Coordinate System选项卡中,选择确定投影的大地基准面或椭球体。在该对话框中我们可以在Name文本框中输入定义的名称,如我们可以输入beijing;在Projection选项卡中单击Name文本框下拉箭头,选择投影类型为Gauss_Kruger;在Geographic Coordinate System选项卡中,选择GCS_beijing54。
Ⅴ 单击OK按钮,完成了投影坐标系统的定义。
图4-22 定义投影坐标系统对话框
(4)基于ArcGIS或ArcView的淤积面积和体积的计算
确定库区范围及淹没深度分布在进行流域蓄水工程规划中,特别是在进行大坝选址时,需要进行流域面积、不同蓄水位条件下水库淹没范围、库容等参数。现行的方法是找到适当比例尺的地图,人工在地图上量算。这是一项非常繁琐的工作,不仅效率低下,而且计算结果不够精确,会因人而异,借助于ArcView和DEM,这些工作变得简单多了。首先可以根据上述的流域识别方法快速地计算流域的各项特征值,结合水文模型计算后得到不同保证率下的坝顶高程。有了坝顶高程和流域范围后,在ArcView中就可以很方便得到淹没范围及淹没深度分布情况。步骤如下:
ⅰ 利用本文所开发的一小段程序,从原始DEM中裁减流域所覆盖范围下的栅格数据,我们称之为流域DEM。
ⅱ 利用空间分析模块中的MapCalculate功能菜单,在对话框中填入“流域DEM≤坝顶高程”,就生成了包含两个字段(value、count)及两条记录的矢量文件,其中value为0或1,分别表示不符合条件和符合条件。Count表示符合或不符合条件的单元数目。很明显符合条件的单元即为水库淹没范围,其面积是count值乘以单元面积。从上述矢量文件中挑选符合条件的记录,形成库区范围矢量文件。
ⅲ 以库区范围矢量文件边界,从流域DEM裁下新的栅格数据,形成库区DEM。再次利用空间分析模块中的MapCalculate功能菜单,在对话框中填入“坝顶高程——库区DEM”,就得到了库区淹没深度分布栅格图。
库容计算通过对库区DEM的分析和计算,很容易地得到不同坝高的库容。这部分工作需要编程来完成。本文利用在工作中编制的两个小程序,进行坝高和库容的计算。第1个程序是生成“平面DEM”。所谓“平面DEM”就是将库区DEM的有效单元值全部换成要计算的高程值。第2个程序库容计算。由于“平面DEM”和库区DEM的单元地理位置分布完全一样,将两个DEM中的对应单元高差乘以单元面积即得到该单元的所能储存的水量,遍历所有单元得到某一高程下的库容。计算多个高程值,就可以得到水库水位——库容关系曲线。
1)基于ArcGIS的淤地坝淤积面积和体积的计算:对于淤积面积的计算,在通过ArcMap自带功能,无需编程就可以从低于一特定高程区域中确定真正的淤积区,并可以直接计算出面积。
利用ArcMap的3D Analysis工具可以直接在DEM上,已知淤积厚度就可以直接计算出淤地坝的淤积面积和淤积的库容容积。其具体方法如下:
a .启动ArcMap模块,点击Tools下的Extensions…中的“3D Analyst”复选框来添加三维分析模块(图4-23)。
b .单击File菜单中的Add Data菜单添加用于分析的DEM层,在相应的路径下找到要计算流域的DEM数据,如我们找到王茂沟流域的数据wmgdem(图4-24)。
图4-23 Extehsions对话框
图4-24 添加wagdem进入ArcMap软件
c .然后启动3DAnalyst按钮,进而进入Surface Analysis菜单,然后单击Area and Volume Statistics按钮,显示Area and Volume Statistics对话框(图4-25)。
图4-25 计算面积与体积的属性图
d .从对话框的,Reference parameters属性表下,起始计算高程为(Height of plan)935.04m。Input height range下表示的是流域内最小和最大高程分别为935.04m和1188.15m。选择Calculate statistics above plane为从最低点向上计算。选择Calculate statistics below plane为从最高点向下计算。
e . Output statistics属性表表示的是2D area流域面积、Surface area表面积和Volume体积。每次变化Height of plan时,可以按Lculate statisti按钮重新计算。需要保存计算结果时,可以选择Save/append statistics to text file,将计算结果保存在文件名为areavol. txt的文本文件中。
2)基于ArcView测量面积和体积:使用三维分析可以测量表面面积和体积,同时还可以测量两个表面之间容积的差异——进行剪切、填充分析。
表面面积是沿表面的坡度进行测量,并将高度考虑在内,计算出的面积总要大于二维平面测量的面积。
体积是计算TIN表面和位于任何指定高程的水平面之间的立体空间,可以是平面之上的,也可以是平面之下的测量体积,在实际应用中一般用来计算土方量。
三维分析测量面积和体积,要在激活TIN主题下,调用Theme下的Convert Grid to TIN才可以实现(图4-26)。如果是Grid数据,则必须转换成TIN主题:
图4-26 Converf Grid to TIN
a .激活TIN主题;
b .点击Surface下的Area and Volume Statistics;
c .输入一个基础高程,这将作为测量面积和容积的水平平面;
d .指定要测量的是平面上还是平面下的面积和容积;
e .按OK确认;
则会出现一个对话框,显示二维面积、表面积及体积(见图4-27)。
图4-27 面积和体积统计的对话框
有时,还会分析比较两个表面模型前后的变化,这时需要调用Surface下的Cut Fill。其前提是激活前、后两个主题。
(5)利用ArcGIS软件实现王茂沟小流域三维虚拟景观
虚拟现实(Virtual Reality,简称VR)技术是一个图像技术、传感器技术、计算机技术、网络技术以及人机对话技术相结合的产物,它以计算机技术为基础并综合利用了计算机的立体视觉、触觉反馈、虚拟立体声等技术,高度逼真地模拟人在自然环境中的视、听、动等行为的人工虚拟环境。这种虚拟环境是通过计算机生成的一种环境,它既可是真实世界的模拟体现,也可以是构想中的世界。
小流域虚拟现实技术的重点研究和难点之处是具有高度真实感的虚拟三维图像的生成,其中包括三维地形的构造、纹理映射和实时动态立体显示等关键技术。
本次研究采用ArcGIS软件的三维可视化功能组件ArcScene程序的来实现王茂沟流域的三维景观浏览。具体制作过程见以下各节。
1)创建三维场景:在三维场景中浏览数据,将给人以耳目一新的感觉。在浏览相同的数据时,三维浏览可提供一些平面地图所无法直接提供的知识。如在三维场景中察看淤地坝坝系,可以清楚得看到各个淤地坝,在平面看影像时就达不到这个效果;在察看峡谷,并感觉山脊与山底在高程上的差异,而不必像看二维地形图时,要从等高线的分布情况进行推断。
a .添加数据:
ⅰ 添加DEM数据:首先打开ArcScene软件,单击Standard工具条上的Add Data按钮,选择王茂沟流域的DEM文件,点击Ddd(图4-28)。
图4-28 显示DEM三维显示
ⅱ 添加QuickBird影像数据:同样单击Standard工具条上的Add Data按钮,选择王茂沟流域的QuickBird影像文件,点击Ddd(图4-29)。
图4-29 显示QuickBird影像数据三维图
b .设置基准高程:
ⅰ 察看DEM设置基准高程:在ArcScene中,右击wmgdem图层,并选择Properties命令,弹出Properties对话框(图4-30),点击Base Heights标签。在Height中显示为0,可以看出该图层为没有基准高程。
图4-30 wmgdem图层属性对话框
ⅱ 设置影像基准高程:右击ji. img图层,并选择Properties命令,弹出Por perties对话框(图4-31),点击Base Heights标签。
图4-31 ji. img图层属性对话框
选择Height下的Obtain heights for layer from surface:选项,选择wmgdem为影像的基准高程,这样就在图层中将影像以三维的形式绘制出来了。Raster Resolution按钮可以改变栅格图的分辨率。Z Unit Conversion可以改变土层中高程和X、Y的单位的转换。为了让多个图层同时显示在一个三维可视化上,因此使用Offset(偏移图层中的高度值),这样在不同高度上看到不同属性值。
2)设置场景属性:在ArcScene中,可以设置某些属性,如垂直拉伸、动画旋转、背景颜色、场景范围、照明情况。这些属性对场景及其内的所有图层有效。用户还可以添加场景的注释以及设置其坐标系统。如果场景使用了多个浏览器,这些属性对所有的浏览器均有效。
a .改变场景的垂直拉伸:垂直拉伸可以用来改变表面的细微变化。在创建地形可视化时,如果表面的水平范围远远大于表面垂直变化的范围,此功能会非常有用。当垂直变化非常厉害时,可以使用分数作为垂直拉伸系数,是表面看上去稍平坦。通俗的说法就是当系数越大时山越高、沟越深,反之山越低、沟越浅。
垂直拉伸对场景内的所有图层均有效。用户也可以改变图层的高程转换系数,对单个图层进行垂直拉伸。
ⅰ 右击Scene layers,选择Scene Properties命令,在弹出的Scene Properties 对话框中,点击General标签(图4-32)。
图4-32 场景(图层)属性
ⅱ 点击Vertical Exaggeration下拉列表箭头,选择垂直拉伸系数,可以输入任意数字,我们这里就选择了3。
ⅲ 点击Calculate From Extent按钮,让系统根据场景范围与高程变化,自动计算垂直拉伸系数。
b .三维场景动画旋转:旋转场景是获取表面总体概况的一种好方法。使用动画旋转时,可以让场景围绕其中心旋转。可以调整旋转的速度与察看角度,并在场景旋转时,对其进行缩放。
ⅰ 激活动画旋转功能:同前选择Scene Properties对话框,点击General 标签。选择Enale Animated Rotation的复选框,激活动画旋转功能。此时,动画旋转功能激活后,场景浏览(Navigate)功能的鼠标指针周围增加了一个小圆圈。
ⅱ 动画旋转和旋转速度:单击场景浏览工具,点击场景,按住鼠标左键,移动鼠标,场景就会开始旋转。
改变旋转速度,当场景在旋转时,按Page Up键加快旋转速度;按Page Down键减慢旋转速度。
c .改变背景颜色:默认场景的背景色是白色的。用户可以改变背景色,以满足视觉的需要。各种色调的可以模仿自然界中的天空颜色,黑色可用来模仿黑夜。所以我们可以使用一种预设颜色值或自定义的颜色定义背景,也可以使用当前背景色为所有场景的默认背景色。
ⅰ 设置背景色:右击Scene layers,选择Scene Properties命令,弹出Scene Properties对话框,点击General标签。点击Background color彩色下拉列表的箭头选择适当的颜色,作为背景色。
ⅱ 设置默认背景色:右击Scene layers,选择Scene Properties命令,弹出Scene Properties对话框,点击General标签。选择Use as default in all new scenes复选框,将当前背景色设置为默认背景色。
d .改变场景的照明:我们可以通过设置光源的高度角、方位角以及对比度来绘制场景的照明情况。场景的照明属性对场景中的所有要素都有效(包括突出的多边形及线要素)。对于各个图层,用户可以在图层的Properties对话框Ren dering标签中,将图层阴影设置选项打开或关闭,以控制各个图层在照明时是否使用阴影。
ⅰ 设置光源的方位角:右击Scene layers,选择Scene Properties命令,弹出Scene Properties对话框,点击Illumination标签(图4-33)。
图4-33 光源属性
在Azimuth的文本框中,输入光源的方位角。方位角是指光线照射场景时,其指南针所显示的方向。
ⅱ 设置光源照明的高度角:在Altitude文本框中,输入光源照明的高度角。高度角是指光源的高度,以光源照身到表面时光线与水平面的夹角来度量。
ⅲ 设置光源照明的对比度:在Contrast文本框中,输入对比度值。用来控制表面上阴影数量的值。
e .改变场景范围:减少场景的范围,可以消除一些无关的信息,增加绘图时的性能。默认情况下,场景的范围为场景中所有图层的范围。用户可以改变场景的范围,使其与某个图层的范围一致,或通过X、Y坐标的最大最小值来指定。处于场景范围之外的数据将不会被显示。
ⅰ 设置场景范围来为某图层范围:右击Scene layers,选择Scene Properties命令,弹出Scene Properties对话框,点击Extent标签(图4-34)。
图4-34 场景范围属性
选择Layer(s)选项,点击下拉列表箭头,选择用来定义场景范围的定义。
ⅱ 使用坐标设置场景范围:选择Custom相项,在文本框中输入X、Y坐标的最大最小值,定义场景。
3)场景输出:输出功能可以将场景的二维影像输出到外部图形文件中,也可以输出为三维虚拟现实建模语音(VRML)。场景的影像可用几种文件格式存储,并在其他文件中使用,在地图或报告中。
a .输出场景的二维图形:
ⅰ 点击File菜单,将鼠标指向Export Scene并选择2D。浏览到要存储场景的影像文件夹。
ⅱ 点击文件类型,选择适合于应用的输出文件格式。指定输出的栅格像元的宽度。指定输出影像的文件名,点击Export。
另外,可以通过使用Edit菜单中的Copy scene to clipboard命令,将场景复制到剪贴板上,直接可以用于文件。
b .输出场景到三维VRML模型中
ⅰ 点击File菜单,将鼠标指向Export Scene并选择3D。浏览到要存储VRML模型的文件夹。
ⅱ 输入VRML模型的文件名,点击Export,场景输出到VRML文件中。
c .打印场景
当需要场景的硬拷贝时,可以打印场景。我们可以将打印到默认打印机中,也可以改变页面设置,选择其他打印机,并指定打印的具体位置。
ⅰ 点击Print工具,单击Setup。在Name下拉列表中选择打印机。
ⅱ 点击Printer Page Size下拉列表,选择页面尺寸。点击Portrait或Land scape选项,选择页面方向。点击Printer Engine下拉列表,选择打印引擎。点击OK。
(1)信息源的选择
1)QuickBird信息源:本项目在影像选择是拟采用QuickBird或IKONOS影像作为信息源,在开展工作后了解到两种数据价格差异不大且QuickBird空间分辨率达到0.61m,而IKONOS空间分辨率为1m,所以我们选择了QuickBird数据。
QuickBird(快鸟)是美国Digital Globel公司于2001年10月18日发射Earth Watch卫星搭载的传感器,它共有5个波段,其中4个多谱波段,1个全色波段,光谱范围为450nm~900nm,空间分辨率最高可达0.61m。其遥感数据是商用卫星最高空间分辨率的一种影像。QuickBird传感器主要成像参数如表4-18。
表4-18 QuickBird传感器主要成像参数表
续表
2)课题涉及的流域范围:本课题涉及研究范围为王茂沟流域,王茂沟流域仅为5.97km 2 ,但购买的影像源最小定购面积为:64 km 2 ,所以在信息源购买的时候同时考虑了今后研究需要,将我站的王茂沟流域、韭园沟流域、桥沟流域及试验场的辛店沟流域的三沟一场都考虑在其中,共计面积为127.6 km 2 。其范围为:
NW Lat=37.65677600;
NW Long=110.23860900;
SE Lat=37.49250000;
SE Long=110.41239900;
拍摄时间为2004年9月9日。
(2)QB影像处理及纠正
1)图像分幅裁剪:由于王茂沟小流域的影像在购买的影像一部分,所以为了提供软件的运行速度,将王茂沟小流域的影像裁剪下来,组成一块小的影像,便于运行和管理。
图像的分幅我们应用ERDAS IMAGINE软件中不规则分幅裁剪功能,即裁剪图像的边界范围是任意多边形(接近小流域的形状),这样就可以事先沿小流域的流域界生成的一个完整闭合多边形区域裁剪,也可以是ARCINFO的一个Polygon Coverage流域界,针对不同的情况采用不同裁剪过程。
a .在ERDAS IMAGINE中生成多边形裁剪:在本次遥感处理中需要将中科院卫星地面站提供的影像进行裁剪,选择04SEP09033059-S2AS_R1C2-000000138971_01_P003. TIF文件将王茂沟流域分割出来,并保存为输出数据类型为Unsigned 8 Bit的ERDAS IMAGINE格式文件wmgz. img。
图4-35 用AOI工具绘制图像
在ERDAS IMAGINE软件菜单中打开要裁剪的文件:g:04SEP09033059-S2AS _R1C2-000000138971 01 _P003. TIF,并应用AOI的多边形工具沿小流域的分水岭绘制小流域界线,如图4-35。然后将多边形AOI保存在文件中(*. aoi)。或可以暂时不退出Viewer窗口,将图像与AOI多边形保留在窗口中。
在ERDAS图标面板工具条中单击DataPrep图标,打开Data Preparation弹出框,单击Subset Image按钮,打开Subset弹出框。
需要设置弹出框的内容如下(见图4-36):
图4-36 Subset对话框
ⅰ 输入原影像文件名称(Input File)为:04SEP09033059-S2AS_R1C2-000000138971_01_P003. TIF。
ⅱ 输出裁剪后文件名称(Output File)为:wmgz. img。
ⅲ 单击AOI按钮确定裁剪范围。打开Choose AOI弹出框,在其中确定为AOI的来源(AOI Source)为File,选择*. aoi文件。不退出Viewer窗口的可以选择Viewer。
ⅳ 输出数据类型(Output Data Type)为Unsigned 8 Bit。
Ⅴ 输出象元波段(Select Layers)为1∶3(表示选择1、2、3三个波段)。
ⅵ 单击OK按钮完成裁剪工作。
b . ArcInfo的一个Polygon Coverage裁剪:在水土保持工作中,一般情况下按照流域界线、行政区域或分区划分对图像进行分割,所以往往是过去就有现成的边界线,如ArcInfo或其他矢量的多边形图,然后可以利用这些图来裁剪图像,实现同形状、同大小的影像和区域。这里就我们使用ArcInfo的Polygon为边界进行影像的裁剪。用王茂沟流域ArcInfo格式界线wmg1对影像wmg. img进行裁剪。
ⅰ 将ArcInfo多边形的矢量文件转换为栅格文件格式
在REDAS的图标面板工具条单击Interpreter图标,在Interpreter弹出框中点击Utilities菜单下的Vector to Raster,打开Vector to Raster弹出框(图4-37)。
图4-37 Vector to Rafter对话框
弹出框需要设置的内容如下:
●输入ArcInfo矢量文件名称(Input Vector File)为wmg1。
●输出栅格文件名称(Output Image File)为wmgsg. img。
●确定矢量文件类型(Vector Type)为Polygon。
●栅格数据类型(Data Type)为Unsigned 8 bit。
●使用矢量属性值(Use Attribute As Value)为wmgl-id。
●栅格文件类型(Layer Type)为Thematic。
●转换范围大小(Size Definition),在ULX、ULY、LRX、LRY微调框中输入需要的值。
●坐标单位(Units)为Meters。
●单击OK按钮。
ⅱ 利用掩膜运算(Mask)实现图像不规则裁剪
在ERDAS图标面板工具条单击Interpreter图标,在Interpreter弹出框中点击Utilities菜单下的Mask,打开Mask弹出框(见图4-38)。
图4-38 掩膜对话框
弹出框需要设置的内容如下:
●输入图像文件名称(Input File)为wmg. img。
●输入掩膜文件名称(Input Mask File)为wmgsg. img
●输出图像文件名称(Output File)为wmgz. img。
●单击Setup Recode设置裁剪区域内新值(NewValue)为1,区域外取0值。
●确定掩膜区域做交集运算为Intersection。
●输出数据类型(Output Data Type)为Unsigned 8 bit。
●输出统计忽略零值即选中Ignore zero In Output Stats复选框。
●单击OK按钮(关闭Mask弹出框,执行掩膜运算)。
2)图像拼接:由于QuickBird卫星影像的成像模式为单景16.5km×16.5km。所以在研究区域内大于此面积和跨几景的情况下,需要对这些具有地理参考的若干相邻图像进行合并。这些要拼接的图像必须具有地图投影信息,或者说这些图像必须经过几何校正处理。同时它们必须具有相同的波段数。在进行图像拼接时,需要确定一幅参考图像将作为输出拼接图像的基准,决定拼接图像的对比度匹配、以及输出图像的地图投影、象元大小和数据类型。
a .图像拼接的加载:
ⅰ 图像的加载:首先启动ERDAS图标面板工具单击DataPrep图标,打开Data Preparation弹出框,单击Mosaic Images按钮,打开Mosaic Tool窗口。
其次加载需要拼接的图像,这里我们需要加载wmg1. img和wmg2. img两格图像进行拼接。
下面实现加载的过程,在Mosaic Tool菜单条单击Edit下的Add Images命令,打开Add Images for Mosaic弹出框。
弹出框需要设置的内容如下(见图4-39):
图4-39 Add Images for Mosaic对话框
●输入拼接图像文件名称(Image File Name)为wmg1. img。
●设置图像拼接区域(Image Area Option)为Use Entire Image。
●单击Add按钮(图像wmg1. img被加载到窗口中)。
●重复前三步操作,加载wmg2. img。
●单击Close按钮(关闭Add Images for Mosaic弹出框)。
ⅱ 图像的叠放次序调整:进行图像叠放次序,单击Mosaic Tool工具条Set Input Mode图标,并在窗口内选择需要调整的图像,进入设置输入图像模式的状态,Mosaic Tool工具条中会出现与该模式对应的调整图像叠放次序的编辑图标,充分利用系统所提供的编辑工具,根据需要进行上下层的调整。其中Send Image to Top选择置图像于最上层,Send Image Up One选择置图像上移一层,Send Image to Bottom选择置图像最下层,Send Image Down One选择置图像下移一层,Reverse Image Order选择图像次序颠倒。调整完成后,在Mosaic Tool 窗口单击就退出图像叠置组合状态。
b .图像的匹配及拼接:
ⅰ 图像匹配:在Mosaic Tool窗口的菜单条单击Edit后,在点击Image Matching命令,打开Matching Options弹出框。选择弹出框中的内容中匹配方式(Matching Method)为Overlap Areas(重叠范围匹配)。然后单击OK按钮(保存设置,关闭Matching Options弹出框)。
ⅱ 设置重叠功能:在Mosaic Tool窗口的菜单条单击Edit后,在点击SetOverlap Function命令,打开Set Overlap Function弹出框。选择弹出框内容中的设置交集类型(Intersection Method)为No Cutline Exists(没有裁割线),设置重叠区象元灰度计算(Select Function)为Average(平均值)。单击Apply按钮(保存设置),然后单击Close按钮(关闭Matc hing Options弹出框)。
ⅲ 图像拼接:在Mosaic Tool窗口的菜单条单击Process后,在点击Image Matching命令,打开Matching Options弹出框。
图像拼接的弹出框内容设置为(见图4-40):
图4-40 Run Mosaic对话框
●确定输出的拼接文件名称(Output File Name)为wmgxin. img。
●确定输出图像区域(Witch Outputs)为All。
●忽略输入图像值(Ignore Input Value)为0。
●输出图像背景值(Output Background Value)为0。
●忽略输出统计值,即选择Stats Ignore Value复选框。
●单击OK按钮(关闭Run Mosaic弹出框,运行图像拼接)。
3)图像几何精校正:遥感图像的几何校正可分两阶段实现:系统校正(几何粗校正),即把遥感传感器的校准数据、传感器的位置、卫星姿态等测量值代入理论校正公式进行几何畸变校正;几何精校正,即利用地面控制点GCP(GroundControl Point,遥感图像上易于识别,并可精确定位的点)对因其他因素引起的遥感图像几何畸变进行纠正。几何粗校正的由卫星接收地面站负责。因此下面重点讨论几何精校正。
几何精校正原理:从物理上看,畸变就是像素点被错误放置,即本该属于此点的像素值却在彼处。因此可用两种方法实现畸变图像的校正:一是把被错置的像素点搬运到该在的位置,此方法被称为直接变换法;二是取回属于该位置的像素值,此方法被称为重采样法。从数学上说,几何精校正就是通过选取一组GCP建立原始的畸变图像空间与校正图像空间的坐标变换关系:
或
式中 x , y ——畸变图像空间中的象元坐标;
ξ , η —— x , y 在校正图像空间中对应的象元坐标,称作 x , y 的共轭点。
下面以较常采用的重采样法为例,说明几何精校正过程。重采样法的几何精校正过程包含如下两方面工作。
1.确定 p -1、 q -1的函数形式。
一般作法是用二维 m 阶多项式逼近 p -1、 q -1即:
然后,根据一组已知像素点(即一组GCP)的对应坐标( ξi , ηi )和( xi , yi )确定 ajk 、 bik 。
2.确定( ξ , η )的灰度值
校正空间像元( ξ , η )的灰度值 g ( ξ , η )等于原始空间共轭点( x , y )的灰度值 f ( x , y ),通常不是象元点,因此( x , y )的灰度值需由( x , y )邻近像元点的灰度值内插获得。通常采用的方法有:最近邻点法、双线性内插法和三次褶积法。
由于我们使用的信息源来自于多个渠道,尤其是图像数据来源于卫星,DEM依据的数据是20世纪50年代1∶10000地形图,GPS数据的地形数据也有所不同,所以图像纠正是至关重要的一个环节。
我们采用的纠正是依据DEM数据和GPS监测数据对我们的QuickBird影像进行纠正,选择的计算模型为多项式变换(Polynomial),在调用多项式模型时,需要确定多项式的次方数(Order),我们选择图像为3次方,共需要10个控制点。
运用ArcView软件对已有的DEM进行1m等高线提取,保存为shape格式即可。
(1)显示图像及DEM数据
1)打开ERDAS,并在图标面板中双击Viewer图标,同时打开两上窗口(Viewer#1/Viewer#2),并将两上窗口平铺放置(图4-41)。
2)在Viewer#1窗口中打开需要校正的QB图像wmg. img。
3)在Viewer#2窗口中打开作为地理参考的shape文件dgx. shape。
(2)启动几何纠正功能模块
1)在Viewer#1菜单条上单击Raster Geometric Correction命令(图4-42)。
图4-41 打开ERDAS后需要校正的QB图像
图4-42 ERDAS IMAGINE窗口
2)打开Set Geometric Model对话框,选择多项式几何校正计算模型:Polynomial,单击OK按钮,即打开Polynomial Model Properties窗口,定义多项式模型参数及投影参数(多项式次方[Polynomial Order]及投影参数[projection]均可根据具体要求所定义),单击Apply按钮应用或单击Close按钮关闭图4-43。
(3)采集DEM控制点
1)打开GCP Tool Reference Setup对话框选择采点模式(图4-43),即选择Existing Viewer单选按钮,用现有的数据进行校正。
图4-43 Polynomial Model properties窗口
图4-43 GCP Tool Reference对话框
2)单击OK按钮,关闭GCP Tool Reference Setup对话框。
3)打开Viewer Selection Instructions指示器(图4-44)。
图4-44 Viewer指示器
4)在显示作为地理参考图像wmgxin. img的Viewer#2中单击。
5)打开Reference Map Information提示框(图4-45)。
图4-45 显示参考图像的投影信息
6)单击OK按钮,关闭Reference Map Information对话框。
7)整个屏幕将自动变化为如图4-46所示的状态,表明控制点工具已被启动,进入控制点采集状态。
图4-46 控制点进入采集状态实况显示
(4)采集点地面检查点
1)在GCP工具对话框中单击Select GCP图标 ,进入GCP选择状态。
2)在GCP数据表中将输入GCP的颜色Color设置为比较明显的黄色。
3)在Viewer#1中移动关联方框位置,寻找明显的地物特征点,作为输入的GCP。
4)在GCP工具对话框中单击Create GCP图标 ,并在Viewer#3中单击定点,GCP数据表将记录一个输入GCP,包括其编号、标识码、X坐标、Y坐标。
5)在GCP工具对话框中单击Select GCP图标 ,重新进入GCP选择状态。
6)在GCP数据表中将参考GCP颜色设置为比较明显的红色。
7)在Viewer#2中移动关联方框位置,寻找明显的地物特征点,作为参考的GCP。
8)在GCP工具对话框中单击Create GCP图标 ,并在Viewer#4中单击定点,系统将自动把参考点的坐标显示在GCP数据表中。
9)在GCP工具对话框中单击Select GCP图标 ,重新进入GCP选择状态;并将光标移回到Viewer#1,准备采集另一个输入控制点。
10)不断重复步骤(1)~(9),采集若干个GCP,直到满足所选定的几何校正模型为止。而后,每采集一个Input GCP,系统就自动产生一个Ref. GCP,通过移动Ref. GCP可以逐步优化校正模型。
采集GCP以后,GCP数据表如图4-47所示。
图4-47 采集点的编号、标识码及其相关坐标
采集GCP点具体位置图如图4-48所示。
图4-48 采集点具体位置图
(3)建立解译标志及图像解译
本项目采用的遥感数据是0.61m分辨率的QuickBrid卫星影像,使用的数字高程模型(DEM)数据我们制作的1∶1万Grid。根据卫星影像纹理及色彩,提取王茂沟流域的土地利用现状及植被覆盖度信息,利用DEM数据提取坡度信息,依据前三类信息综合判断土壤侵蚀等级,以此调查王茂沟流域的水土保持现状。
1)技术路线:以王茂沟流域的数字高程模型(DEM)及地形图为基础,利用遥感影像处理软件对影像进行纠正、调色、裁剪、拼接等处理;通过外业调查,建立影像与实地对应的解译标志;依据解译标志在影像上提取土地利用及植被覆盖度信息,并建立相关矢量图层;利用DEM数据根据栅格数据空间分析获得坡度信息,并生成坡度矢量图层;结合土壤侵蚀分级指标,在已有三类信息的基础上,进行矢量图层叠加,并计算各划分单元的土壤侵蚀强度等级(如图4-49所示)。
图4-49 遥感解译技术流程框图
2)信息源:本次遥感调查的主要信息源有王茂沟流域的QuickBrid卫星影像、1∶1万的DEM和地形图等。
QuickBrid卫星影像拍摄于2004年9月9日,是分辨率为0.61m的真彩色影像,可直观地从影像中解译出各种土地利用类型、植被覆盖度等信息,能够较好地满足小流域水土保持遥感调查的要求。
QuickBrid为8-bit真彩色卫星影像,在解译前进行了文件格式转换、色彩调整等处理,然后在Arc/Info中进行解译。
1∶1万DEM数据栅格分辨率为5m,利用Arc/Info软件的空间分析生成栅格数据形式的坡度图,再将栅格数据转为矢量数据,利用它执行Arc/Info的Union命令,将坡度数据与其他数据进行叠加。
3)影像解译标志的建立:影像解译标志是内业解译的依据,在解译工作开始之前,必须结合影像到实地建立影像色调、纹理结构、颜色、形状与土地利用类型、植被覆盖度的对应关系,这样才可在室内准确地从影像中进行判读,保证解译数据精度(见表4-19)。
表4-19 王茂沟流域Quickbird卫星影像解译标志
本次外业工作的主要内容有:
a .通过外业调查,建立影像解译标志。拍摄相应的野外实况照片,利用GPS准确定位,建立影像特征与实地的对应关系,供内业解译参考。
b .利用1∶1万地形图进行适量的调绘,为内业解译建立判读的认识基础。
c .进行重点淤地坝调查,采用 M s = V s /( F × N )公式计算侵蚀模数( M s 为平均侵蚀模数, V s 为淤积量, F 为控制面积, N 为淤积年限),为内业判读提供参考数据。
d .通过收集资料和走访当地群众进行社会经济状况调查。
4)图斑编码:根据王茂沟流域的实际情况,本次解译采用五位编码,第一、二位码为土地利用类型,第三位码为坡度,第四位码为植被覆盖度,第五位码为侵蚀强度。
5)影像解译:采用人机交互解译法。质量控制要求如下:
a .影像处理后,图像清晰,色调均一,反差适中,阴影、云彩或其他斑痕不影响应用。影像精细纠正后,误差控制在3个像元内。
b .解译图斑属性的判对率大于90%。
c .图斑边界线的走向和形状与影像特征之间的允许误差为1个像元。
d .农田、林地等最小图斑面积为100个像元,荒草地为400个像元。条状图斑短边长度不小于5个像元。
ⅰ 综合判读:根据野外建立的解译标志、影像以及相关资料在Arc/Info下综合判读,按照土地利用类型和植被覆盖度分级情况勾绘图斑,并赋前三位码。农田的植被覆盖度定为4级。
ⅱ 添加坡度属性码:在Arc/Info的空间分析模块下,利用DEM生成矢量坡度图,再与综合判读图叠加,添加坡度属性码。水面、梯田、坝地以及川台地坡度定为0 °~5 °。
ⅲ 添加侵蚀强度属性码:根据前四位码(即土地利用类型、植被覆盖度和坡度),人工交互为每个图斑添加侵蚀强度属性码。人工判别的标准见表4-20。
表4-20 水力侵蚀强度分级参考指标
6)侵蚀强度等级判读:小流域的土壤侵蚀模数可以通过解译成果中所有划分单元(图斑)的面积和土壤侵蚀强度推算出来。推算出来的土壤侵蚀模数与野外实地调查获得的土壤侵蚀模数比较,即可验证解译成果中划分单元(图斑)的土壤侵蚀强度赋值的总体精度。
推算过程中,各侵蚀强度等级对应的侵蚀模数参考指标采用部颁标准——土壤侵蚀分类分级标准(SL 190-96),通过验算,解译成果侵蚀模数计算结果为6206 t/km 2 ·a,而1999年黄河流域遥感普查项目对王茂沟流域进行野外实地调查获得的侵蚀模数为8700 t/km 2 ·a。经过2000年至2005年韭园沟示范区建设王茂沟流域的坡面治理度达到73.2%,又经过沟口监测2000年到2005年5年内流域共计输出泥沙10t/km 2 。因此我们认为土壤侵蚀模数为6206 t/km 2 ·a是比较准确的(见表4-21)。而现在采用的1973年出的《榆林地区实用水文手册》(内部资料)由于时间较早,地形地貌、植被、气候等影响水土流失的因素发生了变化,因此,应该对已经使用了30多年的水文手册进行必要的修订。
表4-21 遥感解译侵蚀强度等级判读精度验算表
本次遥感影像解译中,不仅影像分辨率高,而且采用了叠加DEM生成坡度图和人机交互判定侵蚀强度等技术,使影像判读过程中误差大大降低,客观性增强,判读准确率有很大提高,增强了王茂沟小流域淤地坝建设的科技含量,具有一定的示范作用。
通过本次遥感调查,建立了一套QuickBrid影像解译标志和王茂沟流域遥感解译图形、图像属性库,并制作了王茂沟流域影像图、水土保持措施现状图、坡度图和骨干坝、中型淤地坝现状及规划布设图等水土保持系列图件。从解译结果可以看出,流域右岸植被覆盖较好(为流域的阴坡)、土壤侵蚀相对较轻,其他区域的土壤侵蚀严重,小流域的土壤侵蚀模数为6206 t/km 2 ·a,基本符合流域的实际。
(4)RS与DEM的淤地坝淤积量方法
用RS采集坝系的面积图,基于ERDAS IMAGINE的淤地坝淤积面积计算。
1)淤积范围图的制作:淤积范围的确定,(监测时间前后)是暴雨洪水发生前后的淤积对比过程,它包括三个步骤:
a .在淤积监测前淤地坝地界定。
b .淤积后地坝地面积地界定。
c .将上述两面积图叠加求面积范围与面积。
2)数字高程模型(DEM)的制作与淤积厚度的计算:利用该区1∶10000的地形图制作的数字高程模型(DEM),与淹没范围图叠加,可以获取淹没水深分布图。也可用来进行洪水演进,得到淹没水深分布与淹没范围图。
利用ERDAS软件的Virtual GIS分析功能中的Water Layer层(洪水淹没分析)来计算淤地坝的淤积面积和淤积体积。
在Virtual GIS窗口中可以叠加洪水层(Overlay Water Layers),进行洪水淹没状况分析,即亦可以利用此功能进行淤地坝的淤积面积和淤积体积的计算。系统提供了两种分析模式(Fill Entire Scene和Create Fill Area)进行操作。