CN105954764A - 一种基于椭球的gnss电离层层析投影矩阵获取方法 - Google Patents

Abstract

本发明公开了一种基于椭球的GNSS电离层层析投影矩阵获取方法，包括以下的步骤：S1：确定电离层层析待反演区域及待反演时段；S2：将待反演区域离散化为一系列三维像素；S3：准备电离层层析待反演区域及待反演时段内的GNSS卫星和GNSS接收机的大地空间直角坐标；S4：计算待反演时段内GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点；S5：计算三维像素内GNSS卫星信号的传播路径的截距长度；S6：构成投影矩阵。本发明通过基于椭球的GNSS电离层层析获取投影矩阵，克服了目前基于圆球提取投影矩阵方法时将地球椭球近似圆球的问题，保证了相对于地球椭球的投影矩阵能真实反映出GNSS信号传播路径在反演像素内的真实截距。

Description

一种基于椭球的GNSS电离层层析投影矩阵获取方法

技术领域

本发明涉及电离层层析投影矩阵获取方法，特别是涉及一种基于椭球的GNSS电离层层析投影矩阵获取方法。

背景技术

在高度约为60～1000km范围内的大气层中，由于太阳的紫外线、高能粒子和X射线等作用下，形成的存在着大量自由电子和正离子的区域，这个区域就是电离层。当电磁波在电离层中传播时会发生反射、折射、散射和吸收导致信号的路径产生弯曲、传播速度发生变化等现象。为了解释这些现象的产生需要对电离层的结构进行探测，电离层探测主要是为了获取表征电离层特征的参数(主要研究电子密度)。

随着全球卫星导航系统(Global Navigation Satellite System，GNSS)技术的发展，为电离层的研究带来新的技术和方法。利用GNSS观测数据对电离层探测，主要是获得沿卫星信号传播路径上的总电子含量(total electron content,TEC)，通常TEC被认为是投影函数与相应垂直总电子含量(Vertical Total Electronic Content，VTEC)的函数。电离层TEC的空间上的分布只能反映出电离层某一层的水平结构，却难以描述其垂直结构。层次技术(Computerized Tomography,CT)在电离层探测中的应用形成了电离层层析成像技术(Computerized Ionospheric Tomography,CIT)，该技术对电离层电子密度的时空变化进行分层，克服了单层模型的局限性，且适合于大尺度的电离层电子密度垂直分布及其扰动监测，实现了对电离层垂直结构的研究。目前，利用GNSS的电离层层析成像技术取得了很大的发展，但由于GNSS的电离层层析复杂性，电离层层析还有很多问题需要进一步的研究。

发明内容

发明目的：本发明的目的是提供一种有利于分析高分辨率电离层时空变化的基于椭球的GNSS电离层层析投影矩阵获取方法。

技术方案：为达到此目的，本发明采用以下技术方案：

本发明所述的基于椭球的GNSS电离层层析投影矩阵获取方法，包括以下的步骤：

S1：确定电离层层析待反演区域及待反演时段；

S2：将待反演区域离散化为一系列三维像素；

S3：准备电离层层析待反演区域及待反演时段内的GNSS卫星和GNSS接收机的大地空间直角坐标；

S4：计算待反演时段内GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点：首先确定折半搜索算法在三维像素各个面的目标位置，然后计算GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点；

S5：计算三维像素内GNSS卫星信号的传播路径的截距长度；

S6：构成投影矩阵。

进一步，所述步骤S1包括以下的步骤：

S1.1：确定电离层层析待反演区域；

确定电离层层析待反演区域的经度范围[L₁ L_n]、纬度范围[B₁ B_m]和椭球高范围[H₁ H_q]；

S1.2：确定电离层层析待反演时段[T₁ T_t]，电离层层析待反演的时刻为[T₁ T₂... T_t]。

进一步，所述步骤S2包括以下的步骤：

S2.1：根据电离层层析待反演区域的经度范围、纬度范围和椭球高范围，将待反演区域离散化为一系列的三维像素，即：

经度被离散化为[L₁ L₂ ... L_n]；

纬度被离散化为[B₁ B₂ ... B_m]；

椭球高被离散化为[H₁ H₂ ... H_q]；

S2.2：对三维像素按照经度从西到东、纬度从南到北和椭球高从下到上的方式进行编号：当三维像素的左下角的经度为L_i、纬度为B_j和椭球高为H_k，根据式(1)确定这个三维像素编号f：

f＝(k-1)×((m-1)×(n-1))+(j-1)×(n-1)+i (1)

其中i＝1,2,…,n、j＝1,2,…,m、k＝1,2,…,q。

进一步，所述步骤S4中确定折半搜索算法在三维像素各个面的目标位置包括以下的步骤：

S4.1：确定折半搜索算法在三维像素椭球高面的目标位置：

折半搜索算法在三维像素椭球高面的目标位置分别为S2.1中被间隔的椭球高[H₁H₂ ... H_q]；

S4.2：判断待反演时段内GNSS卫星到GNSS接收机的信号传播路径是否属于电离层层析待反演区域；

根据折半搜索算法计算出T_s时刻GNSS信号传播路径与椭球高面H₁、H_q的交点的大地空间直角坐标分别为[X_h1 Y_h1 Z_h1]^T、[X_hq Y_hq Z_hq]^T，将大地空间直角坐标分别转换成大地坐标[L_h1 B_h1 H_h1]^T、[L_hq B_hq H_hq]^T；根据式(2)判断T_s时刻GNSS卫星到GNSS接收机的信号传播路径是否属于电离层层析待反演区域：若满足式(2)，则T_s时刻GNSS卫星到GNSS接收机的信号传播路径属于电离层层析待反演区域；若不满足式(2)，则T_s时刻GNSS卫星到GNSS接收机的信号传播路径不属于电离层层析待反演区域；

$\left\{\begin{array}{c}{L}_1<\&lsqb;L_{h1},L_{hq}\&rsqb;<L_n\\ B_1<\&lsqb;B_{h1},B_{hq}\&rsqb;<B_m\end{array}\right.---\left(2\right)$

其中s＝1,2,…,t；

S4.3：确定折半搜索算法在三维像素经度面和纬度面的目标位置：

根据S4.2，若T_s时刻GNSS卫星到GNSS接收机的信号传播路径属于电离层层析待反演区域，将S2.1中满足式(3)的经度L_i和纬度B_j分别作为折半搜索算法在三维像素经度面和纬度面的目标位置；

$\left\{\begin{array}{c}{L}_{h1}<L_i<L_{hq}\\ B_{h1}<B_j<B_{hq}\end{array}\right.---\left(3\right).$

进一步，所述步骤S3包括以下的步骤：

S3.1：在电离层层析待反演时段[T₁ T_t]内，准备待反演区域内以及周边的GNSS接收机的大地空间直角坐标；

S3.2：获取电离层层析待反演时段[T₁ T_t]内GNSS卫星星历数据，并根据星历计算待反演时段内GNSS卫星不同时刻的大地空间直角坐标。

进一步，所述步骤S4中计算GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点包括以下的步骤：

S4.4：初始化折半搜索算法的参数：将T_s时刻GNSS卫星的大地直角坐标作为上边界将GNSS接收机的大地直角坐标作为下边界设置目标位置LBH，依据实际电离层层析精度需求设置阀值e，其中s＝1,2,…,t；

S4.5：计算探测位置，并判定是否继续折半搜索：根据式(4)计算探测位置的大地空间直角坐标将探测位置的大地空间直角坐标转换成大地坐标若目标位置LBH与探测位置相应大地坐标值P.LBH的差的绝对值|P.LBH-LBH|＜e，进行步骤S4.7；若|P.LBH-LBH|≥e，进行步骤S4.6；

${\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_P^T=\left({\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_U^T+{\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_L^T\right)/2---\left(4\right)$

S4.6：更新上下边界条件：若P.LBH-LBH＜0，则上边界保持不变，将探测位置的大地空间直角坐标作为下边界；若P.LBH-LBH＞0，则下边界保持不变，将探测位置的大地空间直角坐标作为上边界，然后回到步骤S4.5；

S4.7：停止折半搜索，则当前的探测位置被认为是T_s时刻GNSS信号传播路径与目标位置的交点。

进一步，所述步骤S5包括以下的步骤：

S5.1：对T_s时刻GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点进行排序：依据交点的坐标大小进行降序或者升序排列；

S5.2：根据S5.1的排序结果，确定相邻交点所穿过的三维像素编号；

将相邻交点的大地空间直角坐标[X_fl Y_fl Z_fl]^T、[X_fu Y_fu Z_fu]^T分别转换成大地坐标[L_fl B_fl H_fl]^T、[L_fu B_fu H_fu]^T，计算相邻交点的中点的大地坐标[L B H]^T，若L、B和H满足式(5)，则所穿过的三维像素左下角的经度为L_i、纬度为B_j和椭球高为H_k，根据式(1)确定这个三维像素编号f；

$\left\{\begin{array}{c}{L}_i<L<L_{i+1}\\ B_j<B<B_{j+1}\\ H_k<H<H_{k+1}\end{array}\right.---\left(5\right)$

其中i＝1,2,…,n-1、j＝1,2,…,m-1、k＝1,2,…,q-1。

S5.3：根据S5.1的排序结果，计算GNSS信号传播路径在三维像素f的截距：

相邻交点的大地空间坐标分别为[X_fl Y_fl Z_fl]^T、[X_fu Y_fu Z_fu]^T，根据式(6)计算GNSS信号传播路径在三维像素f的截距d_f；

$d_f=\sqrt{{(X_{fu}-X_{fl})}^2+{(Y_{fu}-Y_{fl})}^2+{(Z_{fu}-Z_{fl})}^2}---\left(6\right).$

进一步，所述步骤S6包括以下的步骤：

S6.1：构造T_s时刻的投影向量A_s：

若GNSS信号传播路径穿过三维像素f，则GNSS信号传播路径在三维像素f的截距为d_f，即A_s(f)＝d_f；若GNSS信号传播路径没有穿过三维像素f，则A_s(f)＝0，f＝1,2,…,(n-1)·(m-1)·(q-1)，s＝1,2,…,t；

S6.2：构造在电离层层析待反演时段[T₁ T_t]内待反演区域的投影矩阵A：

构造电离层层析待反演的时刻[T₁ T₂ ... T_t]相应的投影向量[A₁ A₂ ... A_t],则式(7)就是电离层层析待反演时段[T₁ T_t]内待反演区域的投影矩阵A；

$A=\left[\begin{array}{c}{A}_1\\ A_2\\ .\\ .\\ .\\ A_n\end{array}\right]---\left(7\right)$

有益效果：本发明通过基于椭球的GNSS电离层层析获取投影矩阵，克服了目前基于圆球提取投影矩阵方法时将地球椭球近似圆球的问题，保证了相对于地球椭球的投影矩阵能真实反映出GNSS信号传播路径在反演像素内的真实截距。并且，本发明提出的基于椭球的GNSS电离层层析获取投影矩阵有利于分析高分辨率电离层时空变化。

附图说明

图1为本发明的流程图；

图2为本发明的电离层层析中GNSS信号传播路径穿过三维像素的二维示意图。

具体实施方式

下面结合具体实施方式对本发明的技术方案作进一步的介绍。

本发明公开了一种基于椭球的GNSS电离层层析投影矩阵获取方法，如图1所示，包括以下的步骤：

S1：确定电离层层析待反演区域及待反演时段；

步骤S1包括以下的步骤：

S1.1：确定电离层层析待反演区域；

确定电离层层析待反演区域的经度范围[116.0° 123.0°]；纬度范围[30.0°36.0°]；椭球高范围[70km 1000km]；

S1.2：确定电离层层析待反演时段[05:01:15UT 05:01:45UT]，电离层层析待反演的时刻为[05:01:15UT 05:01:30UT 05:01:45UT]。

S2：将待反演区域离散化为一系列三维像素；

步骤S2包括以下步骤：

S2.1：根据电离层层析待反演区域的经度范围、纬度范围和椭球高范围，将待反演区域离散化为一系列的三维像素，即：

经度被离散化为(116.0+(i-1))°E，(i＝1,2,…,8)；

纬度被离散化为(30.0+(j-1))°N，(j＝1,2,…,7)；

椭球高被离散化为(70+30·(k-1))km，(k＝1,2,…,32)；

S2.2：对三维像素按照经度从西到东、纬度从南到北和椭球高从下到上的方式进行编号：当三维像素的左下角的经度为117.0°E、纬度为31.0°N和椭球高为100km，可知i＝2，j＝2，k＝2，根据式(1)确定这个三维像素编号f为：f＝(2-1)×((7-1)×(8-1))+(2-1)×(8-1)+2＝51。

S3：准备电离层层析待反演区域及待反演时段内的GNSS卫星和GNSS接收机的大地空间直角坐标；

步骤S3包括以下步骤：

S3.1：在电离层层析待反演时段[05:01:15UT 05:01:45UT]内，准备待反演区域GNSS接收机的坐标如表1；

表1 GNSS接收机坐标

S3.2：获取电离层层析待反演时段[05:01:15UT 05:01:45UT]内GNSS卫星不同时刻的大地空间直角坐标如表2。

表2 卫星不同时刻的大地空间直角坐标

S4：计算待反演时段内GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点：根据图2所示的GNSS信号传播路径穿过三维像素的二维示意图，首先确定折半搜索算法在三维像素各个面的目标位置，然后计算GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点。

步骤S4中确定折半搜索算法在三维像素各个面的目标位置包括以下的步骤：

S4.1：确定折半搜索算法在三维像素椭球高面的目标位置：

折半搜索算法在三维像素椭球高面的目标位置分别为S2.1中被间隔的椭球高(70+30·(k-1))km，(k＝1,2,…,32)；

S4.2：判断待反演时段内GNSS卫星到GNSS接收机的信号传播路径是否属于电离层层析待反演区域；

根据折半搜索算法计算出T_s时刻GNSS信号传播路径与椭球高面70km、1000km的交点的大地空间直角坐标如表3所示；

表3 GNSS信号传播路径与椭球高面70km、1000km的交点大地空间直角坐标和其大地坐标

根据式(2)判断T_s时刻GNSS卫星到GNSS接收机的信号传播路径是否属于电离层层析待反演区域；

S4.3：确定折半搜索算法在三维像素经度面和纬度面的目标位置：

将在S2.1中满足式(3)的被间隔经度和纬度分别作为折半搜索算法在三维像素经度面和纬度面的目标位置。

5:01:15UT时刻：经度面目标位置为空；纬度面目标位置为[33.0°N 34.0°N 35.0°N]；

5:01:30UT时刻：经度面目标位置为空；纬度面目标位置为[33.0°N 34.0°N 35.0°N]；

5:01:45UT时刻：经度面目标位置为空；纬度面目标位置为[33.0°N 34.0°N 35.0°N]；

S4.4：初始化折半搜索算法的参数：将T_s时刻GNSS卫星的大地直角坐标作为上边界，将GNSS接收机的大地直角坐标作为下边界，设置目标位置LBH，依据实际电离层层析精度需求设置阀值e；

其中，经度面和纬度面的阀值e设置为(5.7296×10^-8)°；椭球高度面的阀值e设置为0.001m；

S4.5：计算探测位置，并判定是否继续折半搜索：根据式(4)计算探测位置的大地空间直角坐标将探测位置的大地空间直角坐标转换成大地坐标若目标位置LBH与探测位置相应大地坐标值P.LBH的差的绝对值|P.LBH-LBH|＜e，进行步骤S4.7；若|P.LBH-LBH|≥e，进行步骤S4.6；

${\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_P^T=\left({\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_U^T+{\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_L^T\right)/2---\left(4\right)$

S4.6：更新上下边界条件：若P.LBH-LBH＜0，则上边界保持不变，将探测位置的大地空间直角坐标作为下边界；若P.LBH-LBH＞0，则下边界保持不变，将探测位置的大地空间直角坐标作为上边界，然后回到步骤S4.5；

S4.7：停止折半搜索，则当前的探测位置被认为是T_s时刻GNSS信号传播路径与目标位置的交点。

S5：计算三维像素内GNSS卫星信号的传播路径的截距长度；

步骤S5包括以下步骤：

S5.1：对GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面交点排序：5:01:15UT时刻依据交点的X轴坐标大小进行升序排列结果如表4。

S5.2：根据S5.1的排序结果，确定相邻交点所穿过的三维像素编号。5：01：15UT时刻的相邻交点所穿过的三维像素编号如表4；

S5.3：根据S5.1的排序结果，计算GNSS信号传播路径在三维像素f的截距。5：01：15UT时刻的GNSS信号传播路径在三维像素f的截距如表4；

表4 (5:01:15)UT时刻交点坐标、相邻交点的像素编号及像素的截距

S6：构成投影矩阵；

步骤S6包括以下步骤：

S6.1：构造T_s时刻的投影向量A_s：

根据S2.1，T_s时刻的投影向量A_s为1×1302的行向量。若GNSS信号传播路径穿过在三维像素f，根据S5.3可知在三维像素f的截距为d_f，则A_s(f)＝d_f；若GNSS信号传播路径没有穿过在三维像素f，则A_s(f)＝0，(f＝1,2,…,1302)；

S6.2：构造在电离层层析待反演时段[05:01:15UT 05:01:45UT]内待反演区域的投影矩阵A：

根据S6.1构造电离层层析待反演的时刻[05:01:15UT 05:01:30UT 05:01:45UT]相应的投影向量[A₁ A₂ A₃],则电离层层析待反演时段[05:01:15UT 05:01:45UT]内待反演区域的投影矩阵为：

$A=\left[\begin{array}{c}{A}_1\\ A_2\\ A_3\end{array}\right]$

因此，A就是在待反演区域的经度范围[116.0° 123.0°]、纬度范围[30.0°36.0°]、椭球高范围[70km 1000km]和电离层层析待反演时段[05:01:15UT 05:01:45UT]内的投影矩阵。通过基于椭球的GNSS电离层层析获取投影矩阵A，保证了GNSS信号传播路径在相对地球椭球的反演像素内真实截距，并能够快速定位截距所在像素的位置。本发明提出的基于椭球的GNSS电离层层析获取投影矩阵有利于分析高分辨率电离层时空变化。

Claims (8)

1.一种基于椭球的GNSS电离层层析投影矩阵获取方法，其特征在于：包括以下的步骤：

S1：确定电离层层析待反演区域及待反演时段；

S2：将待反演区域离散化为一系列三维像素；

S3：准备电离层层析待反演区域及待反演时段内的GNSS卫星和GNSS接收机的大地空间直角坐标；

S4：计算待反演时段内GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点：首先确定折半搜索算法在三维像素各个面的目标位置，然后计算GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点；

S5：计算三维像素内GNSS卫星信号的传播路径的截距长度；

S6：构成投影矩阵。

2.根据权利要求1所述的基于椭球的GNSS电离层层析投影矩阵获取方法，其特征在于：所述步骤S1包括以下的步骤：

S1.1：确定电离层层析待反演区域；

确定电离层层析待反演区域的经度范围[L₁ L_n]、纬度范围[B₁ B_m]和椭球高范围[H₁H_q]；

S1.2：确定电离层层析待反演时段[T₁ T_t]，电离层层析待反演的时刻为[T₁ T₂ ...T_t]。

3.根据权利要求1所述的基于椭球的GNSS电离层层析投影矩阵获取方法，其特征在于：所述步骤S2包括以下的步骤：

S2.1：根据电离层层析待反演区域的经度范围、纬度范围和椭球高范围，将待反演区域离散化为一系列的三维像素，即：

经度被离散化为[L₁ L₂ ... L_n]；

纬度被离散化为[B₁ B₂ ... B_m]；

椭球高被离散化为[H₁ H₂ ... H_q]；

S2.2：对三维像素按照经度从西到东、纬度从南到北和椭球高从下到上的方式进行编号：当三维像素的左下角的经度为L_i、纬度为B_j和椭球高为H_k，根据式(1)确定这个三维像素编号f：

f＝(k-1)×((m-1)×(n-1))+(j-1)×(n-1)+i (1)

其中i＝1,2,…,n、j＝1,2,…,m、k＝1,2,…,q。

4.根据权利要求3所述的基于椭球的GNSS电离层层析投影矩阵获取方法，其特征在于：所述步骤S4中确定折半搜索算法在三维像素各个面的目标位置包括以下的步骤：

S4.1：确定折半搜索算法在三维像素椭球高面的目标位置：

折半搜索算法在三维像素椭球高面的目标位置分别为S2.1中被间隔的椭球高[H₁ H₂... H_q]；

S4.2：判断待反演时段内GNSS卫星到GNSS接收机的信号传播路径是否属于电离层层析待反演区域；

根据折半搜索算法计算出T_s时刻GNSS信号传播路径与椭球高面H₁、H_q的交点的大地空间直角坐标分别为[X_h1 Y_h1 Z_h1]^T、[X_hq Y_hq Z_hq]^T，将大地空间直角坐标分别转换成大地坐标[L_h1 B_h1 H_h1]^T、[L_hq B_hq H_hq]^T；根据式(2)判断T_s时刻GNSS卫星到GNSS接收机的信号传播路径是否属于电离层层析待反演区域：若满足式(2)，则T_s时刻GNSS卫星到GNSS接收机的信号传播路径属于电离层层析待反演区域；若不满足式(2)，则T_s时刻GNSS卫星到GNSS接收机的信号传播路径不属于电离层层析待反演区域；

$\left\{\begin{array}{c}{L}_1<\&lsqb;L_{h1},L_{hq}\&rsqb;<L_n\\ B_1<\&lsqb;B_{h1},B_{hq}\&rsqb;<B_m\end{array}\right.---\left(2\right)$

其中s＝1,2,…,t；

S4.3：确定折半搜索算法在三维像素经度面和纬度面的目标位置：

根据S4.2，若T_s时刻GNSS卫星到GNSS接收机的信号传播路径属于电离层层析待反演区域，将S2.1中满足式(3)的经度L_i和纬度B_j分别作为折半搜索算法在三维像素经度面和纬度面的目标位置；

$\left\{\begin{array}{c}{L}_{h1}<L_i<L_{hq}\\ B_{h1}<B_j<B_{hq}\end{array}\right.---\left(3\right).$

5.根据权利要求1所述的基于椭球的GNSS电离层层析投影矩阵获取方法，其特征在于：所述步骤S3包括以下的步骤：

S3.1：在电离层层析待反演时段[T₁ T_t]内，准备待反演区域内以及周边的GNSS接收机的大地空间直角坐标；

S3.2：获取电离层层析待反演时段[T₁ T_t]内GNSS卫星星历数据，并根据星历计算待反演时段内GNSS卫星不同时刻的大地空间直角坐标。

6.根据权利要求1所述的基于椭球的GNSS电离层层析投影矩阵获取方法，其特征在于：所述步骤S4中计算GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点包括以下的步骤：

S4.4：初始化折半搜索算法的参数：将T_s时刻GNSS卫星的大地直角坐标作为上边界将GNSS接收机的大地直角坐标作为下边界设置目标位置LBH，依据实际电离层层析精度需求设置阀值e，其中s＝1,2,…,t；

S4.5：计算探测位置，并判定是否继续折半搜索：根据式(4)计算探测位置的大地空间直角坐标将探测位置的大地空间直角坐标转换成大地坐标若目标位置LBH与探测位置相应大地坐标值P.LBH的差的绝对值|P.LBH-LBH|＜e，进行步骤S4.7；若|P.LBH-LBH|≥e，进行步骤S4.6；

${\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_P^T=\left({\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_U^T+{\left[\begin{array}{ccc}X& Y& Z\end{array}\right]}_L^T\right)/2---\left(4\right)$

S4.6：更新上下边界条件：若P.LBH-LBH＜0，则上边界保持不变，将探测位置的大地空间直角坐标作为下边界；若P.LBH-LBH＞0，则下边界保持不变，将探测位置的大地空间直角坐标作为上边界，然后回到步骤S4.5；

S4.7：停止折半搜索，则当前的探测位置被认为是T_s时刻GNSS信号传播路径与目标位置的交点。

7.根据权利要求3所述的基于椭球的GNSS电离层层析投影矩阵获取方法，其特征在于：所述步骤S5包括以下的步骤：

S5.1：对T_s时刻GNSS卫星到GNSS接收机的信号传播路径与三维像素各个面的交点进行排序：依据交点的坐标大小进行降序或者升序排列；

S5.2：根据S5.1的排序结果，确定相邻交点所穿过的三维像素编号；

将相邻交点的大地空间直角坐标[X_fl Y_fl Z_fl]^T、[X_fu Y_fu Z_fu]^T分别转换成大地坐标[L_flB_fl H_fl]^T、[L_fu B_fu H_fu]^T，计算相邻交点的中点的大地坐标[L B H]^T，若L、B和H满足式(5)，则所穿过的三维像素左下角的经度为L_i、纬度为B_j和椭球高为H_k，根据式(1)确定这个三维像素编号f；

$\left\{\begin{array}{c}{L}_i<L<L_{i+1}\\ B_j<B<B_{j+1}\\ H_k<H<H_{k+1}\end{array}\right.---\left(5\right)$

其中i＝1,2,…,n-1、j＝1,2,…,m-1、k＝1,2,…,q-1；

S5.3：根据S5.1的排序结果，计算GNSS信号传播路径在三维像素f的截距：

相邻交点的大地空间坐标分别为[X_fl Y_fl Z_fl]^T、[X_fu Y_fu Z_fu]^T，根据式(6)计算GNSS信号传播路径在三维像素f的截距d_f；

$d_f=\sqrt{{(X_{fu}-X_{fl})}^2+{(Y_{fu}-Y_{fl})}^2+{(Z_{fu}-Z_{fl})}^2}---\left(6\right).$

8.根据权利要求7所述的基于椭球的GNSS电离层层析投影矩阵获取方法，其特征在于：所述步骤S6包括以下的步骤：

S6.1：构造T_s时刻的投影向量A_s：

若GNSS信号传播路径穿过三维像素f，则GNSS信号传播路径在三维像素f的截距为d_f，即A_s(f)＝d_f；若GNSS信号传播路径没有穿过三维像素f，则A_s(f)＝0，f＝1,2,…,(n-1)·(m-1)·(q-1)，s＝1,2,…,t；

S6.2：构造在电离层层析待反演时段[T₁ T_t]内待反演区域的投影矩阵A：

构造电离层层析待反演的时刻[T₁ T₂ ... T_t]相应的投影向量[A₁ A₂ ... A_t],则式(7)就是电离层层析待反演时段[T₁ T_t]内待反演区域的投影矩阵A；

$A=\left[\begin{array}{c}{A}_1\\ A_2\\ \begin{array}{c}.\\ .\\ .\end{array}\\ A_n\end{array}\right]---\left(7\right).$