KenanZhu/KNZ_GeoTrackLab
Archived
1
0 Fork 0
Code Issues Actions Packages Projects Releases Wiki Activity

Modify comments and output

Browse Source
This commit is contained in:
KenanZhu111 committed 2024-09-21 09:56:19 +08:00
1 parent e6e8ac6e9a
commit 22cc906c00
16 files changed
+97800 -228
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

2gpst.h
+10
View File
@@ -1,6 +1,16 @@
#ifndef _2GPST_H_
#define _2GPST_H_
/* -------------------------------------------------------------------------- */
/// @brief 将世界协调时转换为GPS周内秒
/// @param y 年
/// @param m 月
/// @param d 日
/// @param h 时
/// @param min 分
/// @param sec 秒
/// @return 返回GPS周内秒
/* -------------------------------------------------------------------------- */
double TimetoGPSsec(int y, int m, int d, double h, int min, double sec)
{
if (m > 2)
CppProperties.json
-21
View File
@@ -1,21 +0,0 @@
{
"configurations": [
{
"inheritEnvironments": [
"msvc_x64"
],
"name": "x64-Debug",
"includePath": [
"${env.INCLUDE}",
"${workspaceRoot}\\**"
],
"defines": [
"WIN32",
"_DEBUG",
"UNICODE",
"_UNICODE"
],
"intelliSenseMode": "windows-msvc-x64"
}
]
}
Draw the Track[by PyGMT]/LLA_result_for_PyGMT.txt
+31613
View File
File diff suppressed because it is too large. Load diff
Observation_Station_Site-Solving_by_Matlab/LLA_result_for_Matlab.txt
+34494
View File
File diff suppressed because it is too large. Load diff
Observation_Station_Site-Solving_by_Matlab/Position_Bias.png
BIN
View File
Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 143 KiB

Observation_Station_Site-Solving_by_Matlab/Position_Bias_TGD.png
BIN
View File
Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 145 KiB

Observation_Station_Site-Solving_by_Matlab/Station_Site_Solving.asv
-170
View File
@@ -1,170 +0,0 @@
% 打开坐标文件
file=fopen('LLA_result_for_Matlab.txt');
sPRN = [];
satX = [];
satY = [];
satZ = [];
P_t1 = [];
P_t2 = [];
H = [];
TGD = [];
Dt = [];
sat_num = [];
C_V = 299792458;
f1 = 1575.42;
f2 = 1227.60;
satnum = 0;
eponum = 0;
% 设置测站粗略坐标和精确的经纬高数据
apX = -4052052.8310;
apY = 4212835.9820;
apZ = -2545104.4760;
apB = -0.4131213;%纬度
apL = 2.3367431;%经度
apH = 603.235;
% 逐行读取文件
while ~feof(file)
line = fgets(file);
flag_epo = extractBetween(line,1 ,1 );
% 判断历元标识
if strcmp(flag_epo, ">") == 1
eponum = str2double(extractBetween(line,2 ,5 ));
while eponum <= 2880
flag_sat = extractBetween(line,1 ,1);
% 判断卫星标识，读取卫星PRN号和XYZ坐标
if strcmp(flag_sat, "G") == 1
satnum = satnum + 1;
P_t1(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16 ,25 + 16 + 16 + 16 ));
if P_t1(eponum, satnum) ~= 0
sPRN(eponum, satnum) = str2double(extractBetween(line,2 ,3 ));
satX(eponum, satnum) = str2double(extractBetween(line,11 ,25 ));
satY(eponum, satnum) = str2double(extractBetween(line,11 + 16 ,25 + 16 ));
satZ(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 ,25 + 16 + 16 ));
P_t1(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16 ,25 + 16 + 16 + 16 ));
P_t2(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16 ,25 + 16 + 16 + 16 + 16));
H(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16 + 16 ,25 + 16 + 16 + 16 + 16 + 16));
TGD(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16 + 16 + 16 ,25 + 16 + 16 + 16 + 16 + 16 + 16));
Dt(eponum, satnum) = str2double(extractBetween(line,11 + 16 + 16 + 16 + 16 + 16 + 16 + 16,25 + 16 + 16 + 16 + 16 + 16 + 16 + 17));
line = fgets(file);
end
if P_t1(eponum, satnum) == 0
satnum = satnum - 1;
line = fgets(file);
end
end
% 开头字符不是G，但是是>,判断为一个历元的起点，进行数据读取
if strcmp(flag_sat, "G") == 0 && strcmp(flag_sat, ">") == 1
satnum = 0;
line = fgets(file);
continue;
end
% 开头字符既不是G，也不是>,判断为历元结束标志，跳出循环
if strcmp(flag_sat, "G") == 0 && strcmp(flag_sat, ">") == 0 && ~feof(file)
sat_num(eponum) = satnum;
satnum = 0;
break;
end
end
end
if strcmp(flag_epo, ">") == 0
continue;
end
end
fclose(file);
% 平差计算测站位置和接收机钟差解算
for i = 1:1:eponum
if sat_num(i) >= 4
X = zeros(4, 1);
P = zeros(sat_num(i), sat_num(i));
B = zeros(sat_num(i), 4);
l = zeros(sat_num(i), 1);
X_R = zeros(3, 1);
DeltaXYZ = zeros(3, 1);
% 每个历元构成系数矩阵
for j = 1:1:sat_num(i)
R = sqrt((satX(i, j)-apX)^2+(satY(i, j)-apY)^2+(satZ(i, j)-apZ)^2);
B(j, 1) = -((satX(i, j) - apX) / R);
B(j, 2) = -((satY(i, j) - apY) / R);
B(j, 3) = -((satZ(i, j) - apZ) / R);
B(j, 4) = 1;
dtrop = 2.47/sind(H(i, j))+0.0121;
% % % % % % % 理想化的Saastamoinen模型/对流层改正
% P0 = 1013.25; e0 = 11.69; T0 = 288.15;
% % 求大气压强
% p = P0 * (1.0 - 0.0068 / T0 / apH)^5;
% % 求开尔文温度
% T = 288.15 - 0.0068 * apH;
% % 水汽压计算
% if apH > 11000
% e = 0;
% end
% if apH <= 11000
% e = e0 * (1.0 - 0.0068 * apH / T0)^4;
% end
% % 静力学延迟
% dtrop = 0.002277 * p / (1.0 - 0.00266 * cosd(2 * H(i, j)) - 0.00028 * apH * 1e-3) ;
% % % % % % % % % % % % % % % % % % % % % % % % %
if P_t2(i, j) ~= 0
% 双频组合消除电离层延迟
P_t = (f1^2) / (f1^2 - f2^2) * P_t1(i, j) - (f2^2) / (f1^2 - f2^2) * P_t2(i, j);
end
if P_t2(i, j) == 0
P_t = P_t1(i, j);
end
l(j, 1) = P_t - R + C_V * Dt(i,j) - dtrop;
if H(i, j)<5;
H(i, j) = H(i, j);
end;
P(j, j) = sind(H(i, j))^2;
end
apX_R = apX;%迭代坐标传参
apY_R = apY;
apZ_R = apZ;
% 平差计算得到改正数并进行迭代计算
while abs(apX_R - X_R(1,1)) > 0.0001 && abs(apY_R - X_R(2,1)) > 0.0001 && abs(apZ_R - X_R(3,1)) > 0.0001
apX_R = X(1,1);
apY_R = X(2,1);
apZ_R = X(3,1);
Q=inv(B.'*P*B);
D_X=Q*B.'*P*l;%坐标改正数
X(1,1) = apX + D_X(1,1);
X(2,1) = apY + D_X(2,1);
X(3,1) = apZ + D_X(3,1);%改正后测站坐标
X_R(1,1) = X(1,1);
X_R(2,1) = X(2,1);
X_R(3,1) = X(3,1);
end
GDOP(1,i) = sqrt(trace(Q));
StationCLK(1,i) = D_X(4, 1);
DeltaXYZ(1,1) = X(1,1) - apX;
DeltaXYZ(2,1) = X(2,1) - apY;
DeltaXYZ(3,1) = X(3,1) - apZ;
S = [
-sin(apL), cos(apL), 0;
-sin(apB)*cos(apL), -sin(apB)*sin(apL), cos(apB);
cos(apB)*cos(apL), cos(apB)*sin(apL), sin(apB)
];
ENU_R = S * DeltaXYZ;
ENU(1,i) = ENU_R(1,1);
ENU(2,i) = ENU_R(2,1);
ENU(3,i) = ENU_R(3,1);
end
end
Observation_Station_Site-Solving_by_Matlab/Station_Site_Solving.m
+13 -2
View File
@@ -123,9 +123,9 @@ for i = 1:1:eponum
P_t = P_t1(i, j);
end
l(j, 1) = P_t - R + C_V * Dt(i,j) - dtrop;
if H(i, j)<5;
if H(i, j)<5
H(i, j) = H(i, j);
end;
end
P(j, j) = sind(H(i, j))^2;
end
@@ -169,3 +169,14 @@ for i = 1:1:eponum
end
end
%绘制结果图
figure("Position",[200, 700,1000,200],"Name", "East component deviation","NumberTitle","off")
bar(ENU(1,:))%E
box off
figure("Position",[200, 400,1000,200],"Name","North component deviation","NumberTitle","off")
bar(ENU(2,:))%N
box off
figure("Position",[200, 100,1000,200],"Name", "Up component deviation","NumberTitle","off")
bar(ENU(3,:))%U
box off
Pos_out/LLA_result_for_read.txt
+31613
View File
File diff suppressed because it is too large. Load diff
README.md
+4 -4
View File
@@ -19,10 +19,10 @@
* Due to my limited programming level, PyGMT drawing used too miscellaneous judgment statements, resulting in slow drawing or debugging errors, please forgive me.
* When writing the code, with the help of the existing open source code results, and different degrees of reference, here one by one mark and express thanks:
| https://blog.csdn.net/why1472587?type=blog
| https://zhuanlan.zhihu.com/p/416072448
| https://www.pygmt.org/latest/index.html
| https://blog.csdn.net/FrankXCR/article/details/135438701
| *https://blog.csdn.net/why1472587?type=blog*
| *https://zhuanlan.zhihu.com/p/416072448*
| *https://www.pygmt.org/latest/index.html*
| *https://blog.csdn.net/FrankXCR/article/details/135438701*
cal_mod.h
+11 -8
View File
@@ -9,6 +9,7 @@
#include"select_epo.h"
#include"xyz2blh.h"
/* -------------------------------------------------------------------------- */
#define a 6378137.0//长半轴
#define f (1 / 298.257223563)//扁率
#define e2 (f*(2-f))//第一偏心率平方
@@ -18,6 +19,7 @@
#define PI 3.141592653589793
#define Earth_e 7.2921151467e-5 //地球自转角速度
#define C1 0
/* -------------------------------------------------------------------------- */
//卫星位置结构体
typedef struct
@@ -69,7 +71,7 @@
double H;
}rah, *prah;
//GPS卫星位置计算函数
/* ------------------------------- GPS卫星位置计算函数 ------------------------------ */
int gps_pos(int eponum, int sPRN, int best_epoch, double GPSsec,
pnav_body nav_b, pobs_head obs_h, ppos_t pos_t)
{
@@ -114,8 +116,10 @@
double rela = 2*sqrt(GM) * nav_b[best_epoch].e * nav_b[best_epoch].sqrtA * sin(Es)/pow(C_V,2);
pos_t[eponum].delta_clk[sPRN] = nav_b[best_epoch].sa0 + nav_b[best_epoch].sa1 * tk + nav_b[best_epoch].sa2 * pow(tk,2) - rela;
}
/* -------------------------------------------------------------------------- */
//GPS卫星位置计算
/* ------------------------------- GPS卫星位置计算 ------------------------------ */
int sat_gps_pos_clac(FILE * result_file,
pnav_body nav_b, pobs_epoch obs_e, pobs_body obs_b,
pobs_head obs_h, pstation station, ppos_t pos_t,
@@ -157,7 +161,7 @@
gps_pos( i, sPRN, best_epoch, GPSsec, nav_b, obs_h, pos_t);
}
//地固坐标系转经纬大地高坐标系
/* ----------------------------- 地固坐标系转经纬大地高坐标系 ----------------------------- */
pxyz2blh tem1 = NULL;
tem1 = (pxyz2blh)malloc(sizeof(pxyz2blh));
tem1 = XYZtoBLH( tem1, pos_t[i].X[sPRN], pos_t[i].Y[sPRN], pos_t[i].Z[sPRN], a, e2);
@@ -172,8 +176,7 @@
station->L = tem2->L / 57.295779513;
station->H = tem2->H;
//大地高坐标系转东北天坐标系*******************************************************//
/* -------------------------------------------------------------------------- */
double sinL = sin(station->L);
double cosL = cos(station->L);
double sinB = sin(station->B);
@@ -181,8 +184,8 @@
enu->E = -sinL*(pos_t[i].X[sPRN] - obs_h->apX) + cosL*(pos_t[i].Y[sPRN] - obs_h->apY);
enu->N = -sinB*cosL*(pos_t[i].X[sPRN] - obs_h->apX) - sinB*sinL*(pos_t[i].Y[sPRN] - obs_h->apY) + cosB*(pos_t[i].Z[sPRN] - obs_h->apZ);
enu->U = cosB*cosL*(pos_t[i].X[sPRN] - obs_h->apX) + cosB*sinL*(pos_t[i].Y[sPRN] - obs_h->apY) + sinB*(pos_t[i].Z[sPRN] - obs_h->apZ);
//*******************************************************************************//
//卫星方位角a，高度角h，向径r计算
/* -------------------------------------------------------------------------- */
/* ---------------------------- 卫星方位角a，高度角h，向径r计算 --------------------------- */
rah->H = atan2(enu->U, sqrt(enu->E * enu->E + enu->N * enu->N)) * 57.295779513;
rah->A = atan2(enu->E, enu->U) * 57.295779513;
if (rah->A < 0)
@@ -190,7 +193,7 @@
if (rah->A > 2 * PI)
rah->A -= 2 * PI;
rah->R = sqrt(enu->E * enu->E + enu->N * enu->N + enu->U * enu->U);
//*******************************************************************************//
/* -------------------------------------------------------------------------- */
if (blh->H < 0 || nav_b[best_epoch].sHEA != 0)
{
main.c
+3 -2
View File
@@ -1,5 +1,6 @@
//仅支持{ Ver 3.04 }版本RENIX观测值文件与广播星历文件读取
/* -------------------------------------------------------------------------- */
/* 仅支持{ Ver 3.04 }版本RENIX观测值文件与广播星历文件读取
/* -------------------------------------------------------------------------- */
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
main.exe
BIN
View File
Binary file not shown.
read.h
+20 -15
View File
@@ -22,7 +22,7 @@ void setstr(char* des, const char* src, int n)
}
}
//获取文件数据块行数，从END OF HEADER后开始起算
/* ---------------------- 获取文件数据块行数，从END OF HEADER后开始起算 --------------------- */
extern int getgpssatnum(FILE* fp_nav)
{
int gps_satnum = 0;
@@ -60,7 +60,7 @@ extern int getgpssatnum(FILE* fp_nav)
return gps_satnum;
}
//将字符串转换为浮点数,i起始位置，n输入多少个字符
/* ------------------------ 将字符串转换为浮点数,i起始位置，n输入多少个字符 ----------------------- */
static double strtonum(const char* buff, int i, int n)
{
double value = 0.0;
@@ -86,7 +86,7 @@ static double strtonum(const char* buff, int i, int n)
return sscanf(str, "%lf", &value) == 1 ? value : 0.0;
}
//导航头文件结构体
/* -------------------------------- 导航头文件结构体 -------------------------------- */
typedef struct nav_head
{
double ver;//rinex 版本号
@@ -109,7 +109,7 @@ typedef struct nav_head
int leap;
}nav_head, *pnav_head;
//导航数据结构体
/* --------------------------------- 导航数据结构体 -------------------------------- */
typedef struct nav_body
{
//数据块第一行内容：
@@ -169,7 +169,7 @@ typedef struct nav_body
}nav_body, *pnav_body;
//观测值头文件
/* --------------------------------- 观测值头文件结构体 --------------------------------- */
typedef struct obs_head
{
double ver;//RINEX文件版本号
@@ -198,7 +198,7 @@ typedef struct obs_head
}obs_head,*pobs_head;
//观测值历元数据结构体
/* ------------------------------- 观测值历元数据结构体 ------------------------------- */
typedef struct obs_epoch
{
//观测历元时刻
@@ -218,16 +218,19 @@ typedef struct obs_epoch
}obs_epoch,*pobs_epoch;
//观测值数据结构体
/* -------------------------------- 观测值数据结构体 -------------------------------- */
typedef struct obs_body
{
double obs_gps[100][24];//GPS观测值
double obs_bds[100][24];//BDS观测值
}obs_body,*pobs_body;
/* -------------------------------------------------------------------------- */
//读取O文件******************************************************************************************************//
//获取O文件历元数
/* -------------------------------------------------------------------------- */
/* 读取O文件
/* -------------------------------------------------------------------------- */
/* -------------------------------- 获取O文件历元数 -------------------------------- */
extern int get_epochnum(FILE* fp_obs)
{
int n = 0;//记录历元数
@@ -246,7 +249,7 @@ extern int get_epochnum(FILE* fp_obs)
return n;
}
//读取O文件数据头
/* -------------------------------- 读取O文件数据头 -------------------------------- */
extern void read_o_h(FILE* fp_obs, pobs_head obs_h)
{
char buff[MAXRINEX] = { 0 };
@@ -357,7 +360,7 @@ extern void read_o_h(FILE* fp_obs, pobs_head obs_h)
}
}
//读取O文件数据块
/* -------------------------------- 读取O文件数据块 -------------------------------- */
extern void read_o_b(FILE* fp_obs, pobs_epoch obs_e, pobs_body obs_b, int type_gps, int type_bds)
{
int n = 0;//历元数
@@ -431,8 +434,10 @@ extern void read_o_b(FILE* fp_obs, pobs_epoch obs_e, pobs_body obs_b, int type_g
}
//读取N文件***************************************************************************************************************//
//读取N文件数据头
/* -------------------------------------------------------------------------- */
/* 读取N文件
/* -------------------------------------------------------------------------- */
/* -------------------------------- 读取N文件数据头 -------------------------------- */
void read_n_h(FILE* fp_nav, pnav_head nav_h)
{
char buff[MAXRINEX] = { 0 };
@@ -505,7 +510,7 @@ void read_n_h(FILE* fp_nav, pnav_head nav_h)
}
}
//读取N文件数据块
/* -------------------------------- 读取N文件数据块 -------------------------------- */
void read_n_b(FILE* fp_nav, pnav_body nav_b)
{
int i_g = 0;//第i_g个GPS卫星数据
@@ -596,6 +601,6 @@ void read_n_b(FILE* fp_nav, pnav_body nav_b)
}
}
}
/* -------------------------------------------------------------------------- */
#endif
select_epo.h
+9 -3
View File
@@ -1,16 +1,22 @@
#ifndef _SELECT_EPO_H_
#define _SELECT_EPO_H_
#include<stdio.h>
#include<math.h>
#include"read.h"
//挑选合适的星历
/* -------------------------------------------------------------------------- */
/// @brief 寻找与观测历元最近的卫星广播历元
/// @param SecofWeek 计算得到的观测历元的周内秒
/// @param sPRN 需要匹配的该历元卫星号
/// @param nav_b 导航文件数据块结构体，与卫星号匹配
/// @param gps_satnum 该历元的GPS卫星数量
/// @return 返回最佳历元的绝对历元序号
/* -------------------------------------------------------------------------- */
extern int select_epoch(double SecofWeek, int sPRN, pnav_body nav_b, int gps_satnum)
{
int best_epoch = 0;
double min=10000;//假设最小值是1万秒
double min = 10000;//初始化最小值
double Min;
int i;
for (i = 0; i < gps_satnum; i++)
xyz2blh.h
+10 -3
View File
@@ -1,7 +1,6 @@
#ifndef _XYZ2BLH_H_
#define _XYZ2BLH_H_
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
@@ -14,6 +13,16 @@ typedef struct
double H;
}xyz2blh, *pxyz2blh;
/* -------------------------------------------------------------------------- */
/// @brief 将ECEF地心地固坐标转换为经纬高坐标
/// @param xyz2blh 自定义结构体类型，用来传参
/// @param X ECEF的X坐标
/// @param Y ECEF的Y坐标
/// @param Z ECEF的Z坐标
/// @param a 椭球的长半轴
/// @param e2 椭球的第一偏心率平方
/// @return 返回自定义结构体类型包含B，L，H成员
/* -------------------------------------------------------------------------- */
pxyz2blh XYZtoBLH(pxyz2blh xyz2blh,
double X, double Y, double Z, double a, double e2)
@@ -42,6 +51,4 @@ pxyz2blh XYZtoBLH(pxyz2blh xyz2blh,
return xyz2blh;
}
#endif