GPU accelerated lattice Boltzmann (LB) simulations of multiphase flow in porous media have become a powerful tool to study fluid displacement process in porous media. For porous structures with a very low porosity, indirect-addressing memory access methods are preferred due to significantly reduction of the memory footprint despite that those methods are more difficult to implement and the resulting performance may be more sensible to the computing architectures, such as cache hierarchy and size. The direct-addressing methods are straightforward to implement and are able to archive high throughput when the porosity is large, while the methods become less efficient when the structures are too sparse. Nevertheless, the direct-addressing methods combined with multi-block grid technique are promising for many applications. In this work, we present a hybrid-way to tackle multiphase flow simulations in porous media, where the direct-addressing method is employed for the main LB evolution while the indirect-addressing method is employed for the complex boundary conditions. The CSF-based LB color-gradient multiphase model and the geometry-based wetting model are employed to increase accuracy and stability. Thanks to the utilization of AA-pattern streaming scheme, non-slip boundary condition of arbitrary orientation can be enforced without extra cost. We perform comprehensive analysis of the computational performance of the present solver on NVIDIA GPUs. For typical sandstones, our results show that the present implementation is able to achieve over 1.5 speedup compared with other direct-addressing schemes on V100 and A100, respectively and, particularly, the computational performance of the boundary kernels is greatly increased thanks to the increased L2 cache size of the latest GPUs.