XiaoGongWei / MG_APP

MG-APP (https://link.springer.com/article/10.1007/s10291-020-00976-1);
https://xiaogongwei10.blog.csdn.net/
95 stars 52 forks source link

problem when the filter fail at an epoch #21

Open LucaFibbi opened 4 years ago

LucaFibbi commented 4 years ago

Thank you for the your attention to my issues. If you try to use the attach test case with the options: TropDeal=Hopfield(GPT2) Method=Kalman Cut-off Angle=0 Kinemati Yes Satellite system GPS PPP NoSmooth PPP Back Yes run PPP at the Epoch Number: 1844 of the backward filter are visible 5 satellites. The Kalman filter fail because one satellite has bad measure and it is not possible to calculate a solution with 4 satellites. The function m_KalmanClass.KalmanforStatic at row 662 of the file QPPPBackSmooth.cpp return false. So it is executed the following statement memset(spp_pos, 0, 3*sizeof(double)); which set to zero the current position.

file QPPPBackSmooth.cpp
function void QPPPBackSmooth::Run(bool isDisplayEveryEpoch)
659         if (!m_Solver_Method.compare("SRIF", Qt::CaseInsensitive))
                is_filter_good = m_SRIFAlgorithm.SRIFforStatic(prevEpochSatlitData,epochSatlitData,spp_pos,X,P);
            else
                is_filter_good = m_KalmanClass.KalmanforStatic(prevEpochSatlitData,epochSatlitData,spp_pos,X,P);
//Save the last epoch satellite data
            if(is_filter_good)
            {
                prevEpochSatlitData = epochSatlitData;
                continue_bad_epoch = 0;
            }
            else
            {
                continue_bad_epoch++;
                memset(spp_vct, 0, 3*sizeof(double));
                memset(spp_pos, 0, 3*sizeof(double)); <---- With an error in filter the position is set to zero, but it is again used in next epoch in the filter
                X.setZero();
675         }

From this moment the position remains equal to zero, until the end of data elaboration. This because the current position (zero) is used in next epoch and the Kalman filter produce wrong data that are delete at the quality control.

To avoid the problem I see the following possibilities:

LucaFibbi commented 4 years ago

Here some code changes that can prevent a situation such as the one reported in this issue from occurring. Are they correct ?

diff --git a/QPPPBackSmooth.cpp b/QPPPBackSmooth.cpp
--- a/QPPPBackSmooth.cpp
+++ b/QPPPBackSmooth.cpp
@@ -572,6 +572,7 @@ void QPPPBackSmooth::Run(bool isDisplayEveryEpoch)
     multepochSatlitData = last_allPPPSatlitData;// get all epoch data
     int all_epoch_len = multepochSatlitData.length();
     int continue_bad_epoch = 0;
+       double spp_vct[3];
 //Back smooth Multiple epoch cycles
         for (int epoch_num = all_epoch_len - 1; epoch_num >= 0;epoch_num--)
         {
@@ -582,9 +583,6 @@ void QPPPBackSmooth::Run(bool isDisplayEveryEpoch)
             if(epochSatlitData.length() != 0)
             {
                 epochTime = epochSatlitData.at(0).UTCTime;//Obtaining observation time(Each satellite in the epoch stores the observation time)
-                //Set the epoch of the satellite
-                for(int i = 0;i < epochSatlitData.length();i++)
-                    epochSatlitData[i].UTCTime.epochNum = epoch_num;
             }
             else
             {
@@ -603,16 +601,14 @@ void QPPPBackSmooth::Run(bool isDisplayEveryEpoch)
             {
                 int a = 0;
             }
-            // use spp compute postion and smooth pesudorange
-            double temp_spp_pos[3] = {last_allReciverPos.at(epoch_num).spp_pos[0],
-                                      last_allReciverPos.at(epoch_num).spp_pos[1],
-                                      last_allReciverPos.at(epoch_num).spp_pos[2]};
+            // use spp compute position and smooth pesudorange
+            memcpy(spp_vct,last_allReciverPos.at(epoch_num).spp_pos,3*sizeof(double));
             //Monitor satellite quality and cycle slip
             getGoodSatlite(prevEpochSatlitData,epochSatlitData, m_CutAngle);

             // The number of skipping satellites is less than m_minSatFlag
             // update at 2018.10.17 for less m_minSatFlag satellites at the begain observtion
-            if(epochSatlitData.length() < m_minSatFlag || temp_spp_pos[0] == 0)
+            if(epochSatlitData.length() < m_minSatFlag || spp_vct[0] == 0)
             {
                 if(m_isKinematic&&continue_bad_epoch++ > 8)
                 {
@@ -633,7 +629,7 @@ void QPPPBackSmooth::Run(bool isDisplayEveryEpoch)
                 autoScrollTextEdit(mp_QTextEditforDisplay, disPlayQTextEdit);// display for QTextEdit
                 // translation to ENU
                 VectorXd ENU_Vct;
-                double spp_vct[3] = {0};
+                memset(spp_vct, 0, 3*sizeof(double));
                 ENU_Vct.resize(32+3*epochSatlitData.length());
                 ENU_Vct.fill(0);
                 saveResult2Class(ENU_Vct, spp_vct, epochTime, epochSatlitData, epoch_num);
@@ -643,7 +639,6 @@ void QPPPBackSmooth::Run(bool isDisplayEveryEpoch)
             // Choose solve method Kalman or SRIF
             MatrixXd P;
             VectorXd X;//dX,dY,dZ,dT(Zenith tropospheric residual),dVt(Receiver clock),N1,N2...Nm(Ambiguity)[dx,dy,dz,dTrop,dClock,N1,N2,...Nn]
-            double spp_vct[3] = {0};// save pos before fillter
             bool is_filter_good = false;
             X.resize(5+epochSatlitData.length());
             X.setZero();
@@ -654,8 +649,6 @@ void QPPPBackSmooth::Run(bool isDisplayEveryEpoch)
                 P = last_fillter_Q;
                 X[0] = 0;X[1] = 0;X[2] = 0;
             }
-            // store spp position
-            spp_vct[0] = temp_spp_pos[0]; spp_vct[1] = temp_spp_pos[1]; spp_vct[2] = temp_spp_pos[2];
             if (!m_Solver_Method.compare("SRIF", Qt::CaseInsensitive))
                 is_filter_good = m_SRIFAlgorithm.SRIFforStatic(prevEpochSatlitData,epochSatlitData,spp_pos,X,P);
             else
@@ -698,6 +691,17 @@ void QPPPBackSmooth::Run(bool isDisplayEveryEpoch)
             //ENU_Vct[0] = P(1,1); ENU_Vct[1] = P(2,2); ENU_Vct[2] = P(3,3);
             ENU_Vct[0] = spp_pos[0]; ENU_Vct[1] = spp_pos[1]; ENU_Vct[2] = spp_pos[2];
             saveResult2Class(ENU_Vct, spp_vct, epochTime, epochSatlitData, epoch_num, &P);
+            if (spp_pos[0]==0.) {
+                // Recover last correct position
+                if (continue_bad_epoch<m_writeFileClass.allReciverPos.length()) {
+                    const RecivePos &correctRecivePos = m_writeFileClass.allReciverPos.at(continue_bad_epoch);
+                    spp_pos[0] =correctRecivePos.dX;
+                    spp_pos[1] =correctRecivePos.dY;
+                    spp_pos[2] =correctRecivePos.dZ;
+                } else {
+                    memcpy(spp_pos,last_allReciverPos.at(epoch_num).spp_pos,3*sizeof(double));
+                }
+            }
         }//End of multiple epochs    for (int n = 0; n < multepochSatlitData.length();n++)

         // clear multepochSatlitData
diff --git a/QPPPModel.cpp b/QPPPModel.cpp
--- a/QPPPModel.cpp
+++ b/QPPPModel.cpp
@@ -862,6 +862,17 @@ void QPPPModel::Run(bool isDisplayEveryEpoch)
             //ENU_Vct[0] = P(1,1); ENU_Vct[1] = P(2,2); ENU_Vct[2] = P(3,3);
             ENU_Vct[0] = spp_pos[0]; ENU_Vct[1] = spp_pos[1]; ENU_Vct[2] = spp_pos[2];
             saveResult2Class(ENU_Vct, spp_vct, epochTime, epochResultSatlitData, epoch_num, &P);
+            if (spp_pos[0]==0.) {
+                // Recover last correct position
+                if (continue_bad_epoch<m_writeFileClass.allReciverPos.length()) {
+                    const RecivePos &correctRecivePos = m_writeFileClass.allReciverPos.at(continue_bad_epoch);
+                    spp_pos[0] =correctRecivePos.dX;
+                    spp_pos[1] =correctRecivePos.dY;
+                    spp_pos[2] =correctRecivePos.dZ;
+                } else {
+                    memcpy(spp_pos, m_ApproxRecPos, 3*sizeof(double));
+                }
+            }
             epoch_num++;//Increase in epoch
         }//End of multiple epochs.  (int n = 0; n < multepochSatlitData.length();n++)
XiaoGongWei commented 4 years ago

dear LucaFibbi, This bug has been fixed. A new commit has been push. c5b639cab2a5c87c755688bfd468fa86f6c232f6 you can test again, please contact me if you have any questions :)