The cross Gramian is a useful tool in model order reduction but only applicable to square dynamical systems. Throughout this paper, time-limited cross Gramians is firstly extended to square bilinear systems that satisfies a generalized Sylvester equation, and then concepts from decentralized control are used to approximate a cross Gramian for non-square bilinear systems. In order to illustrate these cross Gramians, they are calculated efficiently based on shifted Legendre polynomials and applied to dimension reduction, which leads to a lower dimensional model by truncating the states that are associated with smaller approximate generalized Hankel singular values. In combination of the dominant subspace projection method, our reduction procedure is modified to produce a bounded-input bounded-output stable-preserved reduced model under some certain conditions. At last, the performance of numerical experiments indicates the validity of our reduction methods.