Characterizing thermal stability boundary is fundamental to the operation and control of power system, directly impacting its safety and economic efficiency. With the integration of large-scale renewable energy sources, operation modes of power system have become increasingly complex and variable. The influence characteristics of multi-dimensional variables on system boundary tend toward strong nonlinearity and coupling, posing challenges for traditional characterization methods. This paper first analytically derives the thermal stability boundary manifold and its gradient vector trajectory based on a single-machine equivalent system, exploring the nonlinear nature of system boundary manifold. Secondly, considering the active-reactive power coupling characteristics of renewable energy devices, it analyzes the nonlinear and non-convex features and patterns of system boundary manifold under typical scenarios such as weak sending and weak receiving conditions. Finally, a segmented approximation method suitable for high-dimensional nonlinear boundary manifold is proposed. This method utilizes sensitivity analysis and the K-Medoids algorithm to perform clustering analysis of high-dimensional gradient vectors and segmentation of boundary manifold, followed by fitting linear approximation expressions for each segment using the least squares method. The practicality and effectiveness of the proposed method are validated through a case study of the IEEE 3-machine 9-bus system.