I do have an Adalm Pluto SDR that can be used as an improvised spectrum analyzer up to 6GHz, but it is still not high enough I guess.
With a cut-off at 6GHz, you won't get estimates less than roughly 100ps, even if the true edge is faster.
The idea was to look at the spectrum of the pulses with the SA, and to write down the amplitude of each harmonic.
With 10MHz square you had to write down 300 amplitudes of odd harmonics up to 6Ghz - that's hard work if you cannot automate it

Eventually, starting from the amplitude of each harmonic, it should be possible to apply a reverse FFT in order to reconstruct the time domain pulses, and measure the edges raise/falling time on the reconstructed signal in the time domain (TD reconstruction being possible by assuming the phase - which is not measured by the SA - is always according to an input signal with perfect edges).
I tried that (using simulated ideal SA measurements), but without phase it is not as easy as you think. Eventually I got a quite nice reconstruction if the square wave has a rational duty cycle M/N, where M and N are positive integers, M < N, N >= 2, and N is not too large (where "not too large" depends on the slowest rise time you want to measure and the square wave frequency). With these rational duty cycles, I was able to extract a single impulse response from the time domain signal. But for some (arbitrary) duty cycle values, the impulse responses overlap, making separation in the time domain impossible. Smoothing with a window function helps to some extent (but increases the estimate rise time). Using a lower square-wave frequency should also reduce the set of duty cycles that don't work well (but at 1MHz you'd have to write down 3000 measurements and you'd need 75dB dynamic range). I have not yet investigated the sensitivity to noise or measurement errors.
EDIT:
Addedd example diagrams of my reconstruction attempts
1) when the duty cycle happens to be a "good" one
2) and for comparison, with a "bad" duty cycle
The simulated "ground truth" edge was shaped with a 3rd-order ellipic filter with 1Ghz cutoff.
[ The horizontal offset between ground truth and reconstructed edge seems to be a consequence of the extra cutoff at 6Ghz (since you can't provide measurements beyond 6Ghz). It becomes much smaller when I keep the frequency content >= 6Ghz. I don't consider it a big problem, as it does not affect the shape very much. ]
Btw, if you want to continue this discussion, I also suggest to open a separate thread, as it related, but not really on-topic either.