lsq 最小二乘估计器

使用方法：lsq config_file

在残差编辑和坐标解算以及模糊度固定时均用到了lsq估计器，其后面跟的配置文件为中间文件，即pride自行生成的用于lsq的配置文件，但内容与原始的config_template相当。

lsq.f90

此程序时lsq可执行程序的主程序。我对其浅陋的理解如下：

初始化本地全局变量dcb和LCF%prn(i)。dcb为差分码偏差；LCF为整个程序的控制信息结构体；prn为LCF结构体内保存的卫星编号。
调用函数get_lsq_args(LCF, SITE, OB, SAT)对配置文件进行读取并保存到LCF中。

# 函数作用： 获得最小二乘的参数信息
# 参数：输入： LCF   --最小二乘配置信息结构体SITE  --测站相关结构体SAT   --卫星相关信息结构体
subroutine get_lsq_args(LCF, SITE, OB, SAT)* 读取参数，倘若参数小于1，则输出lsq的帮助语句。否则，读取并打开配置文件；* 读取的内容包括：LCF%jd0     --开始简化儒略日LCF%sod0    --开始的儒略日秒LCF%jd1     --结束时间简化儒略日LCF%sod1    --结束时间的儒略日秒LCF%flnorb   --gps卫星轨道文件名LCF%flnpos   --结果位置文件名LCF%flnsck   --卫星钟差文件名LCF%flnrck   --结果接收机钟差文件名LCF%flnztd   --结果天顶对流层延迟文件名LCF%flnhtg   --结果水平对流层梯度文件名LCF%flnamb   --结果浮点模糊度文件名LCF%flnres   --结果残差文件名LCF%flncon   --结果模糊度约束文件名（整数模糊度）LCF%flnneq   --法方称矩阵逆矩阵文件名LCF%flnvmf   --对流层投影函数VMF1文件名LCF%flnfcb   --相位偏差文件名LCF%flnerp   --地球旋转参数文件名LCF%dintv   --采样间隔LCF%lrmbias --是否参数消去（true or false）LCF%ztdmod  --天顶对流层估计模型LCF%htgmod  --水平对流层梯度模型LCF%nprn    --总的卫星个数LCF%prn(LCF%nprn) --卫星编号SAT(LCF%nprn)%prn --卫星编号SAT(LCF%nprn)%sys --卫星系统(只有“G” GPS)SAT(LCF%nprn)%typ --卫星类型(只有“Block”)SITE%name   --测站名SITE%skd    --处理策略(动态K or 静态S)SITE%map    --对流层投影函数SITE%obsfil  --测站观测O文件SITE%rhdfile --处理后的观测文件SITE%lfnjump --钟跳文件SITE%x      --测站位置（xyz）SITE%geod   --测站大地坐标（blh）SITE%rot12f --enu2xyz的旋转矩阵SITE%olc    --海洋潮汐OB          --观测数据结构体* 检查观测文件是否可用

调用read_dcb(dcb)函数读取dcb。

# 函数作用:读取dcb偏差
# 函数位置：src/lib/read_dcb.f90
# 参数：输出：dcb    --dcb偏差(二维数组)
subroutine read_dcb(dcb)* 打开P1C1.dcb文件，将P1C1的dcb保存到dcb(prn0, 1)中，dcb是一个二维数组。* 打开并读取P2C2.dcb文件，将P2C2的dcb保存到dcb(prn0, 2)中

调用read_snx(fcb_file, bias)函数读取相位偏差。这里的fcb_file是武汉大学生成的相位偏差，包括相位钟和fcb相位偏差。

# 函数作用：读取小数偏差用于模糊度解算
# 函数位置：src/lsq/read_snx.f90
# 参数：输出： bias    --小数偏差（二维数组）
subroutine read_snx(flnfcb, bias)* 打开fcb文件，从'+BIAS/SOLUTION'部分以后开始读取偏差值，将偏差保存到dcb(iprn, 4)二维数组中。分别为每颗卫星发射4种GPS信号L1C、L2W，C1W和C2W的伪标定小数相位偏差。

调用get_ant_ipt(fjd_beg, fjd_end, antnam, antnum, iptatx, enu)函数，对卫星天线相位中心进行改正。

# 函数作用：
# 所在位置：../lib/get_antenna_corr.f90
# 参数：输入： fjd_beg, fjd_end  --起止时间antnam, antnum    --天线名和序列号输出： iptatx            --天线指针enum              --天线偏移
subroutine get_ant_ipt(fjd_beg, fjd_end, antnam, antnum, iptatx, enu)* 前面从内存中读取天线名、天线序号的操作也没看懂* 若内存中没有，则从atx文件中读取，调用的是rdatx(fjd_beg, fjd_end, ATX)函数, 将天线相位改正保存在ATX%enu(j,k) 其中j是三个方向上的改正值，k为频率。PCV保存在ATX%pcv()中* 将ATX%enu输出为enu(6)，其中前三个为第一个频率的，后三个为第二个频率的。* 疑问：最终只输出了enu，没有输出pcv。

打开观测值文件，读取测站天线相关信息，并同样调用get_ant_ipt()函数获取测站天线相位改正，并保存在SITE%enu 中。
调用lsq_cnt_prmt(LCF, SITE, NM)函数统计参数的个数。

# 函数作用：统计最小二乘估计器的待估参数的个数
# 所在位置：src/lsq_cnt_prmt.f90
# 参数：输入： LCF  --最小二乘配置相关信息结构体SITE --测站相关信息结构体输出： NM   --参数相关结构体
subroutine lsq_cnt_prmt(LCF, SITE, NM)* 统计的时候NM%nc为常数参数的个数；NM%np为过程参数的个数；   * 当测站处理策略时S时，即静态时，测站三个坐标当作常数估计；当为K时，即动态时，三个测站坐标当过程参数估计加入NM%np中。* 对流层参数包括天顶对流层、2个方向水平对流层梯度，当其处理策略为分段常数时，都看作过程参数，加入NM%np中。* 接收机钟差当作过程参数，加入NM%np中。* 总参数NM%imtx = 过程参数NM%np + 常数参数NM%nc

总的参数的个数即为NM%imtx加上模糊度的个数NM%amb_tot,释放一个PM(NM%npm)空间存放这些参数。其中NM%npm即为所有的待求参数个数。
正规矩阵的大小。如果需要参数消去，那么参数个数需要加OB%amb_epo+1，其中OB%amb_epo为历元内新添加的模糊度参数最大值；如果不需要参数消去，那么参数个数只需要+1。这里不太明白为什么这样！！！！,经过这一步以后，正规矩阵的大小即为NM%nmtx。
释放内存，用NM%norx(1:NM%nmtx,1:NM%nmtx)这么一个二维数组当作设计矩阵。用NM%iptp(NM%nmtx)存放；用NM%iptx(NM%nmtx)存放。
调用lsq_init(LCF, SITE, SAT, OB, NM, PM)函数初始化设计矩阵。包括所有最小二乘相关项的初始化。
开始按历元循环求解：每个历元操作步骤：
（1）读取观测值，通过RINEX版本不同调用不同的函数读取该历元匹配的RINEX观测值。注意这里读取的观测值已经进行了bias和dcb的改正，都是直接在观测值后面加的。
（2）调用read_obsrhd(jd, sod, nprn, prn, OB)函数对rhd_函数进行读取，rhd文件为观测值文件的健康诊断文件，可以确定历元内哪颗卫星观测值需要删除，哪个历元哪颗卫星需要新添加模糊度参数。
（3）如果有接收机钟跳，还需要调用read_clkjmp(jd, sod, clkjmp_file, nprn, OB)函数读取跳秒,这里的钟跳文件由tedit产生。
（4）对于每颗卫星，调用read_satclk函数读取并插值卫星钟差，这里插值采用的是线性插值：。调用everett_interp_orbit函数插值该历元时刻卫星位置和速度。采用的是拉格朗日插值。
（5）调用lsq_add_newamb(jd, sod, name, LCF, OB, NM, OM)函数添加新的模糊度。添加新模糊度以后需要对待估参数个数以及设计矩阵进行更新。
（6）开始按照测站循环。每个测站的操作步骤如下
一、添加观测方程，OB%omc即为方程左边部分o-c。
二、准备先验接收机钟差改正。
三、调用一种模型计算对流层先验改正。
四、检查高度角和方位角
五、又来一个保存先验接收机钟差改正，而且看代码显示和上面那个正好反过来，不太清楚时什么操作。
六、保存先验对流层水平梯度经度到PM%xini，包括n方向和e方向的梯度
七、调用lsq_add_obs函数将观测值添加到观测方程中。
（7）每个历元保存一个先验天顶对流层延迟。
（8）调用lsq_process函数做最小二乘。这里采用的是参数消去求解。
添加模糊度限制，即求解模糊度固定解
保存法方程
生成残差文件头部分
参数恢复，计算残差
输出模糊度文件，此模糊度为加入模糊度限制后的模糊度解。
输出法方程逆矩阵。

存在的问题：

天线改正时只输出enu，而不输出pcv？是pcv没有改正还是。。
设计矩阵的大小，为什么在flrmbias是否为真都要+1？这个1
是什么意思？
分配内存那里NM%iptp和 NM%iptx分别是什么，其长度都为法方程矩阵的维度？
先验接收机钟差是怎么保存的？
最小二乘中，find a bias to be removed 之后的操作是在干什么？

Pride_pppar 源码解析（三）相关推荐

Disruptor源码解析三 RingBuffer解析
目录系列索引前言主要内容 RingBuffer的要点源码解析系列索引 Disruptor源码解析一 Disruptor高性能之道 Disruptor源码解析二 Sequence相关类解析 D ...
OkHttp3源码解析(三)——连接池复用
OKHttp3源码解析系列 OkHttp3源码解析(一)之请求流程 OkHttp3源码解析(二)--拦截器链和缓存策略本文基于OkHttp3的3.11.0版本 implementation 'com ...
ReactiveSwift源码解析(三) Signal代码的基本实现
上篇博客我们详细的聊了ReactiveSwift源码中的Bag容器,详情请参见<ReactiveSwift源码解析之Bag容器>.本篇博客我们就来聊一下信号量,也就是Signal的的几种状 ...
并发编程与源码解析 (三)
并发编程 (三) 1 Fork/Join分解合并框架 1.1 什么是fork/join Fork/Join框架是JDK1.7提供的一个用于并行执行任务的框架,开发者可以在不去了解如Thread.R ...
前端入门之(vuex源码解析三)
上两节前端入门之(vuex源码解析二)我们把vuex的源码大概的撸了一遍,还剩下(插件.getters跟module),我们继续哈~ 插件童鞋们可以去看看vuex在各个浏览器的状态显示插件,小伙伴可以 ...
拆轮子-RxDownload2源码解析(三)
本文为博主原创文章,未经允许不得转载造轮子者:Season_zlc 轮子用法请戳作者链接 ↑ 前言本文主要讲述 RxDownload2 的多线程断点下载技术. 断点下载技术前提服务器必须支持按 ...
Tomcat源码解析三：tomcat的启动过程
Tomcat组件生命周期管理在Tomcat总体结构 (Tomcat源代码解析之二)中,我们列出了Tomcat中Server,Service,Connector,Engine,Host,Context ...
【Vue.js源码解析三】-- 模板编译和组件化
前言笔记来源:拉勾教育大前端高薪训练营阅读建议:建议通过左侧导航栏进行阅读模板编译模板编译的主要目的是将模板 (template) 转换为渲染函数 (render) <div> ...
Cesium源码解析三（metadata元数据拓展中行列号的分块规则解析）
目录 1.前言 2.layer.json中available参数意义 3.EPSG:4626切片及terrain分块原理 4.Cesium的terrain分块规则 5.自定义terrain分块规则 6 ...
AFNetworking2.0源码解析三
本篇说说安全相关的AFSecurityPolicy模块,AFSecurityPolicy用于验证HTTPS请求的证书,先来看看HTTPS的原理和证书相关的几个问题. HTTPS HTTPS连接建立过程 ...

Pride_pppar 源码解析（三）

lsq 最小二乘估计器

Pride_pppar 源码解析（三）相关推荐

最新文章

热门文章