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

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