Abstract: The reasonable selection of GNSS/leveling points is very important to GNSS height fitting, step-by-step elimination method is a good method for the optimization selection of joint-observation points. The traditional elimination method selects the GNSS/leveling points based on the minimum fitting error of height anomaly, which will easily lead to uneven joint-observation points. According to this view, this paper proposes to optimize the joint-observation point based on the area size of Thiessen polygons generated by the GNSS/leveling points, and on this basis to improve the traditional method, namely taking into account both the size of the height anomaly fitting error and the area size of polygons generated by Thiessen method (referred to as the synthesis method). The experimental results show that the synthesis method can improve the distribution of joint-observation points, and get the fitting results of height anomaly with high stability and accuracy.