当前位置:   article > 正文

ddres( ) 组站星双差方程和设计矩阵_rtklib ddres

rtklib ddres

1 ddres( )参数介绍

rtklib中进行的单频解算

双差观测值,单差的模糊度         单频点双差
    DD (double-differenced) phase/code residuals ------------------------------
    x 模糊度
    P 方差-协方差阵
    sat 共识卫星列表
    ns 共识卫星数量
    y   残差阵,数量:ns*nf*2(每颗卫星的每个频率,×2是伪距和载波)
    e   dx dy dz设计矩阵,数量:ns*3
    freq    频率,数量:ns*nf  每个频点的伪距和载波频率一样,而且伪距也用不到频率,是载波求波长时用的
    v   完整残差阵(包括了模糊度)
    H   dx dy dz N1 N2 N2等模糊度 设计矩阵
    R   方差-协方差阵

  1. static int ddres(rtk_t *rtk, const nav_t *nav, double dt, const double *x,
  2. const double *P, const int *sat, double *y, double *e,
  3. double *azel, double *freq, const int *iu, const int *ir,
  4. int ns, double *v, double *H, double *R, int *vflg)

2 代码解读

计算基线的长度,计算基准站和流动站方差时会用到

  1. bl=baseline(x,rtk->rb,dr);
  2. ecef2pos(x,posu); ecef2pos(rtk->rb,posr);

Ri是参考卫星方差,RJ是非参考卫星方差。用于后面计算双差后的观测值方差--协方差阵

tropu、tropr等是对流层、电离层参数,在此我们不作讨论

  1. Ri=mat(ns*nf*2+2,1); Rj=mat(ns*nf*2+2,1); im=mat(ns,1);
  2. tropu=mat(ns,1); tropr=mat(ns,1); dtdxu=mat(ns,3); dtdxr=mat(ns,3);

初始化伪距和载波相位残差为0,每个卫星每个频点的,用于以后记录双差残差        

rtk->ssat[i];访问到了每颗卫星

double resp[NFREQ]; residuals of pseudorange (m) 
double resc[NFREQ]; residuals of carrier-phase (m) 

MAXSAT;最大卫星数,不管观没观测到,都初始化。这也造成了空间上的浪费

NFREQ;最大频率

  1. for (i=0;i<MAXSAT;i++) for (j=0;j<NFREQ;j++) {
  2. rtk->ssat[i].resp[j]=rtk->ssat[i].resc[j]=0.0;
  3. }

对流层和电离层处理,暂不考虑

  1. for (i=0;i<ns;i++) {
  2. if (opt->ionoopt>=IONOOPT_EST) {
  3. im[i]=(ionmapf(posu,azel+iu[i]*2)+ionmapf(posr,azel+ir[i]*2))/2.0;
  4. }
  5. if (opt->tropopt>=TROPOPT_EST) {
  6. tropu[i]=prectrop(rtk->sol.time,posu,0,azel+iu[i]*2,opt,x,dtdxu+i*3);
  7. tropr[i]=prectrop(rtk->sol.time,posr,1,azel+ir[i]*2,opt,x,dtdxr+i*3);
  8. }
  9. }

之后就进入一个大循环,最外层是先遍历每个卫星系统

for (m=0;m<6;m++)    /* m=0:GPS/SBS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN */

二重循环,根据所选择的定位模式,进行相应的频点遍历。nf*2是因为每个系统的每个频点(*2 包括伪距和载波)

只利用伪距观测值,nf---2*nf,因为伪距残差是在y[]后半

是载波(动态定位模式及以后),那我们从0--2*nf,因为同时使用伪距和载波残差

for (f=opt->mode>PMODE_DGPS?0:nf;f<nf*2;f++)

ns---共识卫星数量

这个就是找参考卫星的,在这里是把高度角最高的那个作为参考卫星。我们也可以改写代码,设置其他的规则。参考卫星下标 i

  1. for (i=-1,j=0;j<ns;j++) {
  2. sysi=rtk->ssat[sat[j]-1].sys;
  3. if (!test_sys(sysi,m)) continue;//当前共识卫星的系统和当前遍历的系统是否一样
  4. if (!validobs(iu[j],ir[j],f,nf,y)) continue; //评价是否可以
  5. if (i<0||azel[1+iu[j]*2]>=azel[1+iu[i]*2]) i=j;//选择排序的那种,i最终是参考卫星的index
  6. }

又是一个大循环,遍历所有的共识卫星,这次是要计算双差残差和设计矩阵。组双差的话,两台接收机是要都观测到相同的卫星的,所以需要共识卫星处理。

以下代码均是在此循环中执行

  1. for (j=0;j<ns;j++) {
  2. }

两个不同的卫星,所以i≠j

freq[f%nf+iu[i]*nf];取频率

f区间为[0,2*nf],所以取余%

freqi,freqj会在电离层、对流层用到

  1. if (i==j) continue;
  2. sysi=rtk->ssat[sat[i]-1].sys;//参考卫星系统
  3. sysj=rtk->ssat[sat[j]-1].sys;//非参考卫星系统
  4. freqi=freq[f%nf+iu[i]*nf];
  5. freqj=freq[f%nf+iu[j]*nf];
  6. if (!test_sys(sysj,m)) continue;//判断和当前遍历系统一样否
  7. if (!validobs(iu[j],ir[j],f,nf,y)) continue;//检验数据有效yy[]

设计矩阵初始化,nx为待估量的个数

  1. if (H) {
  2. Hi=H+nv*rtk->nx;
  3. for (k=0;k<rtk->nx;k++) Hi[k]=0.0;
  4. }

求双差残差,

(y[f+iu[i]*nf*2]-y[f+ir[i]*nf*2]) 共识卫星站间单差,是参考卫星i

y[f+iu[j]*nf*2]-y[f+ir[j]*nf*2]   共识卫星站间单差,是参考卫星j

两个一减,就是双差残差了

 v[nv]=(y[f+iu[i]*nf*2]-y[f+ir[i]*nf*2])-(y[f+iu[j]*nf*2]-y[f+ir[j]*nf*2]);

设计矩阵,e[ ]是dx dy dz的系数,已在站间单差函数中求得,我们只考虑单频,e[ ]中只有单频的信息

e[]是dx dy dz设计矩阵,大小:流动站观测到的卫星+基准站观测到的卫星。 
因为在relpos()函数中,站间单差[先是算基准站,后是流动站],设计矩阵都填到e里面了。所以大小如此

设计矩阵=非参考卫星-参考卫星

每次循环都会产生一组双差dx dy dz的系数。但最最外面还有一个频率循环,那么不同频的dx dy dz的系数是一样的吗?这我还不知道 

这两个会是一样吗?留到以后回答,或者有谁知道

我的猜想是可能一样,因为由计算系数的公式可知,这个与卫星位置和接收机概略位置有关

 

  1. if (H) {
  2. for (k=0;k<3;k++) {
  3. Hi[k]=-e[k+iu[i]*3]+e[k+iu[j]*3];
  4. }
  5. }

IB()不懂是什么。 哪个卫星哪个频点的模糊度位置索引??

看上去像是索引,但为什么x[ ] Hi[ ]都可以用。一个是设计矩阵位置,一个是对应模糊度位置啥的

  1. if (opt->ionoopt!=IONOOPT_IFLC) {
  2. v[nv]-=CLIGHT/freqi*x[IB(sat[i],f,opt)]-
  3. CLIGHT/freqj*x[IB(sat[j],f,opt)];
  4. if (H) {
  5. Hi[IB(sat[i],f,opt)]= CLIGHT/freqi;
  6. Hi[IB(sat[j],f,opt)]=-CLIGHT/freqj;
  7. }
  8. }
  9. else {
  10. //根据站星双差方差,还要减去(模糊度的derT),见理论。而在y[]的计算中,还没有考虑模糊度,在这里给补上
  11. //x[] 模糊度
  12. v[nv]-=x[IB(sat[i],f,opt)]-x[IB(sat[j],f,opt)];
  13. if (H) {
  14. Hi[IB(sat[i],f,opt)]= 1.0; //IB() 哪个卫星哪个频点的模糊度位置索引???不懂
  15. Hi[IB(sat[j],f,opt)]=-1.0;
  16. }
  17. }
  18. }

Hi[IB(sat[i],f,opt)]= 1.0 这里是在补设计矩阵,参考卫星系数为-1。而这里为1,后续应该在某个计算中会改掉

v[ ]是计算好的双差残差,将其存到ssat结构体中

之后检验残差是否在容许范围内 maxinno是设置的标准

  1. if (f<nf) rtk->ssat[sat[j]-1].resc[f ]=v[nv];
  2. else rtk->ssat[sat[j]-1].resp[f-nf]=v[nv];
  3. /* test innovation 检查残差是否符合标准 maxinno */
  4. if (opt->maxinno>0.0&&fabs(v[nv])>opt->maxinno) {
  5. if (f<nf) {
  6. rtk->ssat[sat[i]-1].rejc[f]++;
  7. rtk->ssat[sat[j]-1].rejc[f]++;
  8. }
  9. errmsg(rtk,"outlier rejected (sat=%3d-%3d %s%d v=%.3f)\n",
  10. sat[i],sat[j],f<nf?"L":"P",f%nf+1,v[nv]);
  11. continue;
  12. }

每颗卫星每个频率的观测值噪声(方差)

  1. Ri[nv]=varerr(sat[i],sysi,azel[1+iu[i]*2],bl,dt,f,opt);
  2. Rj[nv]=varerr(sat[j],sysj,azel[1+iu[j]*2],bl,dt,f,opt);

设置标志,这颗卫星这个频点的伪距/载波观测值有效。标记一下

  1. if (opt->mode>PMODE_DGPS) {
  2. if (f<nf) rtk->ssat[sat[i]-1].vsat[f]=rtk->ssat[sat[j]-1].vsat[f]=1;
  3. }
  4. else {
  5. rtk->ssat[sat[i]-1].vsat[f-nf]=rtk->ssat[sat[j]-1].vsat[f-nf]=1;
  6. }

独特的编码,由此码可知,双差的具体信息。哪两个卫星作差,载波还是伪距,第几个频率

vflg[nv++]=(sat[i]<<16)|(sat[j]<<8)|((f<nf?0:1)<<4)|(f%nf);

至此,共识卫星遍历结束


每个小频率有几颗卫星(残差有多少个)

nb[b]++;


至此,频点遍历结束

关于基线约束的

  1. if (opt->mode==PMODE_MOVEB&&constbl(rtk,x,P,v,H,Ri,Rj,nv)) {
  2. vflg[nv++]=3<<4;
  3. nb[b++]++;
  4. }

求双差后的协方差矩阵

ddcov(nb,b,Ri,Rj,nv,R);

3 ddcov( )

非对角线是Ri,对角线是Ri+Rj

  1. static void ddcov(const int *nb, int n, const double *Ri, const double *Rj,
  2. int nv, double *R)
  3. {
  4. int i,j,k=0,b;
  5. trace(3,"ddcov : n=%d\n",n);
  6. //nv*nv 矩阵嘛
  7. for (i=0;i<nv*nv;i++) R[i]=0.0;//初始化
  8. for (b=0;b<n;k+=nb[b++]) {
  9. for (i=0;i<nb[b];i++) for (j=0;j<nb[b];j++) {
  10. R[k+i+(k+j)*nv]=Ri[k+i]+(i==j?Rj[k+i]:0.0); //参考GNSS理论16.2讲。非对角线,参考卫星方差
  11. }
  12. }
  13. trace(5,"R=\n"); tracemat(5,R,nv,nv,8,6);
  14. }

4 validobs()

评价也不太懂,为什么可以不用管伪距

  1. i--基准站共识卫星索引下标
  2. f--频点 0--nf 是载波 nf--2*nf 是伪距
  3. 当f<nf,载波频率残差 有载波,就不管伪距了???解:这是f(0--2*nf)的情况,随着循环进行伪距最终会遍历到(nf--2*nf)检查
  4. y[f+i*nf*2]!=0 检查基准站共识卫星的载波残差
  5. f<nf||(y[f-nf+i*nf*2]!=0.0&&y[f-nf+j*nf*2]!=0.0) f-nf为真,短路原则
  6. 当f>=nf,伪距残差,载波也要有,不然伪距不可用(注释)
  7. y[f+i*nf*2] 每颗共识卫星伪距观测
  8. y[f-nf+i*nf*2] 载波观测
  9. */
  10. static int validobs(int i, int j, int f, int nf, double *y)
  11. {
  12. /* if no phase observable, psudorange is also unusable */
  13. return y[f+i*nf*2]!=0.0&&y[f+j*nf*2]!=0.0&&
  14. (f<nf||(y[f-nf+i*nf*2]!=0.0&&y[f-nf+j*nf*2]!=0.0));
  15. }

5 全部代码

  1. static int ddres(rtk_t *rtk, const nav_t *nav, double dt, const double *x,
  2. const double *P, const int *sat, double *y, double *e,
  3. double *azel, double *freq, const int *iu, const int *ir,
  4. int ns, double *v, double *H, double *R, int *vflg)
  5. {
  6. prcopt_t *opt=&rtk->opt;
  7. double bl,dr[3],posu[3],posr[3],didxi=0.0,didxj=0.0,*im;
  8. double *tropr,*tropu,*dtdxr,*dtdxu,*Ri,*Rj,freqi,freqj,*Hi=NULL;
  9. int i,j,k,m,f,nv=0,nb[NFREQ*4*2+2]={0},b=0,sysi,sysj,nf=NF(opt);
  10. trace(3,"ddres : dt=%.1f nx=%d ns=%d\n",dt,rtk->nx,ns);
  11. //计算基线长度
  12. bl=baseline(x,rtk->rb,dr);
  13. ecef2pos(x,posu); ecef2pos(rtk->rb,posr);
  14. /*
  15. Ri,Rj 基准站和流动站方差,计算以后协方差阵
  16. 对一些中间变量进行初始化,将双差伪距残差和双差载波相位残差初始化为0
  17. */
  18. Ri=mat(ns*nf*2+2,1); Rj=mat(ns*nf*2+2,1); im=mat(ns,1);
  19. tropu=mat(ns,1); tropr=mat(ns,1); dtdxu=mat(ns,3); dtdxr=mat(ns,3);
  20. /*
  21. double resp[NFREQ]; residuals of pseudorange (m)
  22. double resc[NFREQ]; residuals of carrier-phase (m)
  23. 将双差伪距残差和双差载波相位残差初始化为0
  24. 每个卫星每个频点的
  25. 为什么没有j<NFREQ*2,因为伪距和载波分为resp和resc了
  26. */
  27. for (i=0;i<MAXSAT;i++) for (j=0;j<NFREQ;j++) {
  28. rtk->ssat[i].resp[j]=rtk->ssat[i].resc[j]=0.0;
  29. }
  30. /*
  31. compute factors of ionospheric and tropospheric delay
  32. 如果卡尔曼滤波中包含对流层状态量,调用prectrop函数计算基站和移动站的对流层延迟湿分量,存在tropu[i]和tropur[i]中
  33. 如果卡尔曼滤波器包含电离层状态量,调用ionmapf函数分别计算基站和移动站处的投影函数
  34. */
  35. for (i=0;i<ns;i++) {
  36. if (opt->ionoopt>=IONOOPT_EST) {
  37. im[i]=(ionmapf(posu,azel+iu[i]*2)+ionmapf(posr,azel+ir[i]*2))/2.0;
  38. }
  39. if (opt->tropopt>=TROPOPT_EST) {
  40. tropu[i]=prectrop(rtk->sol.time,posu,0,azel+iu[i]*2,opt,x,dtdxu+i*3);
  41. tropr[i]=prectrop(rtk->sol.time,posr,1,azel+ir[i]*2,opt,x,dtdxr+i*3);
  42. }
  43. }
  44. /* 对各个系统以及载波相位、伪距分别进行循环处理
  45. 双重循环,遍历每个系统的每个频点(*2 包括伪距和载波)
  46. 哪种定位模式,是伪距,那我们从nf---2*nf,因为伪距残差是在y[]后半
  47. 是载波(动态定位模式及以后),那我们从0--2*nf,因为同时使用伪距和载波
  48. */
  49. for (m=0;m<6;m++) /* m=0:GPS/SBS,1:GLO,2:GAL,3:BDS,4:QZS,5:IRN */
  50. for (f=opt->mode>PMODE_DGPS?0:nf;f<nf*2;f++) {
  51. /*
  52. search reference satellite with highest elevation
  53. 选择最大的高度角卫星作为参考卫星
  54. 当然,我们也可以自己设置规则
  55. ssat_t ssat[MAXSAT]; satellite status
  56. 所有的共识卫星(全系统)每次都会被遍历一遍,我们选出所需的系统,然后再
  57. */
  58. for (i=-1,j=0;j<ns;j++) {
  59. sysi=rtk->ssat[sat[j]-1].sys;
  60. if (!test_sys(sysi,m)) continue;//当前共识卫星的系统和当前遍历的系统是否一样
  61. if (!validobs(iu[j],ir[j],f,nf,y)) continue; //评价不懂
  62. if (i<0||azel[1+iu[j]*2]>=azel[1+iu[i]*2]) i=j;//选择排序的那种,i最终是参考卫星的index
  63. }
  64. if (i<0) continue;
  65. /*
  66. make DD (double difference) 做双差,遍历共识卫星
  67. 一个历元一个历元处理,输入流动站和基准站的数据
  68. 流动站卫星残差 基准站卫星残差 每颗卫星每个频点的载波和伪距
  69. y[............... ................]
  70. 站间单差,两个接收机(a,b)观测同一颗卫星(p),作差(M1)。卫星p肯定是它们的共识卫星
  71. 星间双差,参考卫星站间单差(M2),M1-M2
  72. */
  73. for (j=0;j<ns;j++) {
  74. if (i==j) continue;
  75. sysi=rtk->ssat[sat[i]-1].sys;//参考卫星系统
  76. sysj=rtk->ssat[sat[j]-1].sys;//流动
  77. freqi=freq[f%nf+iu[i]*nf];
  78. freqj=freq[f%nf+iu[j]*nf];
  79. if (!test_sys(sysj,m)) continue;//判断和当前遍历系统一样否
  80. if (!validobs(iu[j],ir[j],f,nf,y)) continue;//检验数据有效yy[]
  81. //设计矩阵初始化,nx---待解量的个数
  82. if (H) {
  83. Hi=H+nv*rtk->nx;
  84. for (k=0;k<rtk->nx;k++) Hi[k]=0.0;
  85. }
  86. /*
  87. DD residual 作双差残差
  88. iu[i] 共识--参考卫星在基准站残差中的索引
  89. (y[f+iu[i]*nf*2]-y[f+ir[i]*nf*2]) 共识参考卫星站间单差
  90. */
  91. v[nv]=(y[f+iu[i]*nf*2]-y[f+ir[i]*nf*2])-
  92. (y[f+iu[j]*nf*2]-y[f+ir[j]*nf*2]);
  93. /* partial derivatives by rover position
  94. e[]是dx dy dz设计矩阵,大小:流动站观测到的卫星+基准站观测到的卫星。
  95. 因为在relpos()函数中,站间单差[先是算基准站,后是流动站],设计矩阵都填到e里面了。所以大小如此
  96. 设计矩阵=非参考卫星-参考卫星
  97. 因为我们求得是流动站位置,所以只用iu[]即可
  98. 至于不同频率,这里只有单频
  99. */
  100. if (H) {
  101. for (k=0;k<3;k++) {
  102. Hi[k]=-e[k+iu[i]*3]+e[k+iu[j]*3];
  103. }
  104. }
  105. /* DD ionospheric delay term */
  106. if (opt->ionoopt==IONOOPT_EST) {
  107. didxi=(f<nf?-1.0:1.0)*im[i]*SQR(FREQ1/freqi);
  108. didxj=(f<nf?-1.0:1.0)*im[j]*SQR(FREQ1/freqj);
  109. v[nv]-=didxi*x[II(sat[i],opt)]-didxj*x[II(sat[j],opt)];
  110. if (H) {
  111. Hi[II(sat[i],opt)]= didxi;
  112. Hi[II(sat[j],opt)]=-didxj;
  113. }
  114. }
  115. /* DD tropospheric delay term */
  116. if (opt->tropopt==TROPOPT_EST||opt->tropopt==TROPOPT_ESTG) {
  117. v[nv]-=(tropu[i]-tropu[j])-(tropr[i]-tropr[j]);
  118. for (k=0;k<(opt->tropopt<TROPOPT_ESTG?1:3);k++) {
  119. if (!H) continue;
  120. Hi[IT(0,opt)+k]= (dtdxu[k+i*3]-dtdxu[k+j*3]);
  121. Hi[IT(1,opt)+k]=-(dtdxr[k+i*3]-dtdxr[k+j*3]);
  122. }
  123. }
  124. /* DD phase-bias term 双差模糊度 */
  125. if (f<nf) {
  126. if (opt->ionoopt!=IONOOPT_IFLC) {
  127. v[nv]-=CLIGHT/freqi*x[IB(sat[i],f,opt)]-
  128. CLIGHT/freqj*x[IB(sat[j],f,opt)];
  129. if (H) {
  130. Hi[IB(sat[i],f,opt)]= CLIGHT/freqi;
  131. Hi[IB(sat[j],f,opt)]=-CLIGHT/freqj;
  132. }
  133. }
  134. else {
  135. //根据站星双差方差,还要减去(模糊度的derT),见理论。而在y[]的计算中,还没有考虑模糊度,在这里给补上
  136. //x[] 模糊度
  137. v[nv]-=x[IB(sat[i],f,opt)]-x[IB(sat[j],f,opt)];
  138. if (H) {
  139. Hi[IB(sat[i],f,opt)]= 1.0; //IB() 哪个卫星哪个频点的模糊度位置索引???不懂
  140. Hi[IB(sat[j],f,opt)]=-1.0;
  141. }
  142. }
  143. }
  144. //将计算好的残差存到ssat结构体
  145. if (f<nf) rtk->ssat[sat[j]-1].resc[f ]=v[nv];
  146. else rtk->ssat[sat[j]-1].resp[f-nf]=v[nv];
  147. /* test innovation 检查残差是否符合标准 maxinno */
  148. if (opt->maxinno>0.0&&fabs(v[nv])>opt->maxinno) {
  149. if (f<nf) {
  150. rtk->ssat[sat[i]-1].rejc[f]++;
  151. rtk->ssat[sat[j]-1].rejc[f]++;
  152. }
  153. errmsg(rtk,"outlier rejected (sat=%3d-%3d %s%d v=%.3f)\n",
  154. sat[i],sat[j],f<nf?"L":"P",f%nf+1,v[nv]);
  155. continue;
  156. }
  157. /* SD (single-differenced) measurement error variances 每颗卫星每个频率的观测值噪声(方差) */
  158. Ri[nv]=varerr(sat[i],sysi,azel[1+iu[i]*2],bl,dt,f,opt);
  159. Rj[nv]=varerr(sat[j],sysj,azel[1+iu[j]*2],bl,dt,f,opt);
  160. /* set valid data flags 设置标志,这个卫星这个频率有效 */
  161. if (opt->mode>PMODE_DGPS) {
  162. if (f<nf) rtk->ssat[sat[i]-1].vsat[f]=rtk->ssat[sat[j]-1].vsat[f]=1;
  163. }
  164. else {
  165. rtk->ssat[sat[i]-1].vsat[f-nf]=rtk->ssat[sat[j]-1].vsat[f-nf]=1;
  166. }
  167. trace(4,"sat=%3d-%3d %s%d v=%13.3f R=%8.6f %8.6f\n",sat[i],
  168. sat[j],f<nf?"L":"P",f%nf+1,v[nv],Ri[nv],Rj[nv]);
  169. //标记一下,双差的具体信息。哪两个卫星作差,载波还是伪距,第几个频率
  170. vflg[nv++]=(sat[i]<<16)|(sat[j]<<8)|((f<nf?0:1)<<4)|(f%nf);
  171. nb[b]++;//每个小频率有几颗卫星(残差有多少个)
  172. }
  173. b++;//每个小频率
  174. }
  175. /* end of system loop */
  176. /* baseline length constraint for moving baseline 基线约束,以后看 */
  177. if (opt->mode==PMODE_MOVEB&&constbl(rtk,x,P,v,H,Ri,Rj,nv)) {
  178. vflg[nv++]=3<<4;
  179. nb[b++]++;
  180. }
  181. if (H) {trace(5,"H=\n"); tracemat(5,H,rtk->nx,nv,7,4);}
  182. /* DD measurement error covariance 协方差矩阵 */
  183. ddcov(nb,b,Ri,Rj,nv,R);
  184. free(Ri); free(Rj); free(im);
  185. free(tropu); free(tropr); free(dtdxu); free(dtdxr);
  186. return nv;
  187. }

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/小蓝xlanll/article/detail/464546
推荐阅读
相关标签
  

闽ICP备14008679号