Gravimetric and magnetic surveys provide complementary information for subsurface characterization. However, their sparse, irregular sampling complicates accurate interpolation. Classical methods, such as kriging and cokriging, rely on predefined covariance models and fixed spatial neighborhoods, often failing to capture non-stationary structures and complex cross-variable interactions. We therefore propose a joint graph-based interpolation framework that learns end-to-end both the spatial graph topology and the interpolation function directly from the data. Survey locations are modeled as nodes in a learned graph with adaptive connectivity, enabling anisotropic and data-driven neighborhoods. Gravimetric and magnetic data are estimated jointly using a shared graph, allowing cross-field information transfer while preserving variable-specific behavior. Graph-based Dirichlet and biharmonic regularization enforce spatial smoothness and structural consistency and robust loss functions enhance stability under noisy observations. Experiments on field geophysical data demonstrate that the proposed method outperforms state-of-the-art interpolation techniques, particularly in data-limited scenarios.