In　this　paper,　we　study　efficient　numerical　schemes　of　the　classical　phase　field　elastic　bending　energy　model　that　has　been　widely　used　to　describe　the　shape　deformation　of　biological　lipid　vesicles,　in　which　the　free　energy　of　the　system　consists　of　an　elastic　bending　energy,　a　surface　area　constraint　and　a　volume　constraint.　One　major　challenge　in　solving　such　model　numerically　is　how　to　design　appropriate　temporal　discretizations　in　order　to　preserve　energy　stability　with　large　time　step　sizes　at　the　semi-discrete　level.　We　develop　a　first　order　and　a　second　order　time　stepping　scheme　for　this　highly　nonlinear　and　stiff　parabolic　PDE　system　based　on　the　＂Invariant　Energy　Quadratization＂　approach.　In　particular,　the　resulted　semi-discretizations　lead　to　linear　systems　in　space　with　symmetric　positive　definite　operators　at　each　time　step,　thus　can　be　efficiently　solved.　In　addition,　the　proposed　schemes　are　rigorously　proved　to　be　unconditionally　energy　stable.　Various　numerical　experiments　in　2D　and　3D　are　presented　to　demonstrate　the　stability　and　accuracy　of　the　proposed　schemes.　(C)　2016　Elsevier　B.V.　All　rights　reserved.