www.gusucode.com > 《matlab图像处理与界面编程宝典》秦襄培 编著，每章的MATLAB源代码程序 > 第3章/第3.4.5节中的代码.txt

     
% Import Data Files
loc = importdata('loctim.txt');
bound = importdata('border.xy');
lat = loc(:,1);     %纬度
lon = loc(:,2);     %经度
ele = loc(:,3);     %拔高度
tim = loc(:,4);     %地震运动到达时间
% Plot Switzerland
plot(bound(:,1),bound(:,2))
hold on;
% Plot Stations
plot(lon,lat,'*')
time = distance/velocity
t-t0 = sqrt((x-x0)^2 + (y-y0)^2)/v
% 初始猜想
lat0 = 46.9; lon0 = 9;
% 使用红色圆圈标出初始猜想震中
plot(lon0,lat0,'ro')
% 使用绿色圆圈标记地震台G
plot(lon(6),lat(6),'go')
% 将度数转换成公里的参数
latkm = 111.19; lonkm = 75.82; vp = 5.8;
% 初始计算
t0 = sqrt( ((lon(6)-lon0)*lonkm).^2 + ((lat(6)-lat0)*latkm).^2 )/vp
%初始猜想存储到数组m
m(1) = t0; m(2)= lon0; m(3) = lat0;
%循环遍历所有地震台
for i=1:length(lat)
        diffx = (lon(i)-m(2))*lonkm;
        diffy = (lat(i)-m(3))*latkm;
        D0(i) = sqrt(diffx^2+diffy^2);          % 距离
        d(i) = tim(i) - D0(i)/vp - m(1);        % 时间
        % 存储结果到矩阵G
        G(i,3) = -(diffy/D0(i))/vp; G(i,2) = -(diffx/D0(i))/vp; G(i,1) = 1;
end
G*m = d
dm = inv(G'*G)*G'*d'
% 选择迭代次数
for i = 1:6
    for i=1:length(lat)
        diffx = (lon(i)-m(2))*lonkm;
        diffy = (lat(i)-m(3))*latkm;
        D0(i) = sqrt(diffx^2+diffy^2);
        d(i) = tim(i) - D0(i)/vp - m(1);
        G(i,3) = -(diffy/D0(i))/vp; G(i,2) = -(diffx/D0(i))/vp; G(i,1) = 1;
    end
% 最小二乘解
dm = inv(G'*G)*G'*d'
% 将dm变为度数，因为坐标系中用度数
dm(2) = dm(2)/lonkm;
dm(3) = dm(3)/latkm;
% Update 'm' Vector
m = m+dm'
% 使用蓝色圆圈标识每次迭代结果
plot(m(2),m(3),'o')
end
%使用黑色菱形框线标识最终结果
plot(m(2),m(3),'kd')
dm =
1.0e-006 *
-0.0262
0.4140
0.9619
m =
40.2578 7.5580 47.2169