Terahertz (THz) communication and extremely large-scale MIMO (XL-MIMO) are essential for achieving ultra-high data rates in future 6G systems. However, at sub-millimeter wavelengths, typical indoor materials exhibit significant roughness that invalidates conventional ideal smooth surface assumptions, while massive array apertures introduce pronounced near-field effects and spatial non-stationarity. To address these challenges, this paper proposes a hybrid near-field channel model utilizing surface scattering characteristics based on distinct measurement campaigns. First, based on typical indoor materials scattering measurements across the 260-400 GHz band, an improved Beckmann-Kirchhoff (B-K) model is developed to accurately characterize surface roughness and diffuse scattering behavior. The model independently analyzes single-bounce (SB) and multi-bounce (MB) clusters by applying deterministic rough surface scattering theory and geometry-statistical approach, respectively. Then, using near-field spatial non-stationarity measurements from a 630-element virtual array in the 330-360 GHz band, a Dual-Gaussian Mixture Model (DMM) and a Negative Binomial (NB) distribution are adopted to describe the lengths and the number of spatial visibility regions (VRs), respectively. Additionally, a Weibull distribution is employed to model the intra-region power fluctuations. Finally, comprehensive XL-MIMO channel evaluations within the same band demonstrate that the proposed model aligns closely with measured results in terms of the spatial cross-correlation function (SCCF), frequency cross-correlation function (FCF), and channel capacity. By reproducing the spatial sparsity of THz band, the proposed model overcomes the limitation of conventional standard models, such as 3GPP 38.901 and WINNER II, in significantly overestimating channel capacity.