自然水系(或分水岭)自动提取的matlab实现方案

水系提取算法（D8）的matlab实现D8算法是当今非常成熟的提取水系（或分水岭）的计算机实现方法。这里笔者结合自身的编程实践介绍D8算法的matlab实现过程，谨与大家分享。D8算法分为两部分实现：每个栅格点水流方向的计算；每个栅格点汇水面积的计算。1、水流方向的计算在水流方向的计算中，本人结合matlab内嵌的滤波工具，利用一个简单的非线性滤波来实现，具体内容包括一个自己编辑的确定水流方向的函数flowdirection和若干脚本命令，如下：以下是自己编辑的flowdirection函数代码：functionoutput=flowdirection(a)n=size(a);fori=1:n(2)k=0;b=-inf;forj=[1:4,6:9]ifrem(j,2)==0r(j)=a(5,i)-a(j,i);elser(j)=(a(5,i)-a(j,i))/sqrt(2);endifr(j)bb=r(j);k=j;endendoutput(i)=k;end执行的脚本代码如下：uiopen('F:\课件\2009田淑芳\正射校正\DEM.tif',1)%打开DEM高程数据，生成的数据矩阵名可自行定义，这里默认为DEMDEM=mat2gray(DEM);%高程数据归一化处理DEM=padarray(DEM,[11],'replicate');%对数据边缘进行自动填充direc8=colfilt(DEM,[33],'sliding',@flowdirection);%执行非线性滤波，得到方向矩阵，得到的值分别为{1,2,3,4,6,7,8,9}.其方向意义如上图矩阵所示，如1代表北西向，8代表正东1、汇水面积的计算这里笔者根据自己的编程实践经验，给出一个函数来计算不同窗口大小下汇水面积的大小。函数代码如下：functionf=flowareax(a,x)%参数a为水流方向矩阵，x为表征窗口大小的参数，b=a;[mn]=size(b);147258369f=ones([m+2n+2x]);%默认每点有一初始水量fori=1:mforj=1:nifa(i,j)==1f(i,j,1)=f(i,j,1)+1;elseifa(i,j)==2f(i+1,j,1)=f(i+1,j,1)+1;elseifa(i,j)==3f(i+2,j,1)=f(i+2,j,1)+1;elseifa(i,j)==4f(i,j+1,1)=f(i,j+1,1)+1;elseifa(i,j)==6f(i+2,j+1,1)=f(i+2,j+1,1)+1;elseifa(i,j)==7f(i,j+2,1)=f(i,j+2,1)+1;elseifa(i,j)==8f(i+1,j+2,1)=f(i+1,j+2,1)+1;elsef(i+2,j+2,1)=f(i+2,j+2,1)+1;endendend%通过该循环确定经过一步流动时某一点水流量fors=2:xfori=1:mforj=1:nifa(i,j)==1f(i,j,s)=f(i,j,s)+f(i+1,j+1,s-1);elseifa(i,j)==2f(i+1,j,s)=f(i+1,j,s)+f(i+1,j+1,s-1);elseifa(i,j)==3f(i+2,j,s)=f(i+2,j,s)+f(i+1,j+1,s-1);elseifa(i,j)==4f(i,j+1,s)=f(i,j+1,s)+f(i+1,j+1,s-1);elseifa(i,j)==6f(i+2,j+1,s)=f(i+2,j+1,s)+f(i+1,j+1,s-1);elseifa(i,j)==7f(i,j+2,s)=f(i,j+2,s)+f(i+1,j+1,s-1);elseifa(i,j)==8f(i+1,j+2,s)=f(i+1,j+2,s)+f(i+1,j+1,s-1);elsef(i+2,j+2,s)=f(i+2,j+2,s)+f(i+1,j+1,s-1);endendendendf=f(2:m+1,2:n+1,:);该函数的执行结果为一个三维数据体，其中f（:,:,1）、f(:,:,2)、。。。。f(:,:,x)依次代表了不同窗口下每个栅格点的汇水面积，可以根据对比他们的效果并且根据需要选择不同的数据层。需要帮助[e-mail]。以下图片是本人提取的水系的效果图，仅供参考。