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

水系提取算法（D8）的matlab 实现

D8算法是当今非常成熟的提取水系（或分水岭）的计算机实现方法。这里笔者结合自身的编程实践介绍D8算法的matlab 实现过程，谨与大家分享。

D8算法分为两部分实现：每个栅格点水流方向的计算；每个栅格点汇水面积的计算。

1、水流方向的计算

在水流方向的计算中，本人结合matlab 内嵌的滤波工具，利用一个简单的非线性滤波来实现，具体内容包括一个自己编辑的确定水流方向的函数flowdirection 和若干脚本命令，如下：

以下是自己编辑的flowdirection 函数代码：

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); 1 4 7 2 5 8 3 6 9

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

f(i+2,j+1,s)=f(i+2,j+1,s)+f(i+1,j+1,s-1); elseif a(i,j)==7

f(i,j+2,s)=f(i,j+2,s)+f(i+1,j+1,s-1);

elseif a(i,j)==8

f(i+1,j+2,s)=f(i+1,j+2,s)+f(i+1,j+1,s-1); else

f(i+2,j+2,s)=f(i+2,j+2,s)+f(i+1,j+1,s-1); end

end

end

end

f=f(2:m+1,2:n+1,:);

该函数的执行结果为一个三维数据体，其中f（:,:,1）、f(:,:,2)、。。。。f(:,:,x)依次代表了不同窗口下每个栅格点的汇水面积，可以根据对比他们的效果并且根据需要选择不同的数据层。需要帮助1045105061@。

以下图片是本人提取的水系的效果图，仅供参考。