## 2 Methodology
### BSSD PPP/INS Tightly Coupled Integration  
The innovation vector of TCI is based on the difference between GNSS observations (pseudorange, carrier-phase, and doppler) and the corresponding values predicted by INS [20; 21; 28]. The state equation and observation equation can be expressed as  
\[X_{TCI,k}=\phi_{TCI,k-1}X_{TCI,k-1}+\mu_{TCI,k-1},\mu_{TCI,k-1}\sim(0,Q_{TCI,k}) \tag{11}\]  
\[Z_{TCI,k}=H_{TCI,k}X_{TCI,k}+\eta_{TCI,k},\eta_{TCI,k}\sim N(0,R_{TCI}) \tag{12}\]  
\[Z_{TCI,k}= \begin{bmatrix}P_{GNSS,IF}-P_{INS,IF}\\ L_{GNSS,IF}-L_{INS,IF}\\ D_{GNSS,IF}-D_{INS,IF}\end{bmatrix} = \begin{bmatrix}Z_{P_{IF}}\\ Z_{L_{IF}}\\ Z_{D_{IF}}\end{bmatrix}  \tag{13}\]where \(\phi_{TCI,kk-1}\) is the system transform matrix from epoch \(k-1\) to epoch \(k\); \(\mu_{TCI,k-1}\) represent the state noise with the covariance of \(Q_{TCI,k}\); \(Z_{P_{IF}}\), \(Z_{L_{IF}}\), and \(Z_{D_{IF}}\) represent the innovation vector of pseudorange, carrier-phase, and doppler, respectively; \(P_{INS,IF}\), \(L_{INS,IF}\), and \(D_{INS,IF}\) are the INS-predicted values; \(P_{GNSS,IF}\), \(L_{GNSS,IF}\), \(D_{GNSS,IF}\) are the GNSS measurement; \(\eta_{TCI,k}\) represents observation noise with the prior covariance of \(R_{TCI}\). In order to obtain the INS-predicted values corresponding to GNSS measurements, the position and velocity of the receiver updated by INS mechanization are required. However, the reference centers of the INS and GNSS antenna are different from each other, which results in a lever-arm system offset. Therefore, the linearization functions \(Z_{TCI,k}\) after considering the lever arm can be written as [22; 29]  
\[\delta Z_{TCI,k}= \begin{bmatrix}\delta Z_{P_{IF}}\\ \delta Z_{L_{IF}}\\ \delta Z_{D_{IF}}\end{bmatrix} = \begin{bmatrix}C_{1}\Big{(}\delta p_{INS}^{n}+ \Big{(}C_{b}^{n}l^{b}\times\Big{)}\delta\psi\Big{)}+c\delta t_{r}-m_{wet}d_{wet }\\ C_{1}\Big{(}\delta p_{INS}^{n}+\Big{(}C_{b}^{n}l^{b}\times\Big{)}\delta\psi \Big{)}+c\delta t_{r}-m_{wet}d_{wet}+\delta N_{IF}\lambda_{IF}\\ C_{2}D^{-1}\delta p_{INS}^{n}+C_{n}^{*}H_{\psi}\delta\psi+C_{n}^{*}\delta v_{INS }^{n}+C_{n}^{*}C_{b}^{n}\Big{(}l^{b}\times\Big{)}\delta\omega_{lib}^{b}+c\delta t _{r}\end{bmatrix}  \tag{14}\]  
\[D^{-1}= \begin{bmatrix}1/(R_{M}+h)&0&0\\ 0&1/(R_{N}+h)\cos(B)&0\\ 0&0&-1\end{bmatrix}  \tag{15}\]  
where \(\delta Z_{P_{IF}}\), \(\delta Z_{L_{IF}}\), and \(\delta Z_{D_{IF}}\) represent the differential form of innovation vector; \(n\), \(b\), and \(i\) are the navigation frame (\(n\)), the body frame (\(b\)), and the inertial frame (\(i\)); \(C_{k}^{j}\) (\(k=n\) and \(b\), \(j=e\) and \(n\)) represent the rotation matrix from the \(k\)-frame to the \(j\)-frame; \(l^{b}\) is the lever-arm; \(\delta p_{INS}^{n}\), \(\delta v_{INS}^{n}\), and \(\delta\psi\) are the corrections of position, velocity, and attitude at IMU center in \(n\)-frame; \(\delta t_{r}\) and \(\delta t_{r}\) represent the corrections of receiver clock offset and drift; \(C_{1}\) is the rotation matrix of position corrections from \(e\)-frame to \(n\)-frame; \(C_{2}\) is the differential form of \(C_{n}^{*}\); \(\delta\omega_{lib}^{b}\) means gyroscope errors [22; 29]. Based on Equation (12), the satellite single-difference matrix (SSDM) can be expressed as [28]  
\[SSDM= \begin{bmatrix}SSDM_{G}&0\\ 0&SSDM_{B}\end{bmatrix} _{3h\times 2m} \tag{16}\]  
\[SSDM_{B}= \begin{bmatrix}1&0&\cdots&-1&0&0&\cdots&0\\ 0&\ddots&0&-1&\vdots&\vdots&0&\vdots\\ \vdots&0&1&-1&\vdots&\vdots&\vdots&\vdots\\ \vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\ \vdots&\vdots&0&-1&1&0&\cdots&\vdots\\ \vdots&\vdots&\vdots&-1&\vdots&\ddots&0&\vdots\\ \vdots&0&\vdots&-1&0&\cdots&1&0\\ 0&\cdots&0&-1&0&0&\cdots&1\end{bmatrix} _{h_{2}\times m} \tag{17}\]  
where \(SSDM_{G}\) is obtained referring to \(SSDM_{B}\); \(h_{1}\) and \(h_{2}\) (\(h=h_{1}+h_{2}\)) are the total number of GPS and BDS available satellites; \(m\) represents the number of estimated parameters for the BSSD model.