Real-world datasets often contain missing or corrupted values. Completing multidimensional tensor-structured data with missing entries is essential for numerous applications. Smoothness-constrained low-rank factorization models have shown superior performance with reduced computational costs. While effective at capturing global and long-range correlations, these models struggle to reproduce short-scale, high-frequency variations in the data. In this paper, we introduce the \Generalized Least Squares Kernelized Tensor Factorization (GLSKF) framework for tensor completion. GLSKF integrates smoothness-constrained low-rank factorization with a locally correlated residual process; the resulting additive structure can effectively characterize both global dependencies and local variations. In particular, we define the covariance norm to enforce the smoothness of factor matrices in the global low-rank factorization, and use structured covariance/kernel functions to model the local processes. For model estimation, we develop an alternating least squares (ALS) procedure with closed-form solutions for each subproblem. To efficiently handle missing data, GLSKF utilizes projection matrices that preserve the Kronecker structure of covariances, facilitating fast computations through conjugate gradient (CG) and preconditioned conjugate gradient (PCG) algorithms. The proposed framework is evaluated on four real-world datasets across diverse tasks: traffic speed imputation, color image inpainting, video completion, and MRI image reconstruction. Experimental results confirm that GLSKF delivers superior effectiveness and scalability, establishing it as a robust solution for multidimensional tensor completion.