Approximation of interacting kernels by sum of Gaussians (SOG) is frequently required in many applications of scientific and engineering computing in order to construct efficient algorithms for kernel summation or convolution problems. In this paper, we propose a kernel-independent SOG method by introducing the de la Vallee-Poussin sum and Chebyshev polynomials. The SOG works for general interacting kernels and the lower bound of Gaussian bandwidths is tunable and thus the Gaussians can be easily summed by fast Gaussian algorithms. The number of Gaussians can be further reduced via the model reduction based on the balanced truncation based on the square root method. Numerical results on the accuracy and model reduction efficiency show attractive performance of the proposed method.