A memoryless time-invariant non-linearity, a waveshaper, has an input $W$ output $X$ relationship $X = f(W).$ To arbitrarily small error for a given range of input values, the function $f(x)$ can be thought to be a polynomial of sufficient degree. A polynomial $f(x)$ is a weighted sum of powers $x^n$ with nonnegative exponents $n$. The input $W$ to the waveshaper is a sinusoid, assumed to be of unit amplitude for simplicity of presentation:
$$W = \cos(\omega T + \alpha),$$
where $T$ represents time, $\omega$ represents frequency, and $\alpha$ represents phase at time $T = 0$. The output of the waveshaper is a weighted sum of powers of $W$ with exponents $0, 1, 2, 3, 4, 5, \ldots:$
$$1,\\
\color{blue}{\cos(\omega T + \alpha)},\\
\frac{1}{2}\cos(2ωT + 2α) + \frac{1}{2},\\
\frac{1}{4}\cos(3ωT + 3α) + \frac{3}{4}\color{blue}{\cos(ωT + α)},\\
\frac{1}{8}\cos(4ωT + 4α) + \frac{1}{2}\cos(2ωT + 2α) + \frac{3}{8},\\
\frac{1}{16}\cos(5ωT + 5α) + \frac{5}{16}\cos(3ωT + 3α) + \frac{5}{8}\color{blue}{\cos(ωT + α)},\\
\ldots$$
Highlighted in blue, the phase-shifted cosine of frequency $\omega$ appears always with the original phase shift in the different powers of $W.$ If the distortion is mild, then the power $W^1$ with a positive weight dominates the output $X$. Also each $n$th harmonic phase-shifted cosine of frequency $\omega$ appears in the powers of $W$ always with phase shift $n\alpha,$ so it will be in the output $X$ of the waveshaper with that phase shift or the inverted phase shift $n\alpha + \pi$ due to possible negative weights in the waveshaper polynomial.
Things come even easier to see if you construct $f(x)$ as a weighted sum of Chebyshev polynomials of the first kind $T_n(x):$
Figure 1. Chebyshev polynomials of the first kind. (source)
With $n = 0, 1, 2, 3, \ldots,$ the approximation error will be made arbitrarily small for range $-1 \le x \le 1$. The output will be constructed of basis functions that for the sinusoidal input are directly its harmonics:
$$T_n(W) = \cos(n\omega T + n\alpha).$$
If you have some parametric model of the unknown filter, perhaps you can find the parameters based on the phase relationships of the harmonics and construct an inverse filter that recovers the phase alignment of the harmonics and also adjusts their amplitudes so that the waveshaped signal $X$ is recovered.
If you don't have a parametric model for the filter, I don't know if knowing the phase frequency response at a set of discrete frequencies (and with a modulo $\pi$ ambiguity) helps to recover the magnitude frequency response using a relationship between the two for minimum-phase filters.
One parameterized model for a positive-side-soft-clipping waveshaper with sharpness $v > 0$ and maximum value $s$ is:
$$f(x) = x + s - \frac{\ln\left(e^{vx} + e^{sv}\right)}{v}.$$
).svg){kind=link}