Accurate and globally referenced positioning is crucial to autonomous systems with navigation requirements, such as unmanned aerial vehicles (UAV) and autonomous driving vehicles (ADV). GNSS/LiDAR integration is a popular sensor pair that can provide outstanding positioning performance in open areas. However, the accuracy is significantly degraded in urban canyons, due to the excessive unmodeled non-Gaussian GNSS outliers caused by multipath effects and none-line-of-sight (NLOS) receptions. As a result, the violation of the Gaussian assumption can severely distort the sensor fusion process, such as the extended Kalman filter (EKF). To mitigate the effects of these non-Gaussian GNSS outliers, this paper proposes to leverage the Gaussian mixture model (GMM) to describe the potential noise of GNSS positioning and apply it to further sensor fusion. Instead of relying on excessive offline parameterization and tuning, the parameters of the GMM are estimated simultaneously based on the residuals of the GNSS measurements using an expectation-maximization (EM) algorithm. Then the state-of-the-art factor graph optimization (FGO) is applied to integrate the GNSS positioning and LiDAR odometry based on the estimated GMM. The experiment in a typical urban canyon is conducted to validate the performance of the proposed method. The result shows that the GMM can effectively mitigate the effects of GNSS outliers and improves positioning performance.