A computer-implemented method, comprising: a. receiving, at one or more computer processors, a seismic dataset representative of a subsurface volume of interest and an initial earth model, the seismic dataset recorded with seismic sensors in a seismic survey of subterranean geological media; b. performing, via the one or more processors, full waveform inversion using b-splines to generate an updated earth model; c. performing, via the one or more processors, seismic imaging using the updated earth model to generate a seismic image; and d. displaying, via a user interface, the seismic image, characterized in that the performing the full waveform inversion using b-splines comprises: e. computing an objective function, J, according to J V = 1 2 ‖ B d B d * d s V − d 0 ‖ 2, wherein d s denotes synthetic data as a function of velocity V, d 0 denotes the observed data, and B d denotes the data domain b-spline operator B d = B M t, M r c 1, c 2, and wherein M t, M r is a set of b-spline dimension parameters and c1 and c2 are degrees of smoothness associated to each b-spline dimension; f. picking a set of b-spline dimension parameters Mt, M r so that B M t, M r c 1, c 2 d S is expected to match B M t, M r c 1, c 2 d 0 within a quarter scale-length given the initial earth model; g. picking a set of b-spline dimension parameters M z, M x in vertical and horizontal directions, respectively, so that the length-scales, L z / (M x -c1 + 1),L x / (M z -c2 + 1), are approximately the length-scales of velocity variations that match length-scales of B M t, M r c 1, c 2 d S from a forward simulation, where L x and L z are horizontal and vertical sizes of the model; h. solving b * = arg min b ‖ B d B d * B d B d * ( d s B m b − d 0 ‖ 2 given initial b-spline model parameters determined from solving arg min b || B m b - V init || 2, where d S is a function of a set of model space b-spline parameters b, wherein model space b-spline operator B m = B M z, M x c 1, c 2 and data domain b-spline operator B d = B M t, M r c 1, c 2 work in separate domains simultaneously; and i. repeating the steps of computing the objective function, picking the set of b-spline dimension parameters M t, M r. picking the set of b-spline dimension parameters M z, M x, and solving b ∗ with increased Mt, M r, M z, and M x so that matching between length-scales of the b-spline projected data and the model is maintained; j. computing a gradient with respect to velocity, ∇ V J, according to ∇ V J = ∂ d s V ∂ V * B d B d * B d B d * d s − d 0; and k. producing a quasi-Newton's update to generate the updated earth model. a. receiving, at one or more computer processors, a seismic dataset representative of a subsurface volume of interest and an initial earth model, the seismic dataset recorded with seismic sensors in a seismic survey of subterranean geological media; b. performing, via the one or more processors, full waveform inversion using b-splines to generate an updated earth model; c. performing, via the one or more processors, seismic imaging using the updated earth model to generate a seismic image; and d. displaying, via a user interface, the seismic image, e. computing an objective function, J, according to J V = 1 2 ‖ B d B d * d s V − d 0 ‖ 2, wherein d s denotes synthetic data as a function of velocity V, d 0 denotes the observed data, and B d denotes the data domain b-spline operator B d = B M t, M r c 1, c 2, and wherein M t, M r is a set of b-spline dimension parameters and c1 and c2 are degrees of smoothness associated to each b-spline dimension; f. picking a set of b-spline dimension parameters Mt, M r so that B M t, M r c 1, c 2 d S is expected to match B M t, M r c 1, c 2 d 0 within a quarter scale-length given the initial earth model; g. picking a set of b-spline dimension parameters M z, M x in vertical and horizontal directions, respectively, so that the length-scales, L z / (M x -c1 + 1),L x / (M z -c2 + 1), are approximately the length-scales of velocity variations that match length-scales of B M t, M r c 1, c 2 d S from a forward simulation, where L x and L z are horizontal and vertical sizes of the model; h. solving b * = arg min b ‖ B d B d * B d B d * ( d s B m b − d 0 ‖ 2 given initial b-spline model parameters determined from solving arg min b || B m b - V init || 2, where d S is a function of a set of model space b-spline parameters b, wherein model space b-spline operator B m = B M z, M x c 1, c 2 and data domain b-spline operator B d = B M t, M r c 1, c 2 work in separate domains simultaneously; and i. repeating the steps of computing the objective function, picking the set of b-spline dimension parameters M t, M r. picking the set of b-spline dimension parameters M z, M x, and solving b ∗ with increased Mt, M r, M z, and M x so that matching between length-scales of the b-spline projected data and the model is maintained; j. computing a gradient with respect to velocity, ∇ V J, according to ∇ V J = ∂ d s V ∂ V * B d B d * B d B d * d s − d 0; and k. producing a quasi-Newton's update to generate the updated earth model.