Global Navigation Satellite Systems (GNSS) have been greatly developed within recent decades, and they play essential roles in many civil and military applications. Since more usage of them has been taken in urban area, multipath effect then becomes one of the major error sources in positioning. The incident multipath rays are usually generated by the buildings surrounding the receiver, and cause severely auto-correlation function (ACF) distortion and undesired pseudo range measurement bias. The weak received signal strength also increases the difficulty in detecting and distinguishing line-of-sight (LOS) ray from the received signals. There are several researches attempting to reduce the influence of the multipath from different aspects. The nEML and MEDLL are wildly used methods to mitigate the multipath effect in time. The nEML decreases the ACF distortion impact by using narrow discriminator space, and the MEDLL wants to cancel the components from incident multipath rays in ACF by taking iterated parameters estimation for each ray. These two methods could reduce the bias in delay estimate into a low level but cannot eliminate it completely. A class of multi-correlator techniques has been also proposed in recent years. Compared with the nEML and MEDLL, with the use of the temporal ACF properties of different GNSS signals, the multi-correlator techniques could further mitigate the multipath effect, but they still suffer sort of estimation bias in pseudo range measurement. Besides, these multi-correlator techniques are very complicated and hardly used in real-time GNSS receivers. Multi-antenna technique is another promising one to mitigate the multipath effect. With the equipment of multi-antenna, GNSS receiver can not only enhance the signal strength by beamforming, but have the capability in acquiring signals’ spatial properties as well. By using some antennas to form a receiving array, the steering spatial beam is generated by suitable RF front-end designs or digital adaptive beamforming, and the rays from low elevation degrees could be greatly suppressed. Some signal properties, such as polarization and multi-band, are also considered in these antenna array designs. The spatial beam could also be used to form null steering to multipath directions, but it requires too many antennas and large array size. Meanwhile, few of these antenna techniques have taken signal temporal properties into considerations so that their multipath mitigation capabilities could be further improved. Some subspace-based methods with multi-antenna have been provided for the multipath environment applications. By taking advantages of signal characters in time and space, they could distinguish each received rays and estimate their spatial and temporal parameters jointly. Unlike the maximum-like methods, these subspace-based methods could also save much complexity. However, the existing methods does not cover the signal dynamic property cause by the doppler effect. The ACF location is always needed to be held during their estimations but the ACF waveforms are shifting due to the relative movements between the receiver and satellites. Therefore, the existing subspace-based methods are not suitable to be taken to process the real receiving GNSS signals. In this paper, we propose a novel subspace-based spatial and temporal multipath mitigation method, and modify the code tracking result to compensate the multipath bias. This method is cooperating with the delay lock loop (DLL) so that the CPS is followed with the unwanted bias but the subspace-based method can be applied without the code phase dispersion. Because of the BPSK modulation of the GPS signal, we are able to firstly construct a phase extended spatial and temporal snapshots based on the tracking results. By taking the receiving uniform linear array structure into account, a subspace rotation property can be built, and the DOAM-like correlation matrices are also generated. Additional forward and backward spatial smooth (FBSS) is then given on the correlation matrices to solve the coherent problem by the same navigation bits of the multipath rays. After the eigenvalue decomposition on the cross-correlation matrix, the estimates of the angle-of-arrivals (AOA) and the phased extended signatures of the incident rays can be acquired by the eigenvalues and corresponding eigenvectors, which are utilized to provide the relative delays and carrier phases estimation with respect to the DLL and PLL tracking results. The line-of-sight (LOS) ray’s parameters are picked up by its first arrival character, and we can finally modify the DLL tracking result and eliminate its bias. Besides, the application of multi-antenna strengthens the received signal power, which also improves the estimation accuracy. The multipath mitigation performances of the proposed method are evaluated in simulated multipath environments by maximum error envelop (MEE) and the root mean squared errors (RMSE).