In this study, we develop a 3D geostatistical seismic inversion method based on turning bands simulation and co-simulation to estimate acoustic impedance (AI) from seismic and well data. The proposed method uses an iterative approach based on cross-over concept in genetic algorithm optimization to perform global stochastic inversion. The objective function of the optimization algorithm is the measure of correlation coefficients between the modeled and observed seismic data. In the first iteration of the proposed algorithm, the seismic cubes corresponding to the AI realizations generated by the turning bands simulation are compared to the observed seismic data. Subsequently, the local areas of AI models producing the highest correlation coefficients to the observed seismic data are merged to construct a new supplementary cube that is further used as the secondary variable in the turning bands co-simulation for the next iteration. Generation of the seismic cubes and picking the best AI cube via cross-over genetic algorithm approach is performed iteratively until convergence to a constant correlation coefficient between the modeled and observed seismic data. The proposed method is applied to a synthetic dataset that includes 20 wells with known AI logs and 3D stacked seismic data. According to the results, the algorithm converges to a constant correlation coefficient after a few iterations. In addition, it is observed that employing multi-attribute analysis outputs (meta-attributes) during turning bands co-simulation in the initialization step would improve the final global correlation coefficient from 0.774 to 0.815.