Seismic full waveform inversion (FWI) is a powerful method for estimating quantitative subsurface physical parameters from seismic data. As the FWI is a nonlinear problem, the linearized approach updates model iteratively from an initial model, which can get trapped in local minima. In the presence of a high-velocity contrast, such as at Moho, the reflection coefficient and recorded waveforms from wide-aperture seismic acquisition are extremely nonlinear around critical angles. The problem at the Moho is further complicated by the interference of lower crustal (Pg) and upper mantle (Pn) turning ray arrivals with the critically reflected Moho arrivals (PmP). In order to determine velocity structure near Moho, a nonlinear method should be used. We propose to solve this strong nonlinear FWI problem at Moho using a trans-dimensional Markov chain Monte Carlo (MCMC) method, where the earth model between lower crust and upper mantle is ideally parametrized with a 1-D assumption using a variable number of velocity interfaces. Different from common MCMC methods that require determining the number of unknown as a fixed prior before inversion, trans-dimensional MCMC allows the flexibility for an automatic estimation of both the model complexity (e.g. the number of velocity interfaces) and the velocity-depth structure from the data. We first test the algorithm on synthetic data using four representative Moho models and then apply to an ocean bottom seismometer (OBS) data from the Mid-Atlantic Ocean. A 2-D finite-difference solution of an acoustic wave equation is used for data simulation at each iteration of MCMC search, for taking into account the lateral heterogeneities in the upper crust, which is constrained from traveltime tomography and is kept unchanged during inversion; the 1-D model parametrization near Moho enables an efficient search of the trans-dimensional model space. Inversion results indicate that, with very little prior and the wide-aperture seismograms, the trans-dimensional FWI method is able to infer the posterior distribution of both the number of velocity interfaces and the velocity-depth model for a strong nonlinear problem, making the inversion a complete data-driven process. The distribution of interface matches the velocity discontinuities. We find that the Moho in the study area is a transition zone of.7 km, or a sharp boundary with velocities from around 7 km s-1 in the lower crust to 8 km s-1 of the upper mantle; both provide nearly identical waveform match for the field data. The ambiguity comes from the resolution limit of the band-limited seismic data and limited offset range for PmP arrivals.