So I shall answer myself, The main reasin is limited bandwidth of seismic data. So just one person discuss it

I think, the main reason for NMO stretch is the NMO correction itself. It is well explained in the online course  of the Delft University of Technology. Given a certain

stacking velocity and offset, the correction T-T0 becomes smaller when T0 becomes larger. Thus, the correction is not constant along a trace, even if we have a constant offset and constant velocity. Also, we can see from this correction that the effect will become more prominent when the offset becomes larger as well.

In addition, the NMO correction is non-linear with respect to the offset and zero-offset traveltime.

Thus, the NMO stretch increases at greater offset or shallower depth (smaller reflection time), see for example the book of Zhou" Practical Seismic Data Analysis" .

YES, NMO stretch  are due to the nonparallelism of the local traveltime of each reflection event or in the frequency domain to the nonphysical energy changes caused by the nonstationary time shifts applied (see Buchholtz, 1972 and dunkin and levin ,1973). The standard NMO correction got addtional issues that

add to the wavelet bandwidth compression, and those are the partial duplication

of the recorded events and the time inversion of the samples of a reflection. Several methods exist to mitigate and attenuate that effect...depending on the type fo seismic acquisition and/or processing history that affected the final stack.

Thanks David and Djamel, But do you prefer a non-stretch stacking process? i.e. stacking over a phase locked hyperbola?

you can get there by partial correction ...there is a nice paper on the subject :Nonstretch normal moveout through iterative partial

correction and deconvolution, geophysics, 79,131-141.Or , a nonstretch technique based on an iterative-matchingpursuit algorithm, which recovers wavelet distortions caused by NMO correction.