Recently, the factor graph-based global navigation satellite systems (GNSS) positioning methods have attracted much attention for more robust positioning performance in urban canyons compared with the traditional filter-based method. In the existing factor graph optimization (FGO)-based GNSS positioning methods, the pseudorange and Doppler measurements are mainly utilized to construct factors among consecutive GNSS epochs for outlier resisting and robust localization. However, the potential of high-precision positioning by using carrier phase observations is not fully explored. A factor graph optimization-based multi-GNSS real-time kinematic (FGO-RTK) framework is proposed to fill this gap, aiming to realize robust and precise positioning in urban canyons. In our method, a sliding window-based FGO estimator is designed, in which the continuously tracked double difference ambiguities are used to establish the ambiguity constraints of position states that have common-view satellites within the window. A marginalization-based carrier phase ambiguity propagation (AP) method is developed to fully use the information of carrier phase observations for achieving more reliable and continuous ambiguity resolution. Experiments in both medium urban and deep urban environments verified the effectiveness of the proposed method. Results show that FGO-RTK (with or without the AP) can achieve centimeter-level positioning accuracy in medium urban situations, showing comparable performance with the traditional extended Kalman filter (EKF)-based RTK. With the GNSS observation conditions deteriorating in deep urban environments, FGO-RTK without AP cannot outperform EKF-RTK. By contrast, the proposed FGO-RTK, considering the AP, achieves significant improvements in centimeter- to decimeter-level positioning availability, and the 3D positioning accuracy is improved by 69.6%, compared with the EKF-RTK. Furthermore, by analyzing the impact of the window size against the performance of FGO-RTK, we found that in comparison to the FGO-RTK only using the carrier phase observations inside a fixed-size window, the proposed AP method can significantly reduce the dependence on window size, and the optimal window size can be reduced from 10 to 4 epochs with more than 50% optimization time decreasing.