This is a difficult problem that I have been asking myself lately.
A naive approach is to use space filling or random sampling, filter the data so that you are looking at the ones with the lowest function values, and then generate clusters in that region of the performance space. This gives you distinct areas in the design space, of high performance, which can be further optimized. However, this approach doesn't help if you have local minima which are significantly worse than the global minima. It also doesn't work if the minima are in very small regions that will not be found by random sampling.
Another naive approach is find a "global optimum", save it, then penalize that area of the design space, so the next optimum is found on the next run.
Hopefully someone can present a "good" algorithm to directly detect how many modes the problem has. I'd be interested in the response.