RIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B.P. 5 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

<figure><img src="image_1.png" /></figure>

NOUVELLES APPLICATIONS
DE LA MODÉLISATION
DES SUITES GAUSSIENNES
A DEUX INDICES

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

S. Attasi, J. F. Abramatic

Rapport de Recherche N° 157

Mars 1976
NOUVELLES APPLICATIONS DE LA MODÉLISATION
DES SUITES GAUSSIENNES A DEUX INDICES

S. Attasi, J.F. Abramatic

Résumé :

La modélisation stochastique des suites gaussiennes à deux indices avait permis de traiter les pro-
blèmes de filtrage statistique, «ligne par ligne» (Attasi). Un élément nouveau, la possibilité de dé-
terminer l'innovation «point par point» (dans le sens d'un balayage TV) ouvre la voie à l'exten-
sion, au traitement d'images, de nombreux algorithmes récursifs optimaux de traitement statisti-
que des suites temporelles.

Abstract :

In preceeding publication (Attasi), stochastic models for double indexed sequences were introdu-
ced, and «line by line» recursive estimation algorithms were derived. A new result, the possibili-
ty of determining recursively the «innovation» of the scan signal opens the way to various ex-
tensions of classical recursive algorithms used in time series analysis.
TABLE DES MATIERES

Pages

NOUVELLES APPLICATIONS DE LA MODELISATION STOCHASTIQUES DES SUITES A DEUX INDICES.................................................. 1

I. RAPPELS [ATTASI].................................................. 1

I.1. Définition d'une classe de suites gaussiennes.................. 1

I.2. Filtrage statistique aux moindres carrés (dans la suite p = 1).................................................. 2

II. CALCUL DE L'INNOVATION DU SIGNAL D'IMAGE ET APPLICATION AU CALCUL DE PROBABILITES A POSTERIORI.................................. 6

III. APPLICATIONS PRATIQUES........................................ 10

III.1. Détection de signaux........................................ 10

III.2. Codage optimal et application à la compression.................. 13

IV. ESTIMATION DES PARAMETRES DU "FILTRE BLANCHISSANT" PAR MAXIMUM DE VRAISEMBLANCE.................................................. 15

IV.1. Cas général (n > 1 ou R > 0).................................. 16

IV.2. Cas particulier, n=1, R=0.................................... 19

CONCLUSION................................................................. 22

REFERENCES................................................................. 23
# NOUVELLES APPLICATIONS DE LA MODELISATION

## STOCHASTIQUES DES SUITES A DEUX INDICES

### I. RAPPELS [ATTASI]

#### I. 1. Définition d'une classe de suites gaussiennes.

Généralisant la notion de suites gaussiennes à représentation markovienne [FAURRE, KALMAN, KAILATH], une classe de suites gaussiennes (notée Σ), {yᵢ,ⱼ; (i,j) ∈ Z²} à valeurs dans ℝP a été étudiée. Cette classe est constituée des suites admettant des modèles du type :

$$ (1) \qquad \left\{ \begin{array}{l} y_{i,j} = H x_{i,j} \\ x_{i,j} = F_1 x_{i-1,j} + F_2 x_{i,j-1} - F_1 F_2 x_{i-1,j-1} + v_{i,j} \end{array} \right. $$

où

i) {vᵢ,ⱼ} est une suite de variables gaussiennes centrées indépendantes, à valeurs dans ℝⁿ, avec E[vᵢ,ⱼ v'ᵢ,ⱼ] = Q (+).

ii) H, F₁, F₂ sont des matrices respectivement (p × n), (n × n), (n × n) vérifiant F₁F₂ = F₂F₁.

La généralité et "l'identifiabilité" de cette classe Σ ont été démontrées dans les termes suivants :

- Toute suite obtenue par des transformations linéaires, (non nécessairement "causales") sur un bruit blanc peut être approchée par une suite de la classe Σ.

- A partir des données statistiques (Λ(k,1) = E[yᵢ+k,j+1 y'ᵢ,j]) d'une suite {yᵢ,ⱼ} de la classe Σ, on sait déterminer l'ensemble des quadruplets

(+) Pour toute matrice M, M' désigne la transposée de M.
`{H,F₁,F₂,Q} de dimension minimale, associés à cette suite. Cet ensemble est
de plus réduit à un seul élément, modulo la relation d'équivalence :
{H,F₁,F₂,Q} ⇔ {HT⁻¹, TF₁T⁻¹, TF₂T⁻¹, TQT'} où T est une matrice de changement
de base dans Rⁿ.

Remarque :

Dans le cas abordé plus loin, de suites de la classe Σ perturbées par
un bruit blanc additif, le problème de l'identification à partir de la fonction
de corrélation a également une solution unique [ATTASI].

I. 2. Filtrage statistique aux moindres carrés (dans la suite p = 1).

On considère des suites de la classe Σ, perturbées par un bruit
additif, suivant le modèle :

$$
(2) \quad \left\{
\begin{array}{l}
y_{i,j} = H x_{i,j} + w_{i,j} \\
x_{i,j} = F_1 x_{i-1,j} + F_1 x_{i,j-1} - F_1 F_2 x_{i-1,j-1} + v_{i,j}
\end{array}
\right.
$$

avec

$$
E[w_{i+k,j+1} w'_{i,j}] = R \delta_{k,l} , R > 0
$$

$$
(3) \qquad E[v_{i,j} w'_{k,l}] = 0.
$$

A partir de ce modèle on sait déterminer des estimées aux moindres
carrés, dont la plus importante pour la suite est :

$$
(5) \quad E[y_{i,j} | y_{k,l}; 0 \le k < i, l \in \mathbb{Z}] = H E[x_{i,j} | y_{k,l}; 0 \le k < i, l \in \mathbb{Z}].
$$

Posons :

$$
\bullet \quad \hat{x}_{i,j} = E[x_{i,j} | y_{k,l}; 0 \le k < i, l \in \mathbb{Z}]
$$
• $\tilde{\hat{x}}_{i,j} = x_{i,j} - \hat{x}_{i,j}$ (l'erreur d'estimation au point $(i,j)$).

• $S_i(z) = \text{le spectre associé à la suite } \{\tilde{x}_{i,j} ; j \in Z\}$.

• $I_{i,j} = y_{i,j} - H \hat{x}_{i,j} = y_{i,j} - E[y_{i,j}/y_{k,l} ; 0 \le k < i, l \in Z]$.

Ainsi $\{I_{i,j} ; j \in Z\}$ représente "l'innovation" de la suite $\{y_{i,j} ; j \in Z\}$ par rapport aux "observations" $\{y_{k,l} ; 0 \le k < i, l \in Z\}$.

• $\{\hat{\tilde{x}}_{i,j} ; j \in Z\}$ la solution du problème de lissage à un indice :

$$ (6) \qquad \hat{\tilde{x}}_{i,j} = E[\tilde{x}_{i,j}/I_{i,1} ; 1 \in Z] = E[\tilde{x}_{i,j}/H \tilde{x}_{i,1} + w_{i,1}, 1 \in Z]. $$

Alors, les équations donnant $\hat{x}_{i,j}$ sont :

$$ \left\{ \begin{array}{l} \hat{x}_{0,j} = 0 \\[1em] I_{i,j} = y_{i,j} - H \hat{x}_{i,j} \\[1em] \hat{x}_{i,j} = F_1 (\hat{x}_{i-1,j} + \hat{\tilde{x}}_{i-1,j}) \\[1em] \hat{\tilde{x}}_{i,j} = H_i \eta_{i,j} + (G_i' - H_i P_i) \gamma_{i,j} \\[1em] \eta_{i,j} = F_i \eta_{i,j-1} + T_i (I_{i,j-1} - H H_i \eta_{i,j-1}) \\[1em] \gamma_{i,j} = (F_i - T_i H H_i)' \gamma_{i,j+1} + H_i' H' R_i^{-1} (I_{i,j} - H H_i \eta_{i,j}) \end{array} \right. \quad (7) $$

où

• $(H_i, F_i, G_i)$ est une "réalisation" du spectre $S_i(z)$ :

$$ (8) \qquad S_i(z) = H_i (zI - F_i)^{-1} G_i + G_i' (z^{-1} I - F_i') H_i' + H_i G_i $$
lequel est donné par l'équation récurrente :

$$
(9) \quad \left\{
\begin{aligned}
S_{i+1}(z) &= F_1(S_i(z) - S_i(z)H'(HS_i(z)H'+R)^{-1}HS_i(z))F_i' + Q(z) \\
Q(z) &= (zI - F_2)^{-1}Q(z^{-1}I - F_2')^{-1} \\
S_O(z) - F_1S_O(z)F_1' &= Q(z)
\end{aligned}
\right.
$$

• Les matrices $P_i, T_i, R_i$ sont liées à la factorisation du spectre $(HS_i(z)H'+R)$ associé à la suite $\{I_{i,j} ; j \in Z\}$, et sont données par :

$P_i$ est la limite de l'équation de RICCATI :

$$
(10) \quad \left\{
\begin{aligned}
P_{i,0} &= 0 \\
P_{i,1+1} &= F_i P_{i,1} F_i' + (F_i G_i H' - F_i P_i H_i' H') (HH_i G_i H' + R - HH_i P_i H_i')^{-1} \\
&\qquad (F_i G_i H' - F_i P_i H_i' H')'
\end{aligned}
\right.
$$

$$
(11) \quad R_i = HH_i (G_i - P_i H_i') H' + R
$$

$$
(12) \quad S_i = F_i (G_i - P_i H_i') H'
$$

$$
(13) \quad T_i = S_i R_i^{-1}.
$$

Remarques importantes :

1 - Cas où R = 0 : prédiction ligne par ligne sans bruit d'observation. On a toujours :

$$
(14) \quad \hat{x}_{i+1,j} = F_1 (\hat{x}_{i,j} + \hat{\bar{x}}_{i,j})
$$
où $\hat{x}_{i,j}$ est défini par l'équation (6) (avec $w_{i,j} = 0$).

Pour déterminer $\hat{x}_{i,j}$ il faut donc distinguer deux cas :

i) $n > 1$ : le problème de lissage à un indice défini par l'équation (6) garde tout son sens ; $\hat{x}_{i,j}$ se détermine par les mêmes équations que le cas régulier ($R > 0$), en faisant $R = 0$.

ii) $n = 1$ : dans ce cas :

$$ (15) \quad \left\{ \begin{array}{l} y_{i,j} = x_{i,j} \\ \hat{x}_{i,j} = \tilde{x}_{i,j} = I_{i,j} \end{array} \right. $$

et compte-tenu de (14), on a :

$$ (16) \quad \hat{x}_{i,j} = F_1 y_{i-1,j} $$

2 - Si l'on définit :

$$ (17) \quad L_i \triangleq G_i' - H_i P_i $$

alors la posivité du spectre $S_i(z)$ entraîne que :

$$ (18) \quad \left\{ \begin{array}{l} L_i = H_i \Sigma_i \\ \Sigma_i - F_i \Sigma_i F_i' \ge 0 \end{array} \right. $$

($\Sigma_i$ s'interprète en effet comme une covariance d'erreur de filtrage).
## II. CALCUL DE L'INNOVATION DU SIGNAL D'IMAGE ET APPLICATION AU CALCUL DE PROBABILITES A POSTERIORI

En se rappelant que les méthodes de transmission et d'acquisition couramment utilisées, transforment l'image y(i,j) en un signal temporel s(k) dit de "balayage horizontal ", il paraît naturel de définir l'innovation $v_{i,j}$ d'un point y(i,j) de l'image comme l'innovation du point s(k) du signal temporel correspondant, soit :

$$ (19) \qquad v_{i,j} = y_{i,j} - E[y_{i,j} | y_{k,l}; k < i, l \in \mathbb{Z} \text{ et } y_{i,l} ; l < j]. $$

### THEOREME.

(Cas général $R > 0$ ou $n > 1$).

La suite $\{v_{i,j}\}$ définie par l'équation (19) est donnée par :

$$ (20) \qquad v_{i,j} = y_{i,j} - H\hat{x}_{i,j} - HH_i \eta_{i,j} = I_{i,j} - E[I_{i,j} | I_{i,1} ; 1 < j] $$

où $\hat{x}_{i,j}, I_{i,j}, \eta_{i,j}$ sont donnés par les équations (6).

### Démonstration.

Le passage de l'équation (19) à l'équation (20) suit les étapes suivantes :

$$ (21) \qquad \begin{aligned} v_{i,j} = y_{i,j} - & E[y_{i,j} | y_{k,l} ; k < i, l \in \mathbb{Z}] \\ & - E[y_{i,j} | \{y_{i,1} ; 1 < j\} - E[y_{i,1} ; 1 < j | y_{k,l} ; k < i, l \in \mathbb{Z}]] \end{aligned} $$

Soit :

$$ (22) \qquad v_{i,j} = y_{i,j} - H\hat{x}_{i,j} - E[y_{i,j} | I_{i,1} ; 1 \in \mathbb{Z}] $$

$$ (23) \qquad v_{i,j} = y_{i,j} - H\hat{x}_{i,j} - E[I_{i,j} | I_{i,1} ; 1 \in \mathbb{Z}] $$

$$ (24) \qquad v_{i,j} = y_{i,j} - H\hat{x}_{i,j} - HH_i \eta_{i,j}. \quad \blacksquare $$
Dans la suite, on ne s'intéressera qu'aux solutions limites $H_\infty$, $F_\infty$, $P_\infty$, $G_\infty$, $R_\infty$, $S_\infty(z)$ ... Dans le cas de $H_\infty$, $F_\infty$, $G_\infty$, l'équation (8) ne pourra être qu'une approximation du spectre $S_\infty(z)$ puisque celui-ci n'est pas rationnel.

Cas particulier $n = 1, R = 0$ :

L'innovation $\{v_{i,j}\}$ est toujours donnée par (20) avec :

$$ (25) \qquad \left\{ \begin{array}{l} H = H_\infty = 1 \\[1em] \hat{x}_{i,j} = F_1 y_{i-1,j} \\[1em] \eta_{i,j} = F_2 I_{i,j-1} = F_2 (y_{i,j-1} - \hat{x}_{i,j-1}) \end{array} \right. $$

autrement dit :

$$ (26) \qquad \left\{ \begin{array}{l} v_{i,j} = y_{i,j} - F_1 y_{i-1,j} - F_2 y_{i,j-1} + F_1 F_2 y_{i-1,j-1} \\[1em] E[v_{i,j}^2] = Q. \end{array} \right. $$

On sait donc déterminer par un calcul récurrent l'innovation d'une suite quelconque de la classe $\Sigma$, bruitée ($R > 0$) ou non ($R = 0$).

Ce résultat permet déjà de déterminer la probabilité d'occurence de l'échantillon :

$$ (27) \qquad Y_{N,M} = \{y_{i,j} ; 0 \le i \le N, 0 \le j \le M\}, $$

sous l'hypothèse qu'il représente une réalisation partielle d'une suite gaus- sienne à deux indices, de la classe $\Sigma$, modélisée par l'ensemble de paramètres $\Theta = \{H, F_1, F_2, Q, R\}$ donné.
Posons $P(Y_{N,M}/\theta)$ cette probabilité, et soient :

$$ (28) \qquad \left\{ \begin{array}{l} \bar{Y}_{i,j} = \{y_{k,l} ; k < i, 0 \le l \le M \text{ et } y_{i,l} ; l < j\} \\[1em] \bar{Y}_{i,0} = \{y_{k,l} ; k < i, 0 \le l \le M\} \\[1em] \bar{Y}_{0,j} = \{y_{0,l} ; l < j\} \\[1em] \bar{Y}_{0,0} = \text{ensemble vide.} \end{array} \right. $$

On peut alors décomposer, à l'aide de la formule de BAYES, $P(Y_{N,M}/\theta)$ en un produit :

$$ (29) \qquad P(Y_{N,M}/\theta) = \prod_{\substack{i=0, M \\ j=0, N}} P(y_{i,j}/\bar{Y}_{i,j}, \theta). $$

Soit :

$$ (30) \qquad \log P(Y_{N,M}/\theta) = \sum_{\substack{i=1, N \\ j=1, M}} \log P(y_{i,j}/\bar{Y}_{i,j}, \theta). $$

En posant :

$$ (31) \qquad \left\{ \begin{array}{l} v_{i,j} = y_{i,j} - E[y_{i,j}/\bar{Y}_{i,j}, \theta] \\[1em] R_{i,j} = E[v_{i,j}^2] \end{array} \right. $$

l'équation (30) devient (puisqu'il s'agit de lois gaussiennes) :

$$ (32) \qquad \log P(Y_{N,M}/\theta) = -\frac{1}{2} \left( \sum_{i,j} v_{i,j}^2 / R_{i,j} + \log R_{i,j} \right) $$
lorsque N et M sont assez grands l'équation (32) est remplacée par :

$$ (33) \qquad \log P(Y_{N,M}/\theta) = -\frac{1}{2} \left( \sum_{i,j} \frac{v_{i,j}^2}{R_{\infty}} + \log R_{\infty} \right) $$

où $v_{i,j}$ et $R_{\infty}$ sont donnés respectivement par les équations (20) et (11).

Les résultats ci-dessus permettent de traiter le signal d'image comme un signal à un indice, tout en tenant compte de la corrélation bidimensionnelle de l'image qu'il représente.

Dans la suite ces résultats sont utilisés d'une part pour étendre au cas des images des techniques classiques de détection d'hypothèses et de compression de données, d'autre part pour aborder des problèmes d'estimation par maximum de vraisemblance qui fournissent de nouveaux algorithmes d'identification et d'optimisation statistique.
### III. APPLICATIONS PRATIQUES

Complétant les problèmes de filtrage statistique en présence de bruit additif qui représentent la motivation initiale de ces travaux, on présente ci-dessous deux domaines d'applications nouveaux.

#### III. 1. Détection de signaux.

**a) Rappel sur les tests d'hypothèses [VAN TREES]** :

Les problèmes de détection se posent de manière classique sous la forme d'un choix entre diverses hypothèses sur la constitution du signal observé.

Dans le cas du choix entre hypothèses simples (binaire : $H_0, H_1$, resp. multiples : $H_0, \dots, H_{m-1}$) la décision repose sur la valeur prise par le (resp. les) rapport de vraisemblance :

$$ (34) \qquad \left\{ \begin{array}{l} \Lambda_{i,j} = \frac{P\{Y_{i,j} / H_0\}}{P\{Y_{i,j} / H_1\}} \\[1em] \text{(resp. } \Lambda_{i,j}(k,l) = \frac{P\{Y_{i,j} / H_k\}}{P\{Y_{i,j} / H_l\}} \text{)} \\[1em] \text{où } Y_{i,j} \text{ a été défini au chapitre II.} \end{array} \right. $$

Les valeurs de ces rapports de vraisemblance déterminent en effet les zones de décision entre les diverses hypothèses.

Dans le cas d'hypothèses binaires par exemple, une borne $\eta$ sépare les zones de décision entre les deux hypothèses, la valeur de la borne $\eta$ est fixée en tenant compte :

- soit de probabilités à priori d'occurence des hypothèses et de coûts associés aux décisions prises (méthode de BAYES) ;

- soit de probabilités d'erreur qu'il convient de ne pas dépasser ou de minimiser (méthode de NEYMAN-PEARSON).
Dans le cas du choix entre hypothèses composites (par exemple signal à détecter dépendant de paramètres θ) la décision dépend également de valeurs de rapport de vraisemblance généralisés, par exemple :

$$ (35) \qquad \Lambda_{i,j} = \frac{\int_{\theta_0} P(Y_{i,j}/\theta_0, H_0) dp(\theta_0)}{\int_{\theta_1} P(Y_{i,j}/\theta_1, H_1) dp(\theta_1)} $$

qui donne la solution optimale au sens de BAYES dans le cas où les lois de probabilités à priori de $\theta_0$ et $\theta_1$ sont connues.

Une solution sous-optimale pour le cas où les paramètres $\theta_0$ et $\theta_1$ sont déterministes mais inconnus est basée sur le rapport :

$$ (36) \qquad \Lambda_{i,j} = \frac{\max_{\theta_0} P(Y_{i,j} | \theta_0, H_0)}{\max_{\theta_1} P(Y_{i,j} | \theta_1, H_1)} $$

On voit ainsi que la plupart des problèmes de détection se résolvent à l'aide de calculs de rapports de vraisemblance.

Les deux exemples élémentaires qui suivent sont destinés à montrer que les résultats présentés au chapitre II permettent dans le cas de détection de signaux d'image, de réduire la difficulté liée à la corrélation bidimensionnelle de l'image.

b) Détection d'un signal déterministe connu dans un bruit coloré.

Soit à choisir, à partir de données $Y_{i,j}$, entre les deux hypothèses :

$$ H_0 : \quad y_{i,j} = s_{i,j} + z_{i,j} $$

$$ H_1 : \quad y_{i,j} = z_{i,j} $$
où {$s_{i,j}$} est un signal déterministe connu, et {$z_{i,j}$} est une suite gaussienne de la classe $\Sigma$ ($R = 0$ ou $R > 0$).

Il suffit alors de calculer, avec les mêmes paramètres, les innovations {$v_{i,j}^0$} et {$v_{i,j}^1$} de la suite {$z_{i,j}$} sous chacune des hypothèses.

Le logarithme du rapport de vraisemblance s'écrit alors, compte tenu de l'équation (33) :

$$ (37) \qquad \log \Lambda_{i,j} = + \frac{1}{2} \sum_{k,l} \frac{(v_{k,l}^1)^2 - (v_{k,l}^0)^2}{R_{\infty}}. $$

c) Détection d'un signal aléatoire de la classe $\Sigma$, noyé dans un bruit blanc de variance $\sigma_b^2$.

Les hypothèses sont dans ce cas :

$$ H_0 : \quad y_{i,j} = z_{i,j} + w_{i,j} $$

$$ H_1 : \quad y_{i,j} = w_{i,j} $$

Sous l'hypothèse $H_0$, {$y_{i,j}$} est une suite de la classe $\Sigma$, et on a :

$$ (38) \qquad \log P(Y_{i,j}/H_0) = -\frac{1}{2} \left( \frac{v_{k,l}^2}{R_{\infty}} + \log R_{\infty} \right). $$

Sous l'hypothèse $H_1$, on a trivialment :

$$ (39) \qquad \log P(Y_{i,j}/H_1) = -\frac{1}{2} \sum_{k,l} \left( \frac{y_{k,l}^2}{\sigma_b^2} + 2 \log \sigma_b \right). $$

Les exemples ci-dessus ne sont que démonstratifs. On peut résoudre
aussi aisément les problèmes suivants :

i) Choix entre deux signaux déterministes connus noyés dans un bruit coloré.

ii) Choix entre deux signaux aléatoires noyés dans un bruit blanc.

iii) Hypothèses multiples, test séquentiel... [PICINBONO, WALD].

### III. 2. Codage optimal et application à la compression.

Le calcul de l'innovation permet de transformer l'image $\{y_{i,j}\}$ en une image $\{v_{i,j}\}$ équivalente à la première et formée d'échantillons indépendants d'une loi gaussienne de variance minimum $R_\infty$.

Le codage optimal (minimum d'éléments binaires) consiste à générer un code (nécessairement de longueur variable) équivalent à l'image et minimisant la redondance d'information. En négligeant l'effet de la discrétisation initiale, le codage optimal se ramène donc à celui d'une loi gaussienne de variance $R_\infty$.

L'algorithme optimal (†) de compression sans distorsion est donc :

<figure><img src="image_4.png" /></figure>

<footnote>(†) Cet algorithme est optimal au sens statistique sur un ensemble de réalisations d'une suite de la classe Σ.</footnote>
Des schémas pratiques de compression basés sur la non transmission de l'innovation quand celle-ci est inférieure à un seuil sont envisageables.

Dans ce cas on coderait les valeurs de l'innovation supérieures au seuil par un code à longueur fixe suivi de la place de l'échantillon transmis par rapport à l'échantillon précédent.
Soit un échantillon $Y_{N,M}$ d'une suite gaussienne à deux indices, de
la classe Σ (éventuellement bruitée) et soit $\theta$ l'ensemble des paramètres du
"filtre blanchissant" associé à cette suite, c'est-à-dire le système schématisé
par :

<figure><img src="image_5.png" /></figure>

Le problème est de déterminer le paramètre θ qui maximise la proba-
bilité d'occurence de l'échantillon $Y_{N,M}$ observé, soit :

$$
(40) \quad \max_{\theta} P(Y_{N,M}/\theta).
$$

Or, d'après (33) :

$$
(41) \quad \text{Log } P(Y_{N,M}/\theta) = -\frac{1}{2} \sum_{\substack{i=1,N \\ j=1,M}} \left( \frac{v_{i,j}^2}{R_{\infty}} + \text{Log } R_{\infty} \right).
$$

En posant :

$$
(42) \quad J(\theta) = - \text{Log} P(Y_{N,M}/\theta)
$$

le problème posé par (40) est donc :

$$
(43) \quad \min_{\theta} J(\theta) = \min_{\theta} \sum_{\substack{i=1,N \\ j=1,M}} \left( \frac{v_{i,j}^2}{R_{\infty}} + \text{Log } R_{\infty} \right).
$$

Il est nécessaire pour le résoudre de définir les paramètres indé-
pendants du "filtre blanchissant" en précisant les contraintes qui pèsent sur
ces paramètres, contraintes liées d'une part à la positivité, d'autre part à l'existence d'un modèle de la classe Σ sous-jacent.

IV. 1. Cas général (n > 1 ou R > 0).

A. Détermination des paramètres "indépendants".

Formellement, l'innovation $v_{i,j}$ s'obtient à partir des observations
$\{y_{i,j}\}$ par les équations :

$$
\begin{equation} \tag{44}
\left\{
\begin{aligned}
v_{i,j} &= y_{i,j} - H \hat{x}_{i,j} - H H_{\infty} \eta_{i,j} \\
\hat{x}_{i,j} &= F_1 (\hat{x}_{i-1,j} + \hat{x}_{i-1,j}) \\
\hat{x}_{i,j} &= H_{\infty} \eta_{i,j} + L_{\infty} \gamma_{i,j} . \\
\eta_{i,j} &= F_{\infty} \eta_{i,j-1} + T_{\infty} v_{i,j-1} \\
\gamma_{i,j} &= (F_{\infty} - T_{\infty} H H_{\infty})' \gamma_{i,j+1} + H_{\infty}' H_{\infty}' R_{\infty}^{-1} v_{i,j} .
\end{aligned}
\right.
\end{equation}
$$

On voit apparaître ainsi 7 paramètres H, F₁, H∞, F∞, T∞, L∞, R∞. Les équations (9) à (18) permettent déjà de remarquer que T∞ s'exprime en fonction des autres paramètres par :

$$
(45) \qquad T_{\infty} = F_{\infty} L'_{\infty} H' R_{\infty}^{-1}.
$$

On remarquera alors que les 6 paramètres restants H, F₁, H∞, F∞, L∞, R∞
sont algébriquement équivalents aux paramètres H, F₁, H∞, F∞, G∞, R, grâce aux
relations :

$$
\[
(46) \quad \left\{
\begin{aligned}
R &= - H H_{\infty} L_{\infty} H' + R_{\infty} \\
G_{\infty} &= (H_{\infty} P_{\infty} + L_{\infty})' \\
P_{\infty} - F_{\infty} P_{\infty} F_{\infty}' &= T_{\infty} R_{\infty} T_{\infty}'
\end{aligned}
\right.
\]
$$
$P_\infty$ étant également donné par la limite de l'équation (9).

Ces paramètres caractérisent donc le "filtre innovateur".

Il faut maintenant exprimer les contraintes liées à l'existence
d'un modèle $(H, F_1, F_2, Q, R)$.

B. Expression des contraintes sur les paramètres $H_\infty, F_\infty, G_\infty, R_\infty, H, F_1$.

a) **Contraintes de stabilité :**

Les matrices $F_1, F_\infty, F_2$ doivent avoir leurs valeurs propres de module inférieur à l'unité.

b) **Contraintes de positivité :**

- positivité du spectre $(H(S_\infty(z)H' + R))$ :

$$
(47) \qquad R_{\infty} \ge 0 \quad (\text{covariance de l'innovation}).
$$

- positivité de la covariance du bruit d'observation :

$$
(48) \quad R = R_{\infty} - H H_{\infty} L_{\infty} H' \ge 0.
$$

- positivité du spectre $S_\infty(z)$ :

$$
L_{\infty} = H_{\infty} \Sigma_{\infty}
$$

$$
(49) \qquad \Sigma_{\infty} - F_{\infty} \Sigma_{\infty} F_{\infty}' = Q_{\infty} \ge 0.
$$

c) Contraintes directement liées à l'existence de $F_2$ et $Q$.

La relation entre $S_\infty(z)$ et $Q(z)$ donnée par la limite de l'équation (9) ne pouvant être qu'approchée, on ne peut introduire cette contrainte qu'en augmentant $\theta$ de 2 paramètres nouveaux ($Q$ et $F_2$) et en considérant une nouvelle fonction de coût de la forme :

$$
(5C) \quad J_1(\theta) = J(\theta) + \|\hat{Q}(z) - (zI - F_2)^{-1}Q(z^{-1}I - F_2')^{-1}\|_2
$$
où :

$$ (51) \qquad \hat{Q}(z) = S_{\infty}(z) - F_1(S_{\infty}(z) - S_{\infty}(z)H'(HS_{\infty}(z)H' + R)^{-1}HS_{\infty}(z))F_1' $$

|| || est une norme quelconque pour les fonctions définies sur le cercle unité $|z| = 1$, avec les contraintes supplémentaires :

$$ (52) \qquad Q \geq 0 $$

$$ (53) \qquad \hat{Q}(z) \geq 0 \quad \text{pour } |z| = 1. $$

On peut noter toutefois qu'une condition suffisante pour l'équation (53) est :

$$ (54) \qquad S_{\infty}(z) - F_1 S_{\infty}(z) F_1' \geq 0. $$

Cette dernière condition peut se convertir en une condition algébrique de positivité en fonction des matrices $F_1, H_\infty, F_\infty, G_\infty$.

Remarque :

Les contraintes (47) à (49) et (52) à (54) sont définies pour le cas général où les dimensions p, n et $n_\infty$ (dimension de $F_\infty$) sont quelconques. Dans le cas très souvent pratiqué où $p = n = 1$ la condition (53), la plus gênante, est redondante. Si de plus $n_\infty = 1$ (équations récurrentes scalaires) l'ensemble des contraintes se résume à :

$$ (55) \qquad \left\{ \begin{array}{l} |F_1| < 1, |F_2| < 1, |F_\infty| < 1 \\ R_\infty \geq L_\infty \geq 0 \\ Q \geq 0 \end{array} \right. $$

et on peut prendre $H = H_\infty = 1$.
Les résultats précédents sont résumés dans le théorème suivant :

### THEOREME.

Le problème d'estimation par maximum de vraisemblance des paramètres $\theta$ du "filtre blanchissant" revient à minimiser la fonction de coût (50) par rapport aux paramètres $H,F_1,F_2,Q,R,H_\infty,F_\infty,L_\infty,R_\infty$, compte tenu des équations (44) et (45) et des contraintes exprimées dans le paragraphe B.

Dans le cas où tous les paramètres ci-dessus sont scalaires l'ensemble de ces contraintes se réduit aux inégalités (55). ■

### Remarques :

Le problème d'optimisation posé ci-dessus paraît très complexe dans le cas général ($n > 1$). Remarquons toutefois que les méthodes de programmation non linéaire ont déjà été utilisées avec succès [notamment par MEHRA] pour l'estimation des paramètres de signaux à un indice, problème de complexité comparable à celui posé ci-dessus lorsque tous les paramètres sont scalaires.

- Les conditions de positivité et les propriétés de factorisation spectrale associées ont été largement utilisées dans ce paragraphe. Ces notions sont étudiées et détaillées, par exemple dans [FAURRE, CLERGET, GERMAIN], [ANDERSON].

## IV. 2. Cas particulier, $n = 1, R = 0$.

Ce cas particulier qui donne une solution très simple ne doit pas être négligé car il existe de nombreuses situations où l'on a besoin d'extraire simplement et rapidement l'innovation d'une image (ou d'une partie d'image). Il semble raisonnable dans de tels cas d'approcher l'image par un modèle de dimension $n = 1$.

Compte tenu de l'équation (26) le problème de maximisation devient :

$$ (56) \qquad \min_{\Theta} J(\Theta) = \sum_{i=1,N} (Q^{-1} v_{i,j}^2 + \text{Log } Q) \\ \qquad \qquad \sum_{j=1,M} $$
avec la contrainte :

$$ (57) \qquad v_{i,j} = y_{i,j} - F_1 y_{i-1,j} - F_2 y_{i,j-1} + F_1 F_2 y_{i-1,j-1}. $$

L'ensemble des paramètres $\theta$ est évidemment $(F_1, F_2, Q)$.

Les valeurs $(\hat{F}_1, \hat{F}_2, \hat{Q})$ qui annulent le gradient $\frac{\partial J}{\partial \theta}$ sont données par :

$$ (58) \qquad \sum_{\substack{i=1,N \\ j=1,M}} (y_{i,j} - \hat{F}_1 y_{i-1,j} - \hat{F}_2 y_{i,j-1} + \hat{F}_1 \hat{F}_2 y_{i-1,j-1}) (\hat{F}_1 y_{i-1,j-1} - y_{i,j-1}) = 0 $$

$$ (59) \qquad \sum_{\substack{i=1,N \\ j=1,M}} (y_{i,j} - \hat{F}_1 y_{i-1,j} - \hat{F}_2 y_{i,j-1} + \hat{F}_1 \hat{F}_2 y_{i-1,j-1}) (\hat{F}_2 y_{i-1,j-1} - y_{i,j-1}) = 0 $$

$$ (60) \qquad \hat{Q} = \sum_{\substack{i=1,N \\ j=1,M}} (y_{i,j} - \hat{F}_1 y_{i-1,j} - \hat{F}_2 y_{i,j-1} + \hat{F}_1 \hat{F}_2 y_{i-1,j-1})^2. $$

Les équations (58) et (59) se transforment pour donner :

$$ (61) \qquad \hat{F}_1 = \frac{\sum_{\substack{i=1,N \\ j=1,M}} (y_{i,j} - \hat{F}_1 y_{i-1,j} - \hat{F}_2 y_{i,j-1}) (y_{i-1,j} - \hat{F}_1 y_{i-1,j-1} - \hat{F}_2 y_{i,j-1})}{\sum_{\substack{i=1,N \\ j=1,M}} (y_{i,j} - \hat{F}_1 y_{i-1,j} - \hat{F}_2 y_{i,j-1})^2} \approx \frac{\Lambda_{z_1} (1,0)}{\Lambda_{z_1} (0,0)} $$

où

$$ (62) \qquad \left\{ \begin{array}{l} \Lambda_{z_1} (k,1) = E[z_1, i+k, j+1]^2 z_1, i, j \\ \\ z_1, i, j = y_{i,j} - \hat{F}_2 y_{i,j-1} \end{array} \right. $$

et

$$ (63) \qquad \hat{F}_2 = \frac{\sum_{\substack{i=1,N \\ j=1,M}} (y_{i,j} - \hat{F}_1 y_{i-1,j} - \hat{F}_2 y_{i,j-1}) (y_{i,j-1} - \hat{F}_1 y_{i-1,j-1} - \hat{F}_2 y_{i,j-1})}{\sum_{\substack{i=1,N \\ j=1,M}} (y_{i,j} - \hat{F}_1 y_{i-1,j} - \hat{F}_2 y_{i,j-1})^2} \approx \frac{\Lambda_{z_2} (0,1)}{\Lambda_{z_2} (0,0)} $$
$$
\text{où}
\qquad (64) \quad \left\{
\begin{aligned}
\Lambda_{z_2}^{(k,l)} &= E[z_{2,i+k,j+1} z_{2,i,j}] \\
z_{2,i,j} &= y_{i,j}^{-\hat{F}_1} y_{i-1,j}^{\hat{F}_2}
\end{aligned}
\right.
$$

Ces équations montrent déjà que pour N et M assez grands, les valeurs de $\hat{F}_1$ et $\hat{F}_2$ qui annulent les dérivées $\frac{\partial J}{\partial \hat{F}_1}$ et $\frac{\partial J}{\partial \hat{F}_2}$ sont nécessairement de module inférieur à 1.

On n'a donc à résoudre que le système d'équations (61) et (63), (par
substitution de la valeur de $\hat{F}_1$ dans l'expression de $\hat{F}_2$). Le choix entre les
cinq couples $(\hat{F}_1, \hat{F}_2)$ obtenus peut être fait en évaluant pour chacun d'eux la
valeur de $J(\theta)$. Toutefois comme :

$$
(65) \qquad \begin{array}{l} J(\theta) = NM(1 + \text{Log } \hat{Q}). \\ \min \end{array}
$$

Il suffira de comparer les valeurs de $\hat{Q}$ déduites de $(\hat{F}_1, \hat{F}_2)$ par (60)
et de choisir $\hat{F}_1$ et $\hat{F}_2$ qui donnent la valeur minimum de $\hat{Q}$.
L'objectif visé par les développements précédents était de démontrer que les modèles statistiques d'images introduits dans des publications précédentes permettent effectivement de résoudre théoriquement à l'aide d'algorithmes récursifs un grand nombre de problèmes classiques de traitement du signal.

Compte tenu du caractère particulier des images, (non stationnarités...) il apparaît indiqué en pratique d'utiliser des modèles approchés mais simples (modèle d'ordre 1), le caractère approché devant être largement compensé par la possibilité de réaliser des traitements locaux (sur une partition de l'image) avec des paramètres adaptés aux statistiques locales. L'utilisation des techniques d'estimation par maximum de vraisemblance devrait être avantageuse dans ce domaine.
REFERENCES

ANDERSON B.D.O.
HO P.J.
"Spectral Factorization of a Finite-Dimensional Non-stationary Matrix Covariance". IEEE AC19 N° 6. pp.680-692. (December 1974).

ANDERSON B.D.O.
"Algebraic Properties of Minimal Degree Spectral Factors Automatica". Vol.9. pp.491-500. (1973).

ATTASI S.¹
"Introduction d'une classe de systèmes linéaires à deux indices". Comptes rendus, 277, Série A. pp.1135-1138. (1973).

ATTASI S.²
"Introduction d'une classe de suites gaussiennes à deux indices". Comptes rendus, 278, Série A. pp.1375-1378. (1974).

ATTASI S.³
"Filtrage et lissage récursifs de suites gaussiennes à deux indices". Comptes rendus, 279, Série A. pp.903-906. (1974).

FAURRE P.
"Réalisations Markoviennes de Processus Stationnaires". Rapport LABORIA N° 13. (Mars 1973).

FAURRE P.
"Systèmes linéaires positifs". DUNOD, à paraître.

CLERGET M.

GERMAIN F.

KAILATH T.
"An Innovation Approach to Least-Squares Estimation". IEEE Trans AC13 N° 6. pp.646-660. (Décembre 1968). IEEE Trans AC18 N° 6. pp.588-600. (Décembre 1973).

KALMAN R.C.
FALB
"Topics in Mathematical System Theory". Mac Graw Hill. (1969).

ARBIB
MEHRA R.K. "Computational Aspects of Maximum Likelihood Estimation and Reduction in Sensitivity Function Calculations". IEEE Trans AC19 N° 6. (December 1974).

PICINBONO B.
VEZZOSI G. "Extension du critère de Neyman-Pearson en détection d'hy- pothèses multiples". 5ème colloque national sur le traitement du signal et ses applications. (Juin 1975).

VAN TREES H.L. "Detection, Estimation, and Modulation Theory". Part I and II. Wiley. (1971).

WALD A. "Sequential Analysis". Wiley. (1949).