function [pos_station,Deltat_receiver,GXSV,GYSV,GZSV,RXSV,RYSV,RZSV,CXSV,CYSV,CZSV,PDOP,PDOP1,PDOP2,avgP1,avgP2,number,GDOP]=...
singlePointPositionUsingOneEpoch_gps_bds(GpObsHeader,RpObsHeader,CpObsHeader,pObsText_OneEpoch,GpObsText_OneEpoch,RpObsText_OneEpoch,CpObsText_OneEpoch, GpSvHeader,RpSvHeader,CpSvHeader,...
GpSvText,RpSvText,CpSvText, GnumSV,RnumSV,CnumSV,k)
%单历元组合单点定位
%%% X Y Z为解算出的接收机天线相位中心坐标,Deltat_receiver为接收机钟差改正
format long;
Vc=299792458; %光速
X0 = GpObsHeader.X0; Y0= GpObsHeader.Y0; Z0 = GpObsHeader.Z0; %2.测站的近似坐标
GXSV=zeros(1,32);GYSV=zeros(1,32);GZSV=zeros(1,32); %给不同系统卫星位置赋空矩阵
RXSV=zeros(1,24);RYSV=zeros(1,24);RZSV=zeros(1,24);
CXSV=zeros(1,14);CYSV=zeros(1,14);CZSV=zeros(1,14);
if pObsText_OneEpoch.svNumOfEpoch >6
fault=[];prn=[];fault_num=0;
X1 = X0 + 0.1; Y1 = Y0 +0.1 ; Z1 = Z0 +0.1; %3.迭代循环条件准备
%3.1迭代循环计算测站的位置
number=0;
while abs(X0-X1)>0.01||abs(Y0-Y1)>0.01||abs(Z0-Z1)>0.01
number=number+1;
clear A;clear L;clear P;clear V;
X0 = X1;
Y0 = Y1;
Z0 = Z1;
flag=0;%记录参与解算的卫星颗数
Gnum_useful = 0;%可使用的GPS观测值
Cnum_useful = 0;%BDS
for i = 1:pObsText_OneEpoch.svNumOfEpoch
for j=1:length(fault)
if i==fault(j);
prn=0;
else
prn=[];
end
end
if prn==0
continue;
end
% if pObsText_OneEpoch.data(i,pObsHeader.typeObs.C1C)<10000||pObsText_OneEpoch.LLI(i,pObsHeader.typeObs.C1C)==0
% continue; %%%说明伪距观测值不正确
% end
if pObsText_OneEpoch.sys(i) == 'G'; %如果该颗卫星是GPS卫星
%1.常量定义
%GPS与GLONASS/BDS地球引力常数和地球自转角速度不一样,后两者一样。(根据最新接口控制文件和国际天文联合会网站可知)
mu=3.986005e+014; % mu=GM 地球引力常数
a=26560000; %卫星轨道近似圆半径a,用来计算相对论效应误差
prn = GpObsText_OneEpoch.svPrnOfEpoch(i); %卫星号
deltat_GPS = GpObsText_OneEpoch.Rdeltat; %接收机钟差
[X_sv,Y_sv,Z_sv,ts,ineartoe,ek]=G_GetSVPos_RepairedByEarthrotation_UsingReceiverPersuade(GpObsText_OneEpoch,GpSvText,GnumSV,prn);%计算第i颗卫星在该历元的位置
GXSV(prn)=X_sv;
GYSV(prn)=Y_sv;
GZSV(prn)=Z_sv;
%判断是否成功计算卫星prn位置;
if ineartoe ~= 0 %成功计算了卫星位置,并开始构建A.L矩阵的第i行向量
flag=flag+1;
Gnum_useful = Gnum_useful+1;
%计算卫星钟差改正,间p81 4.3节钟误差公式4-20
delta_t_sat = GpSvText(ineartoe).af0 + GpSvText(ineartoe).af1*(ts - GpSvText(ineartoe).toc) +...
GpSvText(ineartoe).af2*(ts - GpSvText(ineartoe).toc)^2 ;
%计算相对论效应改正,间p79 4.2 4-14
delta_rel =(-2*sqrt(mu)/Vc^2)* GpSvText(ineartoe).e * sqrt (a) * sin(ek);%%sqrt (pSvText(ineartoe).a)
%计算卫星硬件延迟改正,间p82 4.3 4-24
delta_tgd =GpSvText(ineartoe).tgd;
%计算卫星总误差改正,间p82 4.3 4-24
delta_t_sv =-Vc *(delta_t_sat + delta_rel- delta_tgd);%卫星钟差总延迟改正
[N,E,U,R,r_Polar,alpha_Sv,Elevation,phi_P,lamda_P,H]=XYZ_to_NEU(X_sv,Y_sv,Z_sv,X0,Y0,Z0,pObsText_OneEpoch.sys(i));
Vion =-Vc*Delay_of_ionosphere( Elevation, alpha_Sv, phi_P, lamda_P, GpSvHeader,GpObsText_OneEpoch.Hour);%Delay_of_ionosphere( Elevation, alpha_Sv, phi_P, lamda_P, pSvHeader,Hour);%计算电离层延迟
% Elevation为卫星的高度角,alpha_Sv为卫星方位角,phi_P和lamda_P为测站P的地心经纬度,输入单位为弧度,下面使用,要将其转换成度
% Hour为观测历元的小时
% 电离层延迟改正公式GPS与GLONASS一样,BDS不同
Vtrop = -Delay_of_troposphere( H, Elevation,25, 1);%Delay_of_troposphere( H, Elevation, model_type)
%-------------------------对流层改正模型三者一样
% H为测站的大地高
% Elevation为卫星的高度角,弧度
% model_type为模型类别,1为萨斯塔莫宁Saastamoinen模型,2为简化霍普菲尔德Hopfield模型
%卫星与测站间的几何距离理论值
p0_i = sqrt((X_sv-X0)^2 +(Y_sv-Y0)^2 + (Z_sv-Z0)^2 );
%常数项
%%使用C1C伪距计算,即L1载波上的C/A码
w0_i =p0_i-GpObsText_OneEpoch.data(i,GpObsHeader.typeObs.C1C)- Vion - Vtrop + delta_t_sv;
%%%%%%%线性化系数
l = (X_sv -X0)/p0_i; m = (Y_sv -Y0)/p0_i; n = (Z_sv -Z0)/p0_i;
A(flag,1)=-l; A(flag,2)=-m; A(flag,3)=-n; A(flag,4)=-1; A(flag,5 ) = 0; %A矩阵,有5个未知数,需要5列
L(flag,1)=-w0_i;%常数项
% %%使用C2P伪距计算,即L2载波上的P码
% w0_i2 = p0_i - pObsText_OneEpoch.data(i,pObsHeader.typeObs.C2P) - Vion - Vtrop + delta_t_sv;
% %%线性化系数
% l2=(X_sv - X0)/p0_i; m2=(Y_sv - Y0)/p0_i; n2= (Z_sv - Z0)/p0_i;
% A2(flag,1)= -l2; A2(flag,2)= -m2; A2(flag,3)= -n2; A2(flag,4)= -1; L2(flag,1)= -w0_i2;
end %if ineartoe ~= 0
end
%%%%%%%%%%%%%%%%%%%%GLONASS卫星系统%%%%%%%%%%%%%%%%%%%%%%%%%%%
if pObsText_OneEpoch.sys(i) == 'C'; %如果该颗卫星是BDS卫星
mu=3.986004418e+014; % mu=GM 地球引力常数
a=26560000; %卫星近似圆轨道半径,计算相对论效应误差
deltat_BDS = CpObsText_OneEpoch.Rdeltat; %接收机钟差
j=i-GpObsText_OneEpoch.svNumOfEpoch-RpObsText_OneEpoch.svNumOfEpoch;
prn = CpObsText_OneEpoch.svPrnOfEpoch(j); %卫星号
[X_sv,Y_sv,Z_sv,ts,ineartoe,ek]=C_GetSVPos_RepairedByEarthrotation_UsingReceiverPersuade(CpObsText_OneEpoch,CpSvText,CnumSV,prn);%计算第i颗卫星在该历元的位置
CXSV(prn)=X_sv;
CYSV(prn)=Y_sv;
CZSV(prn)=Z_sv;
if ineartoe ~= 0 %成功计算了卫星位置,并开始构建A.L矩阵的第i行向量
flag=flag+1;
Cnum_useful = Cnum_useful+1;
%计算卫星钟差改正,间p81 4.3节钟误差公式4-20
delta_t_sat = CpSvText(ineartoe).af0 + CpSvText(ineartoe).af1*(ts - CpSvText(ineartoe).toc) +...
CpSvText(ineartoe).af2*(ts - CpSvText(ineartoe).toc)^2 ;
%计算相对论效应改正,见p79 4.2 4-14
delta_rel =(-2*sqrt(mu)/(Vc^2))* CpSvText(ineartoe).e * sqrt (a) * sin(ek);%%sqrt (pSvText(ineartoe).a)
%计算卫星硬件延迟改正,间p82 4.3 4-24
delta_tgd =CpSvText(ineartoe).tgd1;
%计算卫星总误差改正,间p82 4.3 4-24
delta_t_sv =-Vc *(delta_t_sat + delta_rel- delta_tgd);%卫星钟差总延迟改正
[N,E,U,R,r_Polar,alpha_Sv,Elevation,phi_P,lamda_P,H]=XYZ_to_NEU(X_sv,Y_sv,Z_sv,X0,Y0,Z0,pObsText_OneEpoch.sys(i));
Vion =-Vc*C_Delay_of_ionosphere( Elevation, alpha_Sv, phi_P, lamda_P, CpSvHeader,CpObsText_OneEpoch.Hour);%Delay_of_ionosphere( Elevation, alpha_Sv, phi_P, lamda_P, pSvHeader,Hour);%计算电离层延迟
% Elevation为卫星的高度角,alpha_Sv为卫
评论2