Pride_pppar 源码解析（三）

lsq 最小二乘估计器

使用方法：lsq config_file

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

lsq.f90

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

1. 初始化本地全局变量dcb和LCF%prn(i)。dcb为差分码偏差；LCF为整个程序的控制信息结构体；prn为LCF结构体内保存的卫星编号。
2. 调用函数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          --观测数据结构体
	* 检查观测文件是否可用

2. 调用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)中

3. 调用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的伪标定小数相位偏差。

4. 调用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。

5. 打开观测值文件，读取测站天线相关信息，并同样调用get_ant_ipt()函数获取测站天线相位改正，并保存在SITE%enu 中。
6. 调用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

7. 总的参数的个数即为NM%imtx加上模糊度的个数NM%amb_tot,释放一个PM(NM%npm)空间存放这些参数。其中NM%npm即为所有的待求参数个数。
8. 正规矩阵的大小。如果需要参数消去，那么参数个数需要加OB%amb_epo+1，其中OB%amb_epo为历元内新添加的模糊度参数最大值；如果不需要参数消去，那么参数个数只需要+1。这里不太明白为什么这样！！！！,经过这一步以后，正规矩阵的大小即为NM%nmtx。
9. 释放内存，用NM%norx(1:NM%nmtx,1:NM%nmtx)这么一个二维数组当作设计矩阵。用NM%iptp(NM%nmtx)存放；用NM%iptx(NM%nmtx)存放。
10. 调用lsq_init(LCF, SITE, SAT, OB, NM, PM)函数初始化设计矩阵。 包括所有最小二乘相关项的初始化。
11. 开始按历元循环求解：每个历元操作步骤：
    （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函数做最小二乘。这里采用的是参数消去求解。
12. 添加模糊度限制，即求解模糊度固定解
13. 保存法方程
14. 生成残差文件头部分
15. 参数恢复，计算残差
16. 输出模糊度文件，此模糊度为加入模糊度限制后的模糊度解。
17. 输出法方程逆矩阵。

存在的问题：

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