Applying doubly stochastic filters to evaluate the dynamics of
object sizes on satellite image sequences

                V E Dementyev1, D S Kondratyev1


                1
                 Ulyanovsk State Technical University, ul. Severny Venets, 32, Ulyanovsk, Russia, 432027


                e-mail: dve@ulntc.ru

                Abstract. One of the important tasks facing the regional authorities is to monitor the condition
                of roads and power lines. In the Ulyanovsk region more than 8 thousand km of power lines and
                more than 9 thousand km of roads (including rural). A significant part of these facilities is
                located outside the settlements in places with medium and low availability. In many such
                places there is a problem of uncontrolled forest overgrowth. This work is devoted to solving
                the problem of automated satellite monitoring of such areas. For this purpose, it is proposed to
                use a modified convolutional neural network that processes time sequences of multispectral
                satellite images and allows to allocate territories occupied by forest and undergrowth with high
                accuracy. This approach allows us to assess the dynamics of overgrowth of the territory and
                perform the appropriate forecast with sufficient accuracy for practice.


1. Introduction
One of the important tasks of the satellite image processing is its thematic mapping, i.e. division of
image into identifiable areas containing the similar visual, correlation or texture characteristics of
pixels. The use of standard segmentation algorithms [1-3] for thematic mapping of satellite images
usually leads to significant errors caused by two reasons. Firstly, these algorithms are largely
incapable of taking into account the multi-zonal nature of remote sensing (RS), so each satellite image
contains the results of the Earth's surface registration in different spectral ranges. Some works [4-6]
suggest the possibility of processing hyperspectral images. Thus, the authors based their theory on
criterion of uniformity for reception of connected areas of such hyperspectral image [4], modification
and generalization of algorithm K-means [5] and use of physical properties of a satellite data [6].
Secondly, the existing approaches unable to use data on the observed territory received at previous
points of time for image segmentation. But using such data can significantly improve the quality of
processing at the expense of a fundamentally larger amount of information, but it is fraught with
difficulties of aggregation.
    It is possible to overcome the mentioned disadvantages by using neural network procedures of
segmentation and classification of multidimensional data. In work [9] the variant of the U-NET neural
network with full-connected layers (FCN) modification is presented. Herewith the input layer of the
network consisting from spectral layers of the multizonal image has been extended by three auxiliary
2d halftone images obtained from the original using NDVI, EVI, SAVI transformations and two 2d
arrays , representing the segmentation results at the previous point of time and one year ago. The use
of two such reference markings allows one's to reduce the error of classification in case of rapid
changes of the terrain due to the change of year time and in case of marking array absence at the


                    V International Conference on "Information Technology and Nanotechnology" (ITNT-2019)
Image Processing and Earth Remote Sensing
V E Dementyev, D S Kondratyev



previous point of time in connection, for example, with cloudiness. The quality analysis of such an
algorithm shows the processing comparable accuracy to the qualified operator results.
    It should be noted that the results of the satellite images thematic processing provides information
related to the peculiarities of the Earth's surface in the previous and current moments of time.
However, the necessary element of the RS processing system is the tool that allow to forecast the state
of certain objects and form recommendations for responsible persons. For example, let's consider the
task related to the forecasting about the dangerous convergence possibility of various natural and
technogenic objects. An example of such convergence may be the gradual growth of the forest area
along roads or power lines or landslide processes leading to the destruction of various infrastructures.

2. Satellite Image Processing
Let us formulate the given problem as follows. Let there be an aggregate of segments describing a
certain extended object 𝐓𝐨 (a road, an electric power transmission line etc.). The given aggregate is
usually described by a vector object having geographical reference in absolute coordinates. Let us
separate on this extended object a set of points 𝐓𝐨𝐢 = (𝐱𝟎𝐢 , 𝐲𝐨𝐢 ) having the distance between them
equal to ∆𝐨 . Let us assume that at a certain time instant 𝐭 next to the object there lies a certain extended
domain 𝐆𝐑 𝐭 defined by a set of points each corresponding to a pixel in the original halftone image. For
each point 𝐓𝐨𝐢 let us construct a perpendicular to the segment [𝐓𝐨𝐢 , 𝐓𝐨𝐢+𝟏 ]. Let us find the point
𝐓𝐄𝐢 𝐭 = �𝐳𝐄𝐱𝐢
          𝐭       𝐭
              , 𝐳𝐄𝐲𝐢 � of intersection of this perpendicular and the domain 𝐆𝐑 𝐭. Obviously the set of
points {𝐓𝐄𝐢 𝐭 , 𝐢 = 𝟏, . . , 𝐍𝐨 } describe a conditional boundary of the domain 𝐆𝐑 𝐭 as viewed from the
extended object. It enables to use estimates of the points coordinates {𝐓𝐄𝐢 𝐭 , 𝐢 = 𝟏, . . , 𝐍𝐨 } obtained at
different times 𝐭 = 𝟏, . . , 𝐓 as a source of information on the domain 𝐆𝐑 𝐭 dynamics for the purpose of
forming a prediction about its boundaries at future times 𝐭 > 𝑻 as well. (Figure 1).




   Figure 1. Geometric interpretation of the distance estimate between the stationary object and the
                                            movable area.

   In view of errors arising at satellite images recording and processing we get the following relations
(Eq. 1), in which 𝒏𝒕𝑬𝒙𝒊 and 𝒏𝒕𝑬𝒚𝒊 – white noise samples with zero mean and variance 𝝈𝟐𝒏 .
                            𝑡        𝑡      𝑡     𝑡      𝑡      𝑡
                           𝑧𝐸𝑥𝑖 = 𝑥𝐸𝑖   + 𝑛𝐸𝑥𝑖 , 𝑧𝐸𝑦𝑖 = 𝑦𝐸𝑖 + 𝑛𝐸𝑦𝑖 , 𝑖 = 1, . . , 𝑁𝑜 , 𝑡 = 1, . . , 𝑇, (1)
         𝒕        𝒕                                                                𝟐
where 𝒏𝑬𝒙𝒊 and 𝒏𝑬𝒚𝒊 – white noise samples with zero mean and variance 𝝈𝒏 .
   Direct measurements based on the results of satellite material image-type related mapping show
that MSE �𝛔𝟐𝐧 is approximately equal to 𝟏. 𝟓𝐃𝐱𝐲 , where 𝐃𝐱𝐲 – resolution of the original images.
   Let us suppose that the boundary of the domain 𝑮𝑹 𝒕 can move non-uniformly. Thus, for example,
the area of a precipice or a ravine can increase by tens of centimeters per each year and at a certain
moment its rate might increase exponentially. Then we will use doubly stochastic (DS) model [3,10]
to describe unknown coordinates in the following form (Eq. 2). In this case 𝐫𝐚𝐱 , 𝐫𝐚𝐲 – scalar parameters
determining change potential for the accelerations 𝐚𝐭𝐄𝐱𝐢 and 𝐚𝐭𝐄𝐲𝐢 ; 𝛏𝐭𝐚𝐱𝐢 , 𝛏𝐭𝐚𝐲𝐢 – independent normal
random variables with zero mean and variance 𝛔𝟐𝛏 :
                                𝑡         𝑡−1     𝑡−2    𝑡     𝑡−1      𝑡−2
                               𝑥𝐸𝑖 = 2𝑥𝐸𝑖     − 𝑥𝐸𝑖   + 𝑎𝐸𝑥𝑖 �𝑥𝐸𝑖  − 𝑥𝐸𝑖    �,


V International Conference on "Information Technology and Nanotechnology" (ITNT-2019)                      55
Image Processing and Earth Remote Sensing
V E Dementyev, D S Kondratyev



                                𝑡        𝑡−1       𝑡−2         𝑡     𝑡−1       𝑡−2
                              𝑦𝐸𝑖  = 2𝑦𝐸𝑖     − 𝑦𝐸𝑖      + 𝑎𝐸𝑦𝑖    �𝑦𝐸𝑖   − 𝑦𝐸𝑖    �,                            (2)
                             𝑡            𝑡−1       𝑡        𝑡             𝑡−1      𝑡
                          𝑎𝐸𝑥𝑖 = 𝑟𝑎𝑥 𝑎𝐸𝑥𝑖 + 𝜉𝑎𝑥𝑖 , 𝑎𝐸𝑦𝑖 = 𝑟𝑎𝑦 𝑎𝐸𝑦𝑖 + 𝜉𝑎𝑦𝑖               ,
  The model (2) can be rewritten in the form:
                        𝑋�𝐸𝑖
                           𝑡        𝑡
                               = 𝜑𝐸𝑥𝑖  (𝑋�𝐸𝑖
                                           𝑡−1
                                               ) + 𝜉𝑥𝑖̅𝑡 , 𝑌�𝐸𝑖
                                                             𝑡       𝑡
                                                                 = 𝜑𝐸𝑦𝑖  (𝑋�𝐸𝑖
                                                                            𝑡−1       ̅𝑡 ,
                                                                                ) + 𝜉𝑦𝑖
where
                                               𝑡−1               𝑡
                                             𝑥𝐸𝑖                𝑣𝐸𝑥𝑖              0
                         𝑡      � 𝑡−1
                       𝜑𝐸𝑥𝑖 �𝑋𝐸𝑖 � = � 0                 𝑡−1           𝑡
                                                        𝑣𝐸𝑥𝑖 (1 + 𝑎𝐸𝑥𝑖 )          0 �,
                                                                                𝑡−1
                                               0                 0             𝑎𝐸𝑥𝑖 𝑟𝑎𝑥
                                              𝑡−1                𝑡
                                             𝑦𝐸𝑖                𝑣𝐸𝑦𝑖              0
                                                         𝑡−1           𝑡
                         𝑡
                       𝜑𝐸𝑦𝑖   �𝑌�𝐸𝑖
                                  𝑡−1
                                      �=� 0            𝑣𝐸𝑦𝑖   (1 + 𝑎𝐸𝑦𝑖   )       0 �,
                                                                                    𝑡−1
                                                      0                0           𝑎𝐸𝑦𝑖 𝑟𝑎𝑦
                                               0               0
                                       ̅𝑡 = � 0 �, 𝜉𝑥𝑖
                                     𝜉𝑥𝑖             ̅𝑡 = � 0 �.
                                               𝑡               𝑡
                                              𝜉𝑎𝑥𝑖            𝜉𝑎𝑦𝑖
    The assigned notations enable to apply DS nonlinear filtering to observations and for construction
of predictions of behaviour of the region 𝑹𝟎 in reference to the object𝑻𝒐 . In doing so let us introduce
𝑿               � 𝒕−𝟏
� 𝒕Э𝑬𝒊 = 𝝋𝒕𝑬𝒙𝒊 (𝑿          �𝒕      𝒕    � 𝒕−𝟏
                  𝑬𝒊 ) and 𝒀𝑬𝒊 = 𝝋𝑬𝒚𝒊 (𝑿𝑬𝒊 )– extrapolated predictions for the point 𝐓𝐄𝐢 coordinates at
                                                       𝐭−𝟏
the time point 𝐭 based on the preceding observations 𝐳𝐄𝐱𝐢        𝐭−𝟏
                                                            and 𝐳𝐄𝐲𝐢  . Denote by 𝐏𝐱𝐢 𝐭−𝟏 , 𝐏𝐲𝐢 𝐭−𝟏 – filtering
                                                                           𝐓                          𝐓
error covariance matrices at the time point (𝐭 − 𝟏), 𝐕𝐱𝛏𝐢 𝐭 = 𝐌{𝛏�𝐭𝐱𝐢 𝛏�𝐭𝐱𝐢 }, 𝐕𝐲𝛏 𝐢 𝐭 = 𝐌{𝛏�𝐭𝐲𝐢 �𝛏𝐭𝐲𝐢 } – diagonal
covariance matrices for random increments𝝃�𝒕𝒙𝒊 . Then error covariance matrices for such extrapolation
have the following form:
                                                   𝑇
        𝑃 𝑡 = 𝑀 ��𝑋� 𝑡 − 𝑋� 𝑡 ��𝑋� 𝑡 − 𝑋� 𝑡 � � = 𝜑𝑡 ′(𝑋� 𝑡−1 )𝑃 𝑡−1 𝜑𝑡 ′(𝑋� 𝑡−1 )𝑇 + 𝑉 𝑡 ,
           Э𝑥𝑖              Э𝐸𝑖      Э𝐸𝑖      Э𝐸𝑖         Э𝐸𝑖          𝐸𝑥𝑖    𝐸𝑖    𝑥        𝐸𝑥𝑖   𝐸𝑖     𝑥𝜉 𝑖
                                                                𝑇
         𝑃Э𝑦𝑖 𝑡 = 𝑀 ��𝑌�Э𝐸𝑖
                        𝑡
                            − 𝑌�Э𝐸𝑖
                                𝑡
                                    ��𝑌�Э𝐸𝑖
                                        𝑡
                                            − 𝑌�Э𝐸𝑖
                                                𝑡          𝑡
                                                    � � = 𝜑𝐸𝑦𝑖 ′(𝑌�𝐸𝑖
                                                                   𝑡−1
                                                                       )𝑃𝑦𝑖 𝑡−1 𝜑𝐸𝑦𝑖
                                                                                 𝑡
                                                                                     ′(𝑌�𝐸𝑖 ) + 𝑉𝑦𝜉 𝑖 𝑡 ,
                                                                                         𝑡−1 𝑇


  If we denote by �    𝒙𝒕Э𝑬𝒊 ,𝒚
                              �𝒕Э𝑬𝒊 - the first elements of the vectors 𝑿   � 𝒕Э𝑬𝒊 and 𝒀
                                                                                       � 𝒕Э𝑬𝒊 ; 𝑩𝒙𝒊 𝒕 =
                −𝟏                         −𝟏
𝑷Э𝒙𝒊 𝒕 𝑪𝒙 𝑻 𝑫𝒕𝒙𝒊 ; 𝑩𝒚𝒊 𝒕 = 𝑷Э𝒚𝒊 𝒕 𝑪𝒚 𝑻 𝑫𝒕𝒚𝒊 ; 𝑫𝒕𝒙𝒊 = 𝑪𝒙 𝑷Э𝒙𝒊 𝒕 𝑪𝒙 𝑻 + 𝝈𝟐𝒏 ; 𝑫𝒕𝒚𝒊 = 𝑪𝒚 𝑷Э𝒚𝒊 𝒕 𝑪𝒚 𝑻 + 𝝈𝟐𝒏 ;
С𝒙 = С𝒚 = (𝟏 𝟎 𝟎), then we can write the following relations (Eq. 27) for DS coordinate filters:
                   𝑋� 𝑡 = 𝑋� 𝑡 + 𝐵 𝑡 �𝑧 𝑡 − 𝑥� 𝑡 �, 𝑌� 𝑡 = 𝑌� 𝑡 + 𝐵 𝑡 �𝑧 𝑡 − 𝑦� 𝑡 �.
                           𝐸𝑖       Э𝐸𝑖      𝑥𝑖     𝐸𝑥𝑖         Э𝐸𝑖   𝐸𝑖     Э𝐸𝑖        𝑦𝑖   𝐸𝑦𝑖   (3)
                                                                                                    Э𝐸𝑖
   The filtering error variance at each step is determined by the matrices 𝑷𝒙𝒊 𝒕 = (𝑬 − 𝑩𝒙𝒊 𝒕 )𝑷Э𝒙𝒊 𝒕 ,
𝑷𝒚𝒊 𝒕 = (𝑬 − 𝑩𝒚𝒊 𝒕 )𝑷Э𝒚𝒊 𝒕 .
   Special attention must be given to the fact that due to the necessity for the distance from the
extended object 𝑻𝒐 to the domain 𝑹𝟎 to be surveyed, it is possible to simplify the boundaries
                                                                                              𝐭
coordinates filtering process 𝐑 𝟎 by processing only one coordinate, namely, the distance 𝐱𝐫𝐢   from the
point 𝐓𝐨𝐢 belonging to 𝐓𝐨 up to the boundary 𝐑 𝟎 at the time point 𝐭. An aggregate of similar
                𝐭
observations 𝐳𝐄𝐱  𝐫𝐢
                     can be processed by a technique identical to the above-described one.
   As an illustration of such a technique in figures below series satellite images fragments for the
forest tract in Cherdakly district of the Ulyanovsk region for the period 2001-2017 years (figure 2) and
Milanovsky opencast colliery on riverbank of the Volga in the northern part of the city of Ulyanovsk
for the period 2013-2017 years (figure 3) are presented. Here for convenience of color image
perception and its recovery overlapping of visible spectral bands and superposition of the segmented
image fragment and normals to the object to be monitored is carried out. In the first case the number of
multispectral images to be processed amounted to 42 snapshots, in the second case - 32 snapshots. The
minimal time interval for satellite information production amounts to 14 days.




V International Conference on "Information Technology and Nanotechnology" (ITNT-2019)                             56
Image Processing and Earth Remote Sensing
V E Dementyev, D S Kondratyev




               July 2003                     June 2011                   July 2017
       Figure 2.Forest tract satellite images fragments and the results of these images processing.




         August 2015                       April 2016                         May 2016
   Figure 3. Milanovsky opencast colliery satellite images fragments and the results of these images
                                             processing.

    The above-mentioned groups of multispectral images were separated into a training and an
operating samples. The training samples were used to specify the filtering parameters, in particular, to
estimate the parameters 𝐫𝐚𝐱 and 𝛔𝟐𝛏 . The operating part of the sample was processed by three
                                                              𝐭+𝟏
algorithms enabling to carry out the distance prediction 𝐱𝐫𝐢      basing on the preceding observations.
                                                                      𝐭+𝟏       𝐭     𝐭
The first algorithm (I) involves constructing a simple prediction 𝐱�𝐫𝐢    = 𝟐𝐳𝐫𝐢  − 𝐳𝐫𝐢 wherein only the
variable 𝐱�𝐫𝐢 change speed is taken into account for the time interval (𝐭 − 𝟏, 𝐭). The second algorithm
                                                                𝐭
(II) assumes linear Kalman filtering of the observations 𝐳𝐫𝐢       and construction of the extrapolated
               𝐭+𝟏
             �
prediction 𝐱Э𝐫𝐢 based on the results of the processing. The third algorithm (III) is in the above-
                                                                              𝐭
described doubly stochastic filtering of an aggregate of the observations 𝐳𝐫𝐢     and construction of the
        � 𝐭+𝟏
vector 𝐗 Э𝐄𝐢 . In the table 1 below the values of prediction average errors depending on the object kind
are presented.
    On average the DS filter provides prediction accuracy 6% higher than in case of using Kalman
filter and 58% higher than in case of simple linear predictions. It enables to estimate coordinates and
rate change dynamics for boundaries of the domain to be monitored by using the DS filter. It is


V International Conference on "Information Technology and Nanotechnology" (ITNT-2019)                  57
Image Processing and Earth Remote Sensing
V E Dementyev, D S Kondratyev



essential that the DS filter enables to quicker respond to abrupt rate change of the processes
determining the object behaviour. As an illustration we provide the estimates behaviour for the
                                                                                                   � 𝒕𝒊
distance from opencast colliery to one of the points to be monitored (figure 4a) and the parameter 𝒂
estimate (figure 4b).

                                    Table 1. Prediction average errors.
                                                   Average error Average error                    Average error
                                                       for the          for the                  for the algorithm
                                                     algorithm I     algorithm II                        III
 The forest tract snapshot. October 2014
                                                        6.7 m            2.7 m                        2.6 m
 The forest tract snapshot. May 2015
                                                                    10.7 m              3.9 m         3.7 m
 The forest tract snapshot. June 2016
                                                                     6.2 m              3.6 m         3.3 m
 Milanovsky opencast colliery snapshot. May
 2014                                                                7.1 m              3.8 m         3.6 m

 Milanovsky opencast colliery snapshot. May
 2015                                                                7.3 m              3.9 m         3.8 m

 Milanovsky opencast colliery snapshot. May
 2016                                                                6.9 m              3.9 m         3.7 m

 Milanovsky opencast colliery snapshot. April
                                                                    12.4 m              8.9 m         7.8 m
 2016. Beginning of the avalanche processes.

 Milanovsky opencast colliery snapshot. April                       30.7 m              32.8 m       12.6 m
 2016. Continuation of the avalanche processes.

 Milanovsky opencast colliery snapshot. May                         20.3 m              18.1 m       17.3 m
 2016. Cessation of the avalanche processes.




                       Figure 4. Dependence of the filtering results versus survey time.

3. Conclusion
Direct analysis of the given results in comparison with the data of objective monitoring (solid line)
indicates superiority of the DS filter over conventional linear Kalman filter in filtering accuracy. As it
takes place, this superiority makes itself evident in the most distinct manner in case of abrupt change
of the rock collapsing process (and the corresponding reduction of the distance between the opencast

V International Conference on "Information Technology and Nanotechnology" (ITNT-2019)                           58
Image Processing and Earth Remote Sensing
V E Dementyev, D S Kondratyev



colliery and the point to be monitored). This change corresponds to a significant change of the
parameter 𝒂𝒕𝒊 estimate which enables to register considerable changes in the opencast colliery domain
state basing only on this estimate dependence nature versus survey time.

4. References
[1] Vasiliev K K 2016 Optimum processing of signals in discrete time (Moscow: Radiotehnika)
[2] Krasheninnikov V R and Vasil'ev K K 2018 Multidimensional image models and processing
     Intelligent Systems Reference Library 135 11-64
[3] Vasiliev K K and Dementyev V E 2017 Presentation and processing of satellite multi-zone
     images (Ulyanovsk)
[4] Fursov V A, Bibikov S A and Baida O A 2014 Thematic classification of hyperspectral images
     by conjugacy index Computer Optics 38(1) 154-158
[5] Zimichev E A, Kazansky N L and Serafimovich P G 2014 Spatial classification of hyperspectral
     images using the k-means ++ clustering method Computer Optics 38(2) 281-286
[6] Blokhinov Yu B, Gorbachev V A, Rakutin Yu O and Nikitin A D 2018 Development of an
     algorithm for semantic segmentation of real-time aerial photographs Computer Optics 42(1)
     141-148 DOI: 10.18287/2412-6179-2018-42-1-141-148
[7] Belim S V and Larionov S V 2016 Image segmentation algorithm based on searching
     communities on graphs Computer Optics 40(6) 904-910 DOI: 10.18287/2412-6179-2016-40-6-
     904-910
[8] Dementyev V E, Kondratyev D S and Frenkel A G 2018 Application of convolutional neural
     networks for thematic mapping of sequences of satellite multispectral images (Moscow:
     Radiotehnika)
[9] Dementyev V E 2017 The use of double-stochastic autoregressive models for processing
     satellite images (Moscow: Radiotehnika)
[10] Andriyanov N A and Dement'ev V E 2018 Application of mixed models of random fields for
     the segmentation of satellite images CEUR Workshop Proceedings 2210 219-226




V International Conference on "Information Technology and Nanotechnology" (ITNT-2019)              59