Abstract: Most studies on the selection of mainshock and aftershock ground motions adopt the method of sequential matching to select mainshock and aftershock ground motions, which has difficulties accurately simulating the true correlation between the main and aftershock earthquake. To overcome this limitation, this study extended the basic theory of generalized conditional intensity measure (GCIM), proposing a method for constructing three-dimensional mainshock and aftershock GCIM conditional distributions conditioned on multiple IMs, and developed theoretical details of the joint selection method of three-dimensional mainshock and aftershock ground motions based on GCIM. Then, the main and aftershock earthquake information recorded at WAKC Station in New Zealand was used as the conditional IM to set the earthquake information, and 661 real MS–AS ground-motion pairs were used to calculate interevent and interevent residuals, which can construct correlation coefficient models of horizontal-horizontal mainshock and aftershock IMs and horizontal-vertical mainshock and aftershock IMs. Additionally, existing correlation coefficient models were used to construct three-dimensional mainshock and aftershock GCIM conditional distributions, and joint selection was carried out according to the target MS-AS ground motion database. The results show that the three-dimensional mainshock and aftershock GCIM method performs closer to the real ground motion response spectra than the ASK14 model; that there are some differences between the constructed three-dimensional mainshock and aftershock GCIM conditional distributions considering the correlations between mainshock and aftershock IMs and the "unconditional" distribution without considering the correlations between mainshock and aftershock IMs; and that the joint selection results of three-dimensional mainshock and aftershock ground motions match well with the target theoretical distribution.