Cardiac electrophysiology (EP) is a vital area of research driven by its critical clinical applications, particularly in diagnosing and treating life-threatening arrhythmias. The field relies heavily on biophysical models, such as the two-variable Aliev-Panfilov model, to capture the electrical behavior and pulse dynamics across the myocardium. However, a significant barrier in this domain is the computational intensity of these models; traditional numerical solvers are too slow for real-time clinical scenarios. While Machine Learning alternatives like Physics-Informed Neural Networks (PINNs) have been explored, they are limited by their inability to generalize, often requiring complete retraining if tissue parameters or initial conditions change. To overcome these limitations, this work introduces a novel surrogate modeling framework designed to learn the solution operator for the Aliev-Panfilov model. This method allows for generalization across varying tissue parameters and initial conditions without the need for retraining. The architecture utilizes a hybrid approach, combining Deep Operator Networks (DeepONets) with Convolutional Neural Networks (CNNs). In this system, the CNN acts as an encoder for the initial excitation state, while the DeepONet maps this encoded input and the relevant tissue parameters to the spatiotemporal evolution of the system. A key innovation of this research is the integration of contrastive loss mechanisms and beta-Variational Autoencoders (beta-VAE) to regularize and disentangle the latent space of the model. Dimensionality reduction analysis revealed that these mechanisms effectively group solutions into three distinct physical regimens: wave blockage, non-interactive propagation, and spiral wave generation. Furthermore, the study identifies direct correlations between the learned latent variables and essential clinical biomarkers, specifically APD90. The proposed surrogate model delivers a massive efficiency improvement, achieving a 363x performance gain by reducing dense-solution inference time from 49 seconds in a traditional numerical solver to just 0.135 seconds. This framework offers a highly efficient, interpretable foundation for large-scale parameter sweeps and real-time patient stratification, successfully bypassing the computational constraints of traditional solvers and the rigidity of single-instance neural approximations.