SBAS 电离层修正中的双线性插值方法

文章目录

线性插值
四点双线性插值
三点插值
还有一步

点击这里可以阅读rtklib中的SBAS相关实现解读，其中使用SBAS进行电离层修正时，首先要计算电离层切入点。但是有一个问题，SBAS不可能下发电离层上无限个点的修正信息，而是一个一个的特征点，这样整个电离层球面就会被分成一个个的网格区域。正因如此，电离层切入点的修正值需要使用插值获取。

线性插值

在正式开始我们的插值工作之前，有必要花费半分钟的时间熟悉一下线性插值。如下图，我们如果知道两个点的坐标(x1,y1)(x_1,y_1)(x1,y1)和(x0,y0)(x_0,y_0)(x0,y0)，那么我们就可以通过线性插值计算出横坐标xxx对应的yyy值。

y=x1−xx1−x0y0+x−x0x1−x0y1=αy0+βy1y=\frac{x_1-x}{x_1-x_0}y_0+\frac{x-x_0}{x_1-x_0}y_1=\alpha y_0+\beta y_1 y=x1−x0x1−xy0+x1−x0x−x0y1=αy0+βy1
其中 α+β=1\alpha+\beta=1α+β=1

四点双线性插值

在低纬度地区，我们算出的电离层切入点会位于一个近似于正方形的网格中，如下图。

此时我们已知四个顶点的电离层延迟，要计算切入点的电离层延迟就需要使用双线性插值的方法。

首先沿x方向插值出f(x,y1)f(x,y_1)f(x,y1)和f(x,y2)f(x,y_2)f(x,y2)连个中间点，插值方法采用上边的线性插值，结果如下
f(x,y1)=x2−xx2−x1Q11+x−x1x2−x1Q21=αxQ11+βxQ21f(x,y_1)=\frac{x_2-x}{x_2-x_1}Q_{11}+\frac{x-x_1}{x_2-x_1}Q_{21}=\alpha_{x} Q_{11}+\beta_{x} Q_{21} f(x,y1)=x2−x1x2−xQ11+x2−x1x−x1Q21=αxQ11+βxQ21

f(x,y2)=x2−xx2−x1Q12+x−x1x2−x1Q22=αxQ12+βxQ22f(x,y_2)=\frac{x_2-x}{x_2-x_1}Q_{12}+\frac{x-x_1}{x_2-x_1}Q_{22}=\alpha_{x} Q_{12}+\beta_{x} Q_{22} f(x,y2)=x2−x1x2−xQ12+x2−x1x−x1Q22=αxQ12+βxQ22
现在通过(x,f(x,y1))(x,f(x,y_1))(x,f(x,y1))和x,f(x(x,y2))x,f(x(x,y_2))x,f(x(x,y2))同样适用线性插值的方法，插值出f(x,y)f(x,y)f(x,y)
f(x,y)=y2−yy2−y1f(x,y1)+y−y1y2−y1f(x,y2)=αyf(x,y1)+βyf(x,y2)f(x,y)=\frac{y_2-y}{y_2-y_1}f(x,y_1)+\frac{y-y_1}{y_2-y_1}f(x,y_2)=\alpha_{y} f(x,y_1)+\beta_{y} f(x,y_2) f(x,y)=y2−y1y2−yf(x,y1)+y2−y1y−y1f(x,y2)=αyf(x,y1)+βyf(x,y2)
f(x,y)=[αxβx][Q11Q12Q21Q22][αyβy]f(x,y)=\begin{bmatrix} \alpha_x & \beta_x \end{bmatrix} \begin{bmatrix} Q_{11} & Q_{12} \\ Q_{21} & Q_{22} \end{bmatrix} \begin{bmatrix} \alpha_y\\ \beta_y \end{bmatrix} f(x,y)=[αxβx][Q11Q21Q12Q22][αyβy]
如果网格是正方形的，经过单位化，那么双线性插值的算法简化为下边的样子，
f(x,y)=[1−xx][Q11Q12Q21Q22][1−yy]f(x,y)=\begin{bmatrix} 1-x & x \end{bmatrix} \begin{bmatrix} Q_{11} & Q_{12} \\ Q_{21} & Q_{22} \end{bmatrix} \begin{bmatrix} 1-y\\ y \end{bmatrix} f(x,y)=[1−xx][Q11Q21Q12Q22][1−yy]
下边这两行没有删掉的代码就是单位化的一种体现，

static void searchigp(gtime_t time, const double *pos, const sbsion_t *ion,const sbsigp_t **igp, double *x, double *y)
{*x=(lon-lonp[0])/5.0;*y=(lat-latp[0])/5.0;
}

最终求出的四个点的权值和我们的推导是一致的，

extern int sbsioncorr(gtime_t time, const nav_t *nav, const double *pos,const double *azel, double *delay, double *var)
{/* weight of igps */if (igp[0]&&igp[1]&&igp[2]&&igp[3]) {w[0]=(1.0-x)*(1.0-y); w[1]=(1.0-x)*y; w[2]=x*(1.0-y); w[3]=x*y;}

三点插值

在高纬度地区我们面临的是以下局面。

这时我们就要用三个顶点进行插值，对于三点插值我们采用重心法。

同样，经过代码处理，三角形可退化为等腰直角三角形。三角形的重心坐标与三个子三角形的面积成正比。由于等腰直角三角形的特殊性，那么这个比例又可以退化为切入点到两条直角边的距离。

最后一个三角形的面积（最后一个点的权值）我们是不用计算的，因为三个点的权值之和为1.只需要用1减去已知的两个权值就可以了。

这个等腰直角三角形分为以下四种情况。从第一张图，我们可以看出，w3=xw_3=xw3=x,w0=1−yw_0=1-yw0=1−y,w1=1−w0−w3w_1=1-w_0-w_3w1=1−w0−w3

同样，如果是下图这种情况，w1=1−xw_1=1-xw1=1−x,w2=1−yw_2=1-yw2=1−y,w3=1−w1−w2w_3=1-w_1-w_2w3=1−w1−w2

如果是下图这种情况，w1=yw_1=yw1=y,w2=xw_2=xw2=x,w0=1−w1−w2w_0=1-w_1-w_2w0=1−w1−w2

最后是下图这种情况，w0=1−xw_0=1-xw0=1−x,w3=yw_3=yw3=y,w2=1−w1−w3w_2=1-w_1-w_3w2=1−w1−w3

以下是rtklib中三点插值的处理，和我们推导的结果是一致的。

    else if (igp[0]&&igp[1]&&igp[2]) {w[1]=y; w[2]=x;if ((w[0]=1.0-w[1]-w[2])<0.0) err=1;}else if (igp[0]&&igp[2]&&igp[3]) {w[0]=1.0-x; w[3]=y;if ((w[2]=1.0-w[0]-w[3])<0.0) err=1;}else if (igp[0]&&igp[1]&&igp[3]) {w[0]=1.0-y; w[3]=x;if ((w[1]=1.0-w[0]-w[3])<0.0) err=1;}else if (igp[1]&&igp[2]&&igp[3]) {w[1]=1.0-x; w[2]=1.0-y;if ((w[3]=1.0-w[1]-w[2])<0.0) err=1;}

还有一步

另外，我们发现无论是电离层修正结果还是修正方差，都是要考虑因子fpfpfp的，这又是什么鬼？

*delay*=fp; *var*=fp*fp;

首先，下边是这个参数的计算方法，
Fp=[1−(RecosERe+hI)2]−12F_p=[1-(\frac{R_ecosE}{R_e+h_I})^2]^{-\frac{1}{2}} Fp=[1−(Re+hIRecosE)2]−21
其中EEE是卫星仰角，ReR_eRe 为地球半径，hIh_IhI电离层高度，取350Km.

看出来了吧，延迟修正不仅与电离层切入点的位置有关，还与卫星仰角有关。仰角越小FpF_pFp越大，这是应为仰角越小穿过电力层的距离也就越长。