In this paper, we propose a model order reduction (MOR) method based on general orthogonal polynomials for K-power bilinear systems in the time domain. Constructing proper projection matrices by solving a series of linear equations, a reduced K-power bilinear system is produced, which preserves the original coupled structure. It can match several expansion coefficients of the original output. Then the error bound of our algorithm is also investigated. Moreover, the stability of the reduced system is discussed as well. Finally, two numerical examples are provided to illustrate the effectiveness of our algorithm.