最近研究汇流累积量的计算，下了些资料，不过都是英文的看着真的很吃力。
搞明白了其中一种方法，是用Matlab写的：
78 72 69 71 58 49
74 67 56 49 46 50
69 53 44 37 38 48
64 58 55 22 31 24
68 61 47 21 16 19
74 53 34 12 11 12
1.算法初始值dem高程值矩阵、x坐标值矩阵和y坐标值矩阵。
2.返回汇流累积量flowacc、坡度slope。
首先计算坡度，采用邻接矩阵36*36表示。
[ic1,icd1] = ixneighbors(dem);
e = (dem(ic1)-dem(icd1))./hypot(X(ic1)-X(icd1),Y(ic1)-Y(icd1));
接着计算汇流方向，这里仍采用邻接矩阵表示，没有采用类似ESRI的D8编码方法。single方式，根据计算8个方向的坡度取最大的，另外一种是Multi方式
还有边缘处理、sink点处理等.
最后采用求多元方程组解的方式求得汇流累积向量flowacc，很神奇，原来矩阵运算可以求累积。
switch params.edges
case 'open';
flowacc = (speye(nrc,nrc)-spdiags(edgecorrection,0,nrc,nrc)*M')params.W0(:);
otherwise
flowacc = (speye(nrc,nrc)-M')params.W0(:); %36*1=（36*36-36*36）36*1，表示(speye(nrc,nrc)-M')*flowacc = params.W0(:);
end
flowacc = reshape(flowacc,siz);
这种方法最后得到的数据比正常的结果都增加了1。
1 1 1 1 1 1
1 2 2 3 3 1
1 4 8 6 5 1
1 1 1 21 1 2
1 1 1 2 25 1
1 3 5 8 36 2
参考D8编码方法：
