During the exploration and development of deep coalbed methane (CBM), delineating the thickness of coal seam and lithofacies of the roof and floor is one of the major challenging tasks. In past attempts, the prediction methods of these parameters have been limited to the conventional inversion. However, the effect of coal shielding on adjacent reflecting layers restricts the identification of underlying sand effectively by conventional inversion. Also, the depth at which the deep CBM zone is located (1,500–2000 m) produces a significant overlap of P-wave impedance and Vp/Vs of sands and shale which increases classification uncertainty between these two lithofacies. We proposed a new workflow for high-precision quantitative seismic interpretation of deep CBM reservoir. Not only P-wave impedance but also GR is selected as the optimized attributes for lithofacies classification. To reduce the effect of strong reflection of coal seam and identifying thin coal layers, the seismic waveform indication inversion method is used to obtain high-resolution results of P-wave impedance and GR. It uses horizontal changes in seismic waveforms to reflect lithological assemblage characteristics for facies-controlled constraints. Then, Bayesian classification theory is used to achieve three-dimensional lithofacies classification with multi-source data. To improve the continuity and accuracy of the interpreted results, a Markov chain is applied in the Bayesian rule as the spatial prior constraint. A well-associated synthetic test and field data application in Ordos Basin demonstrates the accuracy of the proposed workflow. Compared with conventional inversion, the results of proposed workflow showed higher resolution and accuracy. By providing a new solution for the identification of roof and floor lithofacies of deep CBM reservoir, this workflow aims to contribute to the better exploration and development of deep CBM.