Geoacoustic inversion can be a computationally expensive task in high-dimensional parameter spaces, typically requiring thousands of forward model evaluations to estimate the geoacoustic environment. We demonstrate Bayesian optimization (BO), an efficient global optimization method capable of estimating geoacoustic parameters in seven-dimensional space within 100 evaluations instead of thousands. BO iteratively searches parameter space for the global optimum of an objective function, defined in this study as the Bartlett power. Each step consists of fitting a Gaussian process surrogate model to observed data and then choosing a new point to evaluate using a heuristic acquisition function. The ideal acquisition function balances exploration of the parameter space in regions with high uncertainty with exploitation of high-performing regions. Three acquisition functions are evaluated: upper confidence bound, expected improvement (EI), and logarithmically transformed EI. BO is demonstrated for both simulated and experimental data from a shallow-water environment and rapidly estimates optimal parameters while yielding results comparable to differential evolution optimization.