希望有人能够帮我看一下这个matlab水准网抗差估计编程的错误之处。顺便改正一下,感激不尽!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%水准网抗差估计
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%闽江学院测绘工程专业《变形监测数据处理》课程专用
%程序功能:读取水准网观测数据,进行抗差估计
%         程序中与普通水准网平差程序不同的地方用############标示出来了
%
%观测数据示例:
% Ha,11               %已知点高程
% Hb,11.5             %已知点高程
% Hc,12.008           %已知点高程
% Ha,P1,1.003,1       %观测值:后视点名,前视点名,高差(m),距离(km)
% P1,P2,0.501,2       %观测值:后视点名,前视点名,高差(m),距离(km)
% Hc,P2,0.503,2       %观测值:后视点名,前视点名,高差(m),距离(km)
% Hb,P1,0.505,1       %观测值:后视点名,前视点名,高差(m),距离(km)
%
%计算结果:Hp1=12.0047    Hp2=12.5083
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%清空已有的变量,清空命令行显示的内容
clear;
clc;%初始化一些数值变量,点名的长度未知,可以在读取文件的时候直接赋值
iObsNum=0;   %观测值的个数
dObsH=[];    %观测高差
dObsL=[];    %测段的距离
iKnown=0;    %已知点的个数
dKnown=[];   %每个已知点的高程%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%读取观测文件
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%观测数据中带有逗号,有字符(点名),又有数字(高差,高程,距离),不能直接用laod命令
fId=fopen('水准观测数据2.txt');  
while 1  %循环,直接读取文件结束为止
   line = fgetl(fId);           %读取观测文件的一行
   if line ==0, break;  end   %如果没有读取到任何内容,说明文件结束,跳出循环
   
   index=[];   %用来保存逗号的位置
   %查找读取的一行数据中的逗号,可能有一个,也可能有三个逗号。这里没考虑更复杂的情况,cosa软件的水准网数据格式更复杂,
   %不仅包含水准测段距离,还包含测站数等内容,具体请查看cosa说明书
   index = findstr(line,',');  %返回逗号的位置。(这里没有考虑中文逗号的问题,实际工作中可能会有中文、英文逗号混合的情况,请同学们自行完善程序)
   if size(index)==1           %逗号的个数只有一个,说明是已知点
       iKnown=iKnown+1;        %统计已知点个数
       temp='';                %
       temp=deblank(line(1:index(1)-1));                    %得到已知点点名。deblank函数去掉点名前后的空格
       strKnown(iKnown,1:length(temp) )=temp;
       dKnown(iKnown)=str2num(line(index(1)+1:end) );       %得到已知点的高程。str2num函数表示将字符转换成数字
   elseif size(index,2)>1      %逗号多于一个,说明是观测高差
       iObsNum=iObsNum+1;      %统计观测值个数
       temp='';          
       temp=deblank(line(1:index(1)-1));                    %后视点点名strHou
       strHou(iObsNum,1:length(temp))=temp;
       temp='';
       temp=deblank(line(index(1)+1:index(2)-1));           %前视点点名strQian
       strQian(iObsNum,1:length(temp))=temp;
       
       dObsH(iObsNum)=str2num(line(index(2)+1:index(3)-1)); %高差
       dObsL(iObsNum)=str2num( line(index(3)+1:end) );      %距离
   end;
   
end;
fclose(fId);  %读取结束后,关闭观测文件%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%搜索未知点
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%在平差之前先要知道未知高程点的个数,具体点名
%搜索方法是从所有观测数据中找出所有的点名(包括已知点),然后再把已知点去掉
iParam=iKnown;            %iParam表示未知高程点的个数。把已知点加进去
strParam=strKnown;        %未知点的点名
%循环所有的观测值,先从后视点名中逐个查找
for i=1:iObsNum
    inew=0;               %表示所查找的点名不在strParam()数组中,那么增加到strParam中去
    for j=1:iParam
        if strncmp(strParam(j,:),strHou(i,:),size(strHou(i,:),2) )  %相同
            inew=0;       %当前的后视点点名已经在strParam中
            break;
        else
            inew=i;       %
        end
    end
    %inew不等于0,那么将inew所对应的那个后视点加入到strParam()中去
    if inew               
       iParam=iParam+1;   %未知点个数增1
       strParam(iParam,1:length(strHou(inew,:)))=strHou(inew,:);        
    end;
end;
%查找前视点,与前面的一样的道理
for i=1:iObsNum
    inew=0;
    for j=1:iParam
        if strncmp(strParam(j,:),strQian(i,:),size(strQian(i,:),2) )  %相同
            inew=0;
            break;
        else
            inew=i;
        end
    end
    if inew
       iParam=iParam+1;
       strParam(iParam,1:length(strQian(inew,:)))=strQian(inew,:);        
    end;
end;
%经过循环,已经找到了观测文件中所有的点名,而且不重复,包含了已知点点名(放在最前面的几个)
%直接把已知点去掉,剩下的就是未知点
iParam=iParam-iKnown;               %
strParam=strParam(iKnown+1:end,:)   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%生成B  L矩阵   (P矩阵要在循环过程中反复计算)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dB=zeros(iObsNum,iParam);
dP=zeros(iObsNum,iObsNum);
dL=zeros(iObsNum,1);
for i=1:iObsNum
    dP(i,i)=1.0/dObsL(i);  %初次计算的时候,权按照默认的与距离成反比计算
    %生成系数阵
    for j=1:iParam
        if strncmp(strParam(j,:),strHou(i,:),size(strParam,2) )  %HA   hAB=HB-HA
            dB(i,j)=-1.0;   %根据hAB=HB-HA,后视A的系数应该是-1
            break;
        end;
    end;
    for j=1:iParam
        if strncmp(strParam(j,:),strQian(i,:),size(strParam,2) ) %HB  hAB=HB-HA
            dB(i,j)=1.0;    %根据hAB=HB-HA,前视B的系数应该是1
            break;
        end;
    end;
    %常数项包含两种情况,如果前视、后视都是未知点,常数项则为观测高差,如果其中某一个是已知点,则是:观测高差 +/- 已知高程
    dL(i,1)=dObsH(i);  %
    for j=1:iKnown
        if strncmp(strKnown(j,:),strHou(i,:),size(strHou,2) ) %已知点
            dL(i,1)=dL(i,1)+dKnown(j);   %观测高差+已知高程
            break;
        end;
    end
    for j=1:iKnown
        if strncmp(strKnown(j,:),strQian(i,:),size(strQian,2) ) %已知点
            dL(i,1)=dL(i,1)-dKnown(j);   %观测高差-已知高程
            break;
        end;
    end
end
%############################################################
%####求解x,需要用到选权迭代法                                  
%####选权迭代法的权函数有各种不同形式,这里只介绍一次范数最小法   
%####其余方法请参考《实用测量数据处理方法》 P56                 
%############################################################
dx=zeros(iParam,1);
dx1=dx;
while 1    
    dx=inv(dB'*dP*dB)*(dB'*dP*dL)   %平差,得到各个未知高程点的高程    
    dXv=dx1-dx;                     %两次平差后的未知参数的变化量
    if(max(abs(dXv))<0.0001)         %如果两次平差计算结果的变化量小于1mm,终止循环
        break;
    end;
    dv=dB*dx-dL;                    %求解改正数v=Bx-L  
    dk=min(abs(dv))/100.0;          %k相对于v是一个微小量,用于避免v=0出现的定权问题
    for i=1:iObsNum
        dP(i,i)=1.0/(abs(dv(i))+dk);%重新定权
    end;
    dx1=dx;                         %保存上一次计算结果,用于跟下一次的平差结果进行比较
end;