Comment calculer numériquement une fonction à partir de son gradient bruyant?

J'ai le modèle .$\ s(x,y)=x^2+y^2, 0 \leq x \leq 1, 0 \leq y \leq 1$

Au lieu d'observer directement le modèle j'observe les dérivées du modèle + un peu de bruit (e):

$\ p(x,y)=s_x+e, q(x,y)=s_y+e$

A partir des mesures de p (x, y et q (x, y) je veux estimer s (x). Disons que je sais que s (0,0) = 0.

Selon le théorème du gradient: $\ s(x,y) = \int_{(0,0)}^{(x,y)} [s_x,s_y]d\vec{r}$

quel que soit le chemin que nous intégrons.

Comme petite expérience (dans Matlab), j'ai ajouté du bruit distribué normal, N (0,1), à p = 2x et q = 2y. Ensuite, j'ai intégré d'abord le long de x puis le long de y: SXY. Ensuite, j'ai intégré d'abord le long de y, puis le long de x: SYX.

Les résultats montrent que le théorème du gradient ne tient pas dans ce cas (à cause du bruit):

Les erreurs quadratiques moyennes par rapport au modèle sont les suivantes:

ErmsXY =
    0.1125
ErmsYX =
    0.0920

Comment puis-je trouver une meilleure estimation (moins d'erreur RMS et plus fluide) de s à partir de p et q?

ÉDITER:

D'après ce que j'ai lu; l'utilisation de l'intégrale de la courbe est appelée intégration locale. Il existe également des méthodes d'intégration globale où l'on essaie plutôt de choisir un S (x, y) qui minimise:

$\ \int_0^1\int_0^1 [|S_x - P|^2 + |S_y - Q|^2] dxdy$

Les méthodes d'intégration globale sont censées donner de meilleurs résultats lorsque le gradient est bruyant, mais comment faire cela en pratique?

EDIT 2:

Une approche que j'ai utilisée est la suivante:

nous introduisons d'abord des opérateurs de dérivation linéaire: .$\ s_x = D_x*s,s_y = D_y*s$

Le résultat est le système d'équations linéaires suivant:

$\ D_x*s=p, D_y*s=q$

Trouvez ensuite une solution Least Square Error à ces équations. Une solution LSE à ces équations est censée être équivalente à la minimisation de l'intégrale par le haut. Comment cela peut-il être démontré?

Les résultats sont bons:

L'erreur RMS représente environ 1/5 de celle de SXY et SYX et la solution est également plus fluide.

Cependant, cette approche présente certains inconvénients:

1. il est difficile à mettre en œuvre; doit utiliser les différences centrales et "aplatir" la matrice 2D en vecteur, etc.

2. Les matrices de dérivation sont très grandes et clairsemées, elles peuvent donc consommer beaucoup de RAM.

Une autre approche qui semble potentiellement à la fois plus simple à coder, moins consommatrice de RAM et plus rapide consiste à utiliser la FFT. Dans l'espace de Fourier, ces pdes deviennent une équation algébrique. Ceci est connu comme l'algorithme de Frankot-Chellappa, mais malheureusement je ne l'ai pas fait fonctionner sur mes données d'exemple.

Vous pouvez filtrer soit le gradient lui - même ou le résultat, . Il faudrait assez bien connaître les caractéristiques des véritables gradients pour savoir quelle est la largeur de bande de fréquence. À ce stade, vous pouvez concevoir un filtre passe-bas qui préservera le signal mais éliminera le bruit de fréquence plus élevée.$s$