A new diabatization method is presented, which is suitable for the development of accurate high-dimensional coupled potential energy surfaces for use in quantum dynamics studies. The method is based on the simultaneous use of adiabatic wave function and energy data, respectively, and combines block-diagonalization and diabatization by ansatz approaches. It thus is called hybrid diabatization. The adiabatic wave functions of suitable ab initio calculations are projected onto a diabatic state space and the resulting vectors are orthonormalized like in standard block-diagonalization. A parametrized diabatic model Hamiltonian is set up as an ansatz for which the block-diagonalization data can be utilized to find the optimal model. Finally, the parameters are optimized with respect to the ab initio reference data such that the deviations between adiabatic energies and eigenvalues of the model as well as projected state vectors and eigenvectors of the model are minimized. This approach is particularly advantageous for problems with a complicated electronic structure where the diabatic state space must be of higher dimension than the number of calculated adiabatic states. This is an efficient way to handle problems with intruder states, which are very common for reactive systems. The use of wave function information also increases the information content for each data point without additional cost, which is beneficial in handling the undersampling problem for high-dimensional systems. The new method and its performance are demonstrated by application to three prototypical systems, ozone (O3), methyl iodide (CH3I), and propargyl (H2CCCH).