As a relative positioning technique, light detection and ranging (LiDAR)-inertial odometry (LIO) is known to suffer from drifting and can only provide local coordinates. To compensate for these shortages of LIO, an effective way is to integrate a global navigation satellite system (GNSS) with LIO. In this contribution, we proposed a tightly coupled GNSS real-time kinematic (RTK)/inertial navigation system (INS)/LiDAR system under the factor graph optimization framework, termed FGO-GIL, to achieve high-precision and continuous navigation in urban environments. This integration system fuses raw GNSS measurements (i.e., pseudorange and carrier phase measurements) with inertial measurement unit (IMU) and LiDAR information at the observation level atop a factor graph. Moreover, a keyframe-based nonlinear optimization scheme is designed to efficiently utilize measurements from the mixed heterogeneous sensors. Specifically, nonkeyframes are united with IMU preintegration for interframe optimization, which can provide accurate and high-frequency state predictions for the scan-matching of keyframes. Sparse keyframes are used to construct LiDAR factors by matching with the submap for sliding window optimization. To evaluate the effectiveness of our approach, real-world experiments are conducted in both campus and urban environments. The results demonstrate that our system can achieve continuous decimeter-level positioning accuracy in these complex environments, outperforming other state-of-the-art frameworks.