多年以后，发现自己的帖子还没回答。于是找到当年自己的成果，把它贴出来，大家指点指点。

——为记

load exampledate1.txt

x_data=exampledate1(:,1);

y_data=exampledate1(:,2);

z_data=exampledate1(:,3);

%% 筛选指定区域的样本点 %%

N=length(x_data);

j=0;

for w = 1:N

if (84740<=x_data(w) & x_data(w)<= 84880)&(89300<=y_data(w) & y_data(w)<=89500)

j=j+1;

%% 将选中的样本点存入数组 %%

x_new(j)=x_data(w);

y_new(j)=y_data(w);

z_new(j)=z_data(w);

end

end

x=x_new()';

y=y_new()';

z=z_new()';

%% 搜索圆内参考点个数num %%

num=7;

%% 采样点总数n %%

n=length(x);

%% 区域面积A %%

A=(max(x)-min(x))*(max(y)-min(y));

%% 搜索圆半径 %%

R=sqrt(num*A/(n*pi));

%% DEM分辨率 %%

delt_x=5;

delt_y=5;

%% 格网点的个数 x轴方向：mm；y轴方向：nn %%

%% ceil函数是向正无穷方向取整 %%

mm=ceil((max(x)-min(x))/delt_x);

nn=ceil((max(y)-min(y))/delt_y);

%画格网

xx = [min(x):delt_x:max(x)];

yy = [min(y):delt_y:max(y)];

[X,Y]=meshgrid(xx,yy);

figure(1)

clf

hold on

plot(x,y,'o')

plot(X,Y,'+')

xlabel('X(m)')

ylabel('Y(m)')

title('Grid and observation points')

axis ([min(x),max(x),min(y),max(y)])

%% 待插值点的点的坐标 %%此处，0.5*delt_x的0.5为自己设定的值。根据实际情况可更改。

x0=min(x)+0.5*delt_x;

y0=min(y)+0.5*delt_y;

for u = 1:mm

for v = 1:nn

%% 计算待插值点坐标x_grid(u,v),y_grid(u,v) %%

x_grid(u,v) = x0 + (u-1)*delt_x;

y_grid(u,v) = y0 + (v-1)*delt_y;

j=0;

%% 所有样本点到待插值点的距离d(w) %%

for w = 1:n

d(w)=sqrt((x_grid(u,v)-x(w))^2+(y_grid(u,v)-y(w))^2);

if d(w)-R<=0  %% 判断取样点是否在搜索圆区域内 %%

j=j+1;

%% 将选中的样本点存入数组 %%

x_select(j)=x(w);

y_select(j)=y(w);

z_select(j)=z(w);

end

end

%% 当取样点个数少于7个时，扩大搜索半径 %%

R_new=R;

while j<7

R_new=R_new+0.5*delt_x;

k=0;

for w = 1:n

if d(w)-R_new<=0

k = k+1;

x_select(k)=x(w);

y_select(k)=y(w);

z_select(k)=z(w);

end

end

j = k;

end

%% 当取样点个数大于10个时，缩小搜索半径 %%

R_new=R;

while j>10

R_new=R_new-0.5*delt_x;

k=0;

for w = 1:n

if d(w)-R_new<=0

k = k+1;

x_select(k)=x(w);

y_select(k)=y(w);

z_select(k)=z(w);

end

end

j = k;

end

% 截除冗余数据 %%

x_select = x_select(1:j);

y_select = y_select(1:j);

z_select = z_select(1:j);