WO2024224134A1

WO2024224134A1 - Method, device and computer program for measuring a difference of fill ratio of a tank

Info

Publication number: WO2024224134A1
Application number: PCT/IB2023/000224
Authority: WO
Inventors: Roland AKIKI; Carlo DE FRANCHIS; Gabriele Facciolo; Raphaël GRANDIN; Jean-Michel Morel
Original assignee: Kayrros SAS
Current assignee: Kayrros SAS
Priority date: 2023-04-28
Filing date: 2023-04-28
Publication date: 2024-10-31
Anticipated expiration: 2025-10-28

Abstract

The invention relates to a method for measuring a difference of fill ratio of a tank over time, comprising the steps of downloading a time series of N+1 SAR images, aligning the SAR images, detecting in each SAR images at least one point (Q) selected among a first point (A) corresponding to an upper reflector (40) located on a circumference (20) of the tank (1) and a second point (B) corresponding to a bottom (21) of the tank, selecting K couples of SAR images, calculating for each couple an interferogram, estimating a topographic phase difference for each couple of SAR images, calculating a differential interferogram by subtracting the topographic phase difference from the interferogram for each couple, estimating the difference of fill ratio of the tank from the K differential interferograms taken at the at least one point (Q).

Description

Method, device and computer program for measuring a difference of fill ratio of a tank

BACKGROUND OF THE INVENTION

The invention concerns a method, a device and a computer program for measuring a difference of fill ratio of a tank.

One field of application of the invention is storage tanks.

A storage tank is a storage medium typically used for volatile liquids, such as crude oil. There are different types of storage tanks, among which we distinguish tanks that have a floating roof, and tanks that have a fixed roof. For floating roof tanks, the roof moves vertically as the volume of the oil changes to reduce evaporation loss.

Document US 10,672,139 B2 relates to a method for remotely measuring the volume of liquid stored in external floating roof tanks. However, the method known by this document does not enable to estimate a difference of fill ratio of fixed-roof tanks.

SUMMARY OF THE INVENTION

An aim of the invention is to obtain a method and device capable of estimating a difference of fill ratio of a fixed roof tank, in order to obtain information about the volume of liquid stored in the tank and/or a fill rate of the tank.

A first subject matter of the invention is a method for measuring a difference of fill ratio of a tank over time, comprising the following steps carried out by at least one calculator: downloading a time series of N+1 SAR images having pixels covering the tank and having been acquired by acquisitions made by a satellite from different positions relatively to the tank, wherein N is a first integer higher than 2 and is prescribed, aligning the SAR images by co-registration, detecting in each SAR images at least one point selected among a first point corresponding to an upper reflector located on a circumference of the tank and a second point corresponding to a bottom of the tank, selecting K couples of SAR images among the time series of SAR images, wherein K is a second integer and is prescribed, calculating for each couple of SAR images an interferogram, estimating a topographic phase difference of for each couple of SAR images from a digital elevation model, from satellite orbits during the acquisitions and from image timing information, calculating a differential interferogram by subtracting the topographic phase difference from the interferogram for each couple of SAR images, estimating the difference of fill ratio of the tank from the K differential interferograms taken at the at least one point.

The interferogram is related to a deformation of the tank, which is itself related to the difference of fill ratio of the tank between the two SAR images of each couple. Consequently, the differential interferogram is related to the difference of fill ratio of the tank between the two SAR images of each couple. The invention enables to estimate a difference of fill ratio of a fixed roof tank.

According to an embodiment of the invention, the upper reflector is an upper platform located on the circumference of the tank.

According to an embodiment of the invention, the at least one point is or comprises the first point corresponding to the upper reflector located on the circumference of the tank and the second point corresponding to a bottom of the tank, wherein estimating the difference of fill ratio of the tank is made from the K differential interferograms calculated for the first point corresponding to the upper reflector located on the circumference of the tank and based on the K differential interferograms calculated for the second point corresponding to the bottom of the tank.

According to an embodiment of the invention, the differential interferogram is calculated between the SAR image I, of the time series and the SAR image IM of the time series for i being a third integer going from 1 to N, wherein the SAR image l,_+i is consecutive to the SAR image I, in the time series, the K couples of SAR images are the N couples of SAR images l„ IM, wherein estimating the difference of fill ratio of the tank from the K differential interferograms is made based on the differential interferograms taken at the at least one point for each SAR image I, for i going from 1 to N.

According to an embodiment of the invention, the method comprises: choosing a search window of pixels surrounding a prescribed guard window containing the tank in each SAR image l„ calculating for each pixel Pc of the search window an amplitude dispersion DA(PC) equal to

wherein m_rA(P_c) is the mean of the amplitude values (c/Zj(P_c)) for the pixel Pc taken along the times series of N+1 SAR images,

<7^(Pc) ^1S the standard deviation of the amplitude values for the pixel Pc taken along the times series of N+1 SAR images, determining stable pixels as the pixels Pc, which are situated in the search window and outside the prescribed guard window and whose amplitude dispersion

is lower than or equal to a prescribed amplitude dispersion threshold, calculating a reference phase for the stable pixels for each SAR image l„ wherein estimating the difference of fill ratio of the tank is made based on a doublephase difference taken at the at least one point, which is equal to the differential interferograms taken at the at least one point to which the reference phase is subtracted for each SAR image I, for i going from 1 to N.

According to an embodiment of the invention, the method comprises: calculating the reference phase <p_ref,i for the stable pixels for each SAR image I, as being the differential interferogram calculated for the pixels Pc whose amplitude dispersion D_rA(P_c) is the lowest.

According to an embodiment of the invention, the method comprises: calculating the reference phase <p_ref,i for the stable pixels U for each SAR image h as being an angle

wherein z. is an angle

V designates the stable pixels U, a>t(U) is a weight for each stable pixel U, <p, is the differential interferogram calculated between the SAR image h of the time series and the SAR image l,_+i of the time series for i going from 1 to N.

According to an embodiment of the invention, (^(LT) = 1.

According to another embodiment of the invention, (^(LT) is calculated depending on the amplitude value c/Zj(U) of the stable pixel U in SAR image I, and on the amplitude value c/Zj₊₁(U) of the stable pixel U in SAR image IM .

According to an embodiment of the invention,

According to another embodiment of the invention,

According to another embodiment of the invention

According to an embodiment of the invention, the method comprises: calculating the reference phase <p_ref,i for the stable pixels for each SAR image I, as being an affine phase <p_re/,i(R) = ax_R + py_R + c for a point R of coordinates (x_K,y_K) in the SAR image l„ where a, ft and c are parameters fitted with where and

^ ^ wherein z. is an angle

Nb is a number of stable pixels in the search window and outside the prescribed guard window.

According to an embodiment of the invention, the double-phase difference 6<pi(Q) taken at the at least one point Q between the differential interferograms cpi(Q) taken at the at least one point Q and the reference phase <Pref,i for each SAR image I, for i going from 1 to N is equal to

where W(<p) is a wrap operator that transforms a phase <p into the ]-

n,n] interval. According to an embodiment of the invention, the method comprises: calculating from each couple of SAR images l„ IM for point Q an orthogonal baseline which is a projection of a line joining positions of at least one satellite having acquired the SAR images l„ IM in a direction, which is orthogonal to a line of sight of point Q defined for SAR image I, and which departs from the position for SAR image l„ wherein the projection of the line is made along the line of sight, estimating by applying to the double-phase difference 6<pi(Q) a linear regression model, which follows 8(pi(Q) = K(Q) • Bi(Q) + v_t(Q), which depends on the orthogonal baseline and in which K(Q) is a first variable parameter and Vi(Q) is a second variable parameter, and by varying K(Q) and v,(Q), the first variable parameter K(Q) being optimized into K(Q), which is equal to

where p is the power of the L^p -norm minimized, calculating an estimated residual double-phase difference 8<p_{i-topo res}{Q') equal to

calculating a compensated double-phase difference <5<Pi-comp(Q) by subtracting the estimated residual double-phase difference §<Pi-topo,res(.Q) from the double-phase difference <5cpi(Q). wherein estimating the difference of fill ratio of the tank is made based on the compensated double-phase difference <5<Pi-comp(Q) •

According to an embodiment of the invention, the method comprises: calculating the difference 8u_{i noisy}(Q) of fill ratio of the tank as a function of the compensated double-phase difference <5<Pi-comp(Q) and of a prescribed or calculated slope

According to an embodiment of the invention, the at least one point is the first point corresponding to the upper reflector located on the circumference of the tank and the second point corresponding to a bottom of the tank, the K differential interferograms are calculated for the first point corresponding to the upper reflector located on the circumference of the tank and based on the K differential interferograms calculated for the second point corresponding to the bottom of the tank, the method comprising estimating the difference Auj _idenoised °f fill ratio of the tank as

wherein C is a set of points, constituted by the first point corresponding to the upper reflector located on the circumference of the tank and by the second point corresponding to the bottom of the tank, qji(Q) is a reliability weight, higher than 0 and lower than 1 , which is associated to each point of the set of points,

and

is a prescribed or calculated slope ??(Q).

According to an embodiment of the invention,

According to another embodiment of the invention, ip; G?) = SRC(Q), wherein

wherein Z is a rectangular zone of pixels P_t of SAR image l„ which is situated around the point Q and from which the point Q, the pixels situated in a line to which the point Q belongs and the pixels situated in a column to which the point Q belongs are removed, c/Zj(<2) is an amplitude value for the point Q, c/Zj(Pt) is an amplitude value for the pixels P_t,

NP is the number of pixels in the zone Z.

According to an embodiment of the invention, the method comprises: estimating, by applying to the compensated double-phase difference

having been calculated for another tank another linear regression model, which follows

which depends on the compensated double-phase difference <5<Pi-_COmp(Q) and in which ??(Q) is a third variable parameter and f(Q) is a fourth variable parameter, and by varying ??(Q), the third variable parameter ??(Q) being optimized into slope which is equal to

where t is the power of the Z^-norm minimized, wherein Auj are other known differences of fill ratio of the other tank,

I is a fifth integer going from 1 to L, wherein L is a sixth prescribed integer higher than or equal to 1 , calculating the difference of fill ratio of the tank as a function of the compensated double-phase difference <5<Pi-_COmp(Q) and of the slope

According to an embodiment of the invention, the slope is estimated based on the method carried out on the other tank having a floating roof.

According to an embodiment of the invention, the slope estimated based on the method carried out on the other tank having a floating roof is used for estimating the difference Au_t .denoised of fill ^rati° of the tank having a fixed roof.

A second subject matter of the invention is a device for measuring a difference of fill ratio of a tank over time, comprising at least one calculator configured to carry out the following steps: downloading a time series of N+1 SAR images having pixels covering the tank and having been acquired by acquisitions made by a satellite from different positions relatively to the tank, wherein N is a first integer higher than 2 and is prescribed, aligning the SAR images by co-registration, detecting in each SAR images at least one point selected among a first point corresponding to an upper reflector located on a circumference of the tank and a second point corresponding to a bottom of the tank, selecting K couples of SAR images among the time series of SAR images, wherein K is a second integer and is prescribed, calculating for each couple of SAR images an interferogram, estimating a topographic phase difference of for each couple of SAR images from a digital elevation model, from satellite orbits during the acquisitions and from image timing information, calculating a differential interferogram by subtracting the topographic phase difference from the interferogram for each couple of SAR images, estimating the difference of fill ratio of the tank from the K differential interferograms taken at the at least one point.

A third subject matter of the invention is a non-transitory computer-readable storage medium having instructions that, when executed by a processor, cause the processor to execute steps for measuring a difference of fill ratio of a tank over time, wherein the steps comprise downloading a time series of N+1 SAR images having pixels covering the tank and having been acquired by acquisitions made by a satellite from different positions relatively to the tank, wherein N is a first integer higher than 2 and is prescribed, aligning the SAR images by co-registration, detecting in each SAR images at least one point selected among a first point corresponding to an upper reflector located on a circumference of the tank and a second point corresponding to a bottom of the tank, selecting K couples of SAR images among the time series of SAR images, wherein K is a second integer and is prescribed, calculating for each couple of SAR images an interferogram, estimating a topographic phase difference of for each couple of SAR images from a digital elevation model, from satellite orbits during the acquisitions and from image timing information, calculating a differential interferogram by subtracting the topographic phase difference from the interferogram for each couple of SAR images, estimating the difference of fill ratio of the tank from the K differential interferograms taken at the at least one point. DESCRIPTION OF THE DRAWINGS

The invention will be more clearly understood from the following description, given solely by way of non-limiting example in reference to the appended drawings, in which :

Figure 1 show is a schematic view of a fixed roof tank or floating roof in vertical cross-section, to which the invention may be applied,

Figure 2 is a diagram showing steps of a method for measuring a deformation of a tank according to embodiments of the invention.

Figure 3 is a schematic view of a deformation of a tank, which may be measured by the invention.

Figure 4 is a schematic view of a deformation of a tank, which may be measured by the invention.

Figure 5 is a schematic view of a device for measuring a difference of fill ratio of a tank according to the invention.

Figure 6 is a schematic view of windows in an image for calculating a reference phase according to embodiments of the invention.

Figure 7 is a diagram showing sub-steps of a method for measuring a difference of fill ratio of a tank according to embodiments of the invention. Figure 8 is a diagram showing sub-steps of a method for measuring the slope relating the compensated double phase difference at the roof or the base of the tank to the difference of the tank fill ratio according to embodiments of the invention.

Figure 9 shows schematically a selection of a zone of pixels which may be used in a method for measuring a difference of fill ratio of a tank according to embodiments of the invention.

DETAILED DESCRIPTION OF THE INVENTION

An example of a tank 1 is shown in Figure 1 .

Such a tank 1 comprises a circumference 20 formed by a cylindrical wall 2 (for example substantially vertical) defining its internal volume and a roof 3. In a first example, the roof 3 may be a floating roof 3 which floats freely on top of the liquid comprised within tank 1 , and therefore rises and falls as tank 1 is filled and emptied. Such tank is called floating roof tank 1 .

In another example, the roof 3 may be a fixed roof which does nor move and does not rise and does not fall as tank 1 is filled and emptied. Such tank is called fixed roof tank 1.

The invention may apply to a floating roof tank 1 or to a fixed roof tank 1 .

The liquid stored in the tank 1 may be crude oil or other.

On Figure 2, the method for measuring a difference of fill ratio of the tank 1 according to the invention has the steps described below. The invention provides also a device 10 for measuring a difference of fill ratio of the tank 1 .

In a first step S10, a time series of N+1 SAR images I, covering the tank 1 from different positions P, S of a satellite (one or more) relatively to the tank 1 are downloaded to a calculator CAL, as shown on figure 5. N is a first integer higher than 2 and is prescribed. Consequently, i is an integer going from 1 to N+1. The integers i correspond to different dates of acquisition of the SAR images I,.

Metadata comprising satellite orbit during the acquisition (time, position, velocity) and image timing information (first line time, first column time, line sampling rate, column sampling rate) are associated to each 1 SAR image and are downloaded to the calculator CAL CAL.

The device 10 may comprise the calculator CAL (one or more) to carry out the steps and operations mentioned of the method of the invention. The calculator CAL can comprise a computer (one or more), a processor (one or more), a microprocessor (one or more), a graphics processing unit (one or more), a memory (one or more), a permanent or non- transitory computer-readable storage medium (one or more), a permanent or non-transitory computer-readable storage memory (one or more), a server (one or more), a database (one or more), a computer program (one or more), a network (one or more), a framework (one or more), a neural network (one or more), or others. The device may comprise a physical input interface (one or more), a physical output interface (one or more), a man-machine interface (one or more), a user interface (one or more) or others. Computer programs according to the invention may comprise instructions for executing the steps and operations of the method. The computer program may be recorded on any storage media, which may be permanent storage memory, i.e. a non-transitory computer-readable storage medium, or a non-permanent storage memory.

The calculator CAL may output the difference of fill ratio on a physical output 11 (interface, user interface, screen or other) to which it is connected.

SAR means “Synthetic Aperture Radar”. The SAR images li have been acquired from one or more satellites 100. SAR satellite constellations, such as COSMO-SkyMed, TerraSAR-X or Sentinel-1 or other, may be used.

In step S15, the calculator CAL aligns the SAR images by a co-registration step.

Since the storage tanks often have large dimensions, they can be imaged on images from Synthetic Aperture Radar (SAR) satellite constellations. Bright pixels can be detected for each tank in the amplitude c/Zj of a SAR image in slant range geometry, corresponding to corner reflections A, B on the tanks.

When looking at a co-registered time series of SAR images, some bright pixels remain in the same position. Whereas for floating roof tanks, the bright spot corresponding to the floating roof corner moves by a few pixels from one date to another, corresponding to a height change. A method (such as document US 10,672,139 B2) were previously developed to convert the floating roof column index at a certain date into a crude oil volume (or a normalized "fill ratio" e [0,1]) for the storage tank.

On the other hand, for fixed roof tanks, there is no movement of bright spots in the SAR image and it is not possible to deduce the crude oil volume directly with the previous technique. According to the invention, Interferometric Synthetic Aperture Radar (InSAR) techniques may be applied on fixed roof tanks. A change in oil quantity inside the storage tank induces a small deformation of the tank. InSAR is applied to measure small phase changes of the reflectors located on the tank, since these measurements can be later used as a proxy of the change of oil volume inside the tank.

In step S20, the calculator CAL identifies in each SAR image I, a point A corresponding to an upper reflector 40 located on the circumference 20 of the tank 1 (floating roof tank 1 or fixed roof tank 1 ).

The upper reflector 40 is situated on the outside of the wall 2 and may surround the wall 2. The upper reflector 40 is situated near the roof 3 and may surround the roof 3. The upper reflector 40 may be metallic. The upper reflector 40 is situated at a distance (not null) from the bottom 21 of the wall 2 and is situated closer from the top 22 of the wall 2 than from the bottom 21 of the wall 2.

The upper reflector 40 may be or comprise the upper platform 4 located on the circumference 20 of the tank 1 or another upper reflector located on the circumference 20 of the tank 1 . The upper platform 4 may be a walking platform on which a user can walk.

The upper reflector 40 produces bright (or brighter) pixels for each tank 1 in the amplitude c/Zj of each SAR image I, in slant range geometry, corresponding to corner reflections on the tank 1 as shown by the arrows on Figure 1. Consequently, the upper reflector 40 may be detected as first point A from the other pixels in each SAR image I,. The point A may be the corner formed between the upper platform 40 on top of the tank 1 and the tank facade 2, i.e., the fixed roof corner A.

Also, the bottom 21 of the wall 2 of the tank produces bright (or brighter) pixels for each tank 1. The bottom 21 of the wall 2 of the tank may be detected as second point B from the other pixels in each SAR image I,. The point B may be the corner formed between the tank facade 40 and its base 41 , i.e., the fixed base corner.

The bright pixels for points A and B may appear aligned on the same image row with increasing column index according to range (distance to the satellite) or may be other.

The first point A and/or the second point B are the at least one point Q, as mentioned below.

The point A may be detected by one or more pixels among the pixels of the SAR image li. In case that there is a change of height of the liquid stored in the tank 1 (floating roof tank 1 or fixed roof tank 1 ) from one SAR image I, to the other, this change of height will cause different deformations of the wall 2 and a different deformation of the upper reflector 40 or upper platform 4 (for point A) and/or of point B.

As shown on Figures 3 and 4, a SAR sensor 101 (on a satellite 100) is an active system, i.e. the antenna of the SAR sensor 101 sends an electromagnetic wave signal and then records the signal of echoes reflected by the tank 1 situated on earth. The SAR sensor 101 produces a SAR image li from the signal of echoes received. The signal of echoes comprises the amplitude and phase for each pixel of the SAR image. The amplitude and phase of the echoes received by the SAR sensor 101 are stored in each pixel of the SAR image li as a complex number made of a real part and of an imaginary part. The phase of each SAR image li contains information about the path delay.

The pixels of tank 1 are called target M or M’ below. Looking at Figure 3 and considering that the satellite 100 is located at primary position P (P for primary acquisition), the phase <I>P(M) acquired by the SAR sensor 101 for target M reflecting the signal back to the SAR sensor can be written as:

where PM is the distance between primary position P of the satellite and target M, A is the carrier wavelength of the signal, 4n is related to the two-way path, <X>_SCatt,p(M) is the phase shift due to the interaction of the microwaves of the signal and the target M.

Then the SAR sensor 101 performs another acquisition at secondary position S (S for secondary acquisition) of the satellite 100. If due to the deformation of the tank 1 the target M has moved to position M’, the phase d>s(M’ ) acquired by the SAR sensor 101 for target M’ can also be written as:

where SM’ is the distance between secondary position S of the satellite and target M’, and <I>_Scatt,s(M’) is the phase shift due to the interaction of the microwaves of the signal and the target M’ .

Two acquisitions of the SAR sensor 101 at position P and S acquire the amplitude and the phase of the target at positions M and M’.

The figures 3 and 4 are not to scale, but have been exaggerated for visualization purposes. The distance M to M’ (deformation of the tank, for example between 0 and 1 meter) 1 is actually many orders of magnitude smaller than the satellite-target distance (distance PM or SM’ of about 800 km). The distance PS between the two positions P and S of the satellite (approximately of several hundred of meters) is actually many orders of magnitude smaller than the satellite-target distance and many orders of magnitude longer than the distance M to M’. Consequently, the direction PM can be considered as being approximately parallel to direction SM’

In step S30, the calculator CAL CAL selects K couples (IP, IS) of SAR images among the time series of SAR images I,. Consequently, P is an integer higher than or equal to 1 and lower than or equal to N+1 . Consequently, S is an integer higher than or equal to 1 and lower than or equal to N+1 . P is different from S. In step S40, the calculator CAL calculates for each couple (IP, IS) of SAR images an interferogram AO,_nt . The interferogram AO,_nt is an image containing the difference of phases modulo 2n of a primary image and of a secondary image, which has been previously aligned with the primary image. The interferogram AO,_nt is calculated for all the pixels of the SAR images (IP, IS).

From the above equations, the interferogram AO,_nt for each couple (IP, IS) of SAR images is equal to

wherein tp is the acquisition time of SAR image IP, and ts is the acquisition time of SAR image Is.

It follows from the above :

The first term <p_topo(t_P, t_s, M) = - y(PM - SM) is a topographic phase difference.

The second term is a deformation phase.

Since the calculator CAL has access to the satellite orbits (time, position and velocity of the satellite during the acquisition) and to image acquisition timing from the metadata of the SAR images, a simulation <I>topo,simu(tp, ts, M) of the topographic phase Oto_Po,simu(tp, ts, M) using a Digital Elevation Model (DEM) is performed in step S50 described below to compensate this term <p_topo(t_P, t_s, M) from the interferogram Ad>_int(tp, ts, M) to calculate the differential interferogram AODmt = A<X>D_int(tp, ts, M) in step S60 as follows :

This is the basis of the Differential InSAR(DlnSAR) technique (or DlnSAR as mentioned on Figure 2), to measure the deformation depending on the differential interferogram AODjnt(tp, ts, M). Looking at the deformation phase </)_ciefo(tp> ts> M') = - y (SM - SM'), the deformation is measured in the Line Of Sight (LOS) to the satellite, i.e. in the direction defined by the line joining the satellite to the target. Since there are K couples of SAR images (IP, IS), there are K differential interferograms

AODint-

In step S50, the calculator estimates the topographic phase difference Oto_Po,simu(tp, ts, M) of for each couple (IP, IS) of SAR images from a digital elevation model DEM, and from data contained in the metadata :satellite orbits during the acquisitions and image timing information.

In step S50, the calculator CAL may choose locations having each a longitude Long, a latitude Lat and an altitude h in the digital elevation model DEM covering the SAR image I,. The calculator CAL may determine from the metadata a geographical extent (minimum/maximum longitude, minimum/maximum latitude) of the SAR image I,. The calculator 101 may download the digital elevation model DEM covering a geographical extent of the SAR image I,. The calculator 101 may download the digital elevation model DEM using external tools. The digital elevation model DEM gives for each of the locations of the digital elevation model DEM the longitude Long of the location, the latitude LAT of the location and the altitude h of the location. The primary coordinates R, COL of one of the pixels in the SAR image I, can be computed by a geolocation operation (in this case, a projection operation). The digital elevation model DEM gives the longitude, latitude and altitude of some or all points of a surface of the ground surface area covered by the SAR image I,. The locations Q are chosen among the points of this surface of the digital elevation model DEM.

The steps S20, S30, S40 and S50 are grouped in a step of image graph generation in figure 2.

In step S60, the calculator CAL calculates the differential interferogram A<X>D_int(tp, ts, M) by subtracting the topographic phase difference <X>to_Po,simu(tp, ts, M) from the interferogram for each couple (IP, IS) of SAR images, i.e.

Consequently, the calculator CAL calculates K differential interferograms AODmt = AOj_nt(tp, ts, M) respectively for the K couples of SAR images (IP, IS). In step S70, the calculator CAL estimates the difference of fill ratio of the tank from the K differential interferograms AODmt = A<Wi_nt(tp, ts, M) taken at the at least one point Q.

In step S70, the calculator CAL may estimate the deformation of the tank from the K differential interferograms AODmt between the images IP, Is. of each couple at point Q. The deformation may be calculated from

In an embodiment, as shown on Figure 4, the K couples (IP, IS) of SAR images are the N couples of SAR images which are successive in the time series of SAR images. Consequently, K = N, P=i and S=i+1. The K couples (IP, IS) of SAR images are the N couples (l„ IM) of SAR images I, and IM for i going from 1 to N. The K couples (IP, IS) of SAR images are the N couples of successive (or consecutive) SAR images (l„ IM ) in the time series of SAR images for i going from 1 to N.

In the following the method and device are described for this embodiment.

Consequently, IP = I, and Is = IM .

Consequently, AP = A, and As = AM .

Consequently, the differential interferogram AODmt is noted

for each couple of SAR images (l„ IM) for i going from 1 to N.

Below are described embodiments of the step S70 of estimating the difference of fill ratio of the tank from the K differential interferograms AODmt = <Pi taken at the at least one point Q, as illustrated on Figures 7 and 8.

According to an embodiment of the invention, in step S70, the calculator CAL chooses a search window SW of pixels surrounding a prescribed guard window GW containing the tank in each SAR image L in sub-step S71 , as shown in figure 6.

In sub-step S72, the calculator CAL calculates for each pixel Pc of the search window SW an amplitude dispersion DA(PC) equal to

wherein m_rA(P_c) is the mean of the amplitude values c/Zj(P_c) for the pixel Pc taken along the times series of N+1 SAR images, <T^(P_C) ^1S the standard deviation of the amplitude values c/Zj(P_c) for the pixel Pc taken along the times series of N+1 SAR images.

In sub-step S73, the calculator CAL determines stable pixels U as the pixels Pc, which are situated in the search window SW and outside the prescribed guard window GW and whose amplitude dispersion DA(PC) is lower than or equal to a prescribed amplitude dispersion threshold T.

If no stable pixel U is found or the amount of stable pixels U is too small, the size of the search window SW is increased until a sufficient amount of stable pixels U is found. The set of stable pixels U found in this window SW is defined S.

In sub-step S74, the calculator CAL calculates a reference phase <p_ref,i for the stable pixels U for each SAR image I, and for each interferogram.

In step S70, the calculator CAL estimates the difference of fill ratio of the tank is made based on a double-phase difference 6<pi(Q) taken at the at least one point Q, which is equal to the differential interferograms A<Wi_nt= cp<(Q) taken at the at least one point Q to which the reference phase <p_ref,i is subtracted for each SAR image I, for i going from 1 to N, i.e. 6<pi(Q) = cpi(Q) - <p_ref,i in sub-step S75. The double-phase difference 6<pi(Q) is related to a difference in fill ratio of the tank between two dates.

Performing the reference phase (p_ref_:i computation in the search window SW and outside the guard window GW mitigates the influence of the tank filling on the reference phase. Indeed, the phase of the pixels on the tank or those that are very close to the tank might be affected by the tank filling. The reference phase is supposed to represent the phase of the stable background pixels.

In sub-step S73, the calculator CAL may calculate the reference phase <p_ref,i for the stable pixels U for each SAR image L as being the differential interferogram <p, calculated for the pixels Pc whose amplitude dispersion DA(PC) is the lowest. , i.e. (p_ref_:i = <?[(&) of the most stable pixel : U having the lowest amplitude dispersion: U = argmin_UeS£>^(U)), wherein at least one stable pixel U is needed.

In sub-step S73, the calculator CAL may calculate the reference phase <p_ref,i for the stable pixels U for each SAR image I, as being an angle being a weighted average phase as follows :

wherein z. is an angle of £ ' UCV <^JJI

in the]-n,n] interval,

V designates the stable pixels U, a>t(U) is a weight for each stable pixel U, <p, is the differential interferogram calculated between the SAR image L of the time series and the SAR image l,_+i of the time series for i going from 1 to N.

In an embodiment, (^(LT) = 1.

In another embodiment, the calculator CAL may calculate (^(LT) depending on the amplitude value c/Zj(U) of the stable pixel U in SAR image I, and on the amplitude value c/Zj₊₁(U) of the stable pixel U in SAR image IM . For example, a)t(PS) =

Or, _in another example, a)_i(U') = c/Z;(U) ■ c/Z;₊₁(U). Or, in another example, a)_i(U')

In sub-step S73, the calculator CAL may calculate the reference phase <p_ref,i for the stable pixels U for each SAR image I, as being an affine phase <p_re/,i(R) = ax_R + j3y_R + c for a point R of coordinates (x_R,y_R) for point R = Q in the SAR image l„ where a, ft and c are parameters fitted on the stable pixels U belonging to V, using the periodogram technique, with

where

and c = Z-Y^a,^ wherein z. is an angle

Nb is a number of stable pixels U in the search window SW and outside the prescribed guard window GW. At least 20 stable pixels U may be needed. The non-linear optimization can be solved using recursive-exhaustive search in the ]-n,n] interval (other optimization schemes can be adopted near the peak).

In sub-step S75, the calculator CAL may calculate the double-phase difference 6<pi(Q) taken at the at least one point Q between the differential interferograms cpi(Q) taken at the at least one point Q and the reference phase <p_ref,i for each SAR image I, for i going from 1 to N is equal to

MQ) = W(<Pi(Q) - (p_ref,t) where W(<p) = <p - is a wrap operator that transforms a phase <p into the ]-

n,n] interval. It is called double-phase difference 6<pi(Q) because <pi is a temporal difference (recall that <pi is a simplified notation of the differential interferogram A<W,_nt between consecutive dates). Then, we are doing a spatial difference with reference to the reference phase deduced from the neighboring area.

In an embodiment of step S70, the calculator CAL calculates from each couple of SAR images l„ IM for point Q an orthogonal baseline

which is a projection of a line PS joining positions P, S of the satellite(s) 100 having acquired respectively the SAR images l„ IM on a direction D, which is orthogonal to the line of sight PM of point Q from position P defined for SAR image I, and which departs from this position P for SAR image l„ in sub-step S76. The projection of line PS is made along the line of sight PM.

In sub-step S77 of step S70, the calculator CAL estimates the first variable parameter

K(Q) being optimized into 7C(Q), which is equal to

AT(Q) = argmin

where p is the power of the L^p-norm minimized. This estimation is made by applying to the double-phase difference <5q>i(Q) a linear regression model, which follows <5</>i(Q) = K(Q) ■

which depends on the orthogonal baseline

and in which K(Q) is a first variable parameter and Vi(Q) is a second variable parameter, and by varying K(Q) and v,(Q) .

In sub-step S78 of step S70, the calculator CAL calculates an estimated residual double-phase difference Sq>i-_topo,res (.Q) equal to ficPi-topo.res G2) = K(Q) • BtW

In sub-step S79 of step S70, the calculator CAL calculates a compensated doublephase difference <5<Pi-comp(Q) by subtracting the estimated residual double-phase difference Sq>i-_topo,res (.Q) from the double-phase difference 6<pi(Q), as follows

Scpi-comp tQ) = ^(Pi (.Q) ~ Sq)i-topo,res (.Q) - This compensates the effect of topography.

The calculator CAL estimates in step S70 the difference of fill ratio of the tank based on the compensated double-phase difference <5<Pi-comp(Q) •

In an embodiment of step S70, the calculator CAL calculates the difference ^^ui, noisy (Q) of fill ratio of the tank as a function of the compensated double-phase difference 8<Pi-comp (Q) and of a slope

sub-step S80. The slope ?7(Q) is prescribed in a sub-step S81 or calculated in a sub-step S83.

The normalized tank fill ratio at date i for SAR image I, is

e [0,1]. For floating roof tanks, Ui can be deduced using the floating roof position (see the document US 10,672,139 B2, which is incorporated herein by way of reference). Thus, from a study we conducted on floating roof tanks for Q, we know that the deformation 8<pi__def₀ between consecutive dates is related to a change in tank fill ratio Au_t e [-1,1] through a linear relationship:

, where we will have in practice a constant for each point Q (i.e. r/_A for the roof and ri_B for the base). Suppose we have L measurements of a tank fill ratio difference Auj between some consecutive dates. They might be measured by the floating roof position (see the document US 10,672,139). They might be measured by in-situ sensors, which also works on fixed roof tanks; Another option would be to use another image acquisition modality, such as thermal images captured by drones, where the oil volume can be seen even in the fixed roof tanks. In any case, we will have the following linear system of L equations to be solved to find ??_R:

In another embodiment of step S70, called aggregation in sub-step S82, the at least one point Q is the first point A corresponding to the upper reflector 40 located on the circumference 20 of the tank 1 and the second point B corresponding to a bottom 21 of the tank.

The calculator CAL calculates in step S70 the K differential interferograms A<W,_nt for the first point A corresponding to the upper reflector 40 located on the circumference 20 of the tank 1 and based on the K differential interferograms (AODmt) calculated for the second point B corresponding to the bottom 21 of the tank.

In sub-step S82, the calculator CAL estimates the difference Au_iidenoised of fill ratio of the tank as

1

6Ui , denoised i 7777

wherein C is the set of points A, B, constituted by the first point A corresponding to the upper reflector 40 located on the circumference 20 of the tank 1 and by the second other point B corresponding to the bottom 21 of the tank, qji(Q) is a reliability weight, higher than 0 and lower than 1 , which is associated to each point Q of the set C of points (A, B),

and ?7(Q) is a prescribed or calculated slope ??(Q).

This reduce noise.

For example, ip;

In another example, ip; (Q) = SRC(Q), wherein

wherein Z is a rectangular zone of pixels P_t of SAR image l„ which is situated around the point Q and from which the pixels PL situated in the line LQ to which the point Q belongs and the pixels PC situated in the column COLQ to which the point Q belongs, as well as the point Q are removed as shown on Figure 9. c/Zj(<2) is the amplitude value for the point Q, c/Zj(Pt) is the amplitude value for the pixels Pt,

NP is the number of pixels in the zone Z.

Or ip; (Q) = — ~ — z may be a combination of — - — _? and of SRC(Q),.

(^(0)) OWQ))

In an embodiment of step S70, the calculator CAL estimates the third variable parameter 77 (Q) being optimized into slope 77 (Q) as follows.

The calculator CAL estimates in sub-step S83 as shown on figure 8, by applying to the compensated double-phase difference <5<p_7-Comp(Q) having been calculated for another tank another linear regression model, which follows Aitj = ritQ) • 8(pj__comp(Q) + f(Q), which depends on the compensated double-phase difference <5<p_7-Comp(Q) and in which ??(Q) is a third variable parameter and f(Q) is a fourth variable parameter, and by varying ??(Q), the third variable parameter 77 (Q) being optimized into slope 77 (Q) which is equal to

??(Q) = argmin p(Q)

I is a fifth integer going from 1 to L, wherein L is a sixth prescribed integer higher than or equal to 1 .

The compensated double-phase difference <5<p_7-Comp(Q) has been calculated by the sub-steps S71 , S72, S73, S74, S75, S76, S77, S78 and S79 described above applied to cp , c/Z₇,

Consequently, these sub-steps S71 , S72, S73, S74, S75, S76, S77, S78 and S79 described above applied to

and then sub-step S83 of slope estimation of Figure 8 are part of a method for measuring the slope. This method for measuring the slope relates the compensated double phase difference at the roof or the base of the tank to the difference of the fill ratio.

The calculator CAL calculates the difference Au_{i noisy} (Q) or Au_t .denoised of fill ratio of the tank as a function of the compensated double-phase difference <5<Pi-_COmp(Q) and of the slope ??(Q).

The calculator CAL may estimate the slope 17 (Q) based on the method carried out on the other tank having a floating roof in sub-step S83.

The calculator CAL may carry out the method carried out on the other tank having a floating roof to estimate the slope 77 (Q) and use this slope 77 (Q) for estimating the difference Aui, denoised of fill ^rati° of the tank having a fixed roof.

Instead of difference of fill ratio of the tank, the method, device, and computer program according to the invention may calculate the deformation of the tank corresponding to the difference of fill ratio of the tank.

The invention relates also to the device for measuring a difference of fill ratio of a tank over time having the calculator CAL, the device comprising the at least one calculator CAL configured to carry out the above-mentioned steps and/or sub-steps.

The invention relates also to the non-transitory computer-readable storage medium having instructions that, when executed by a processor, cause the processor to execute the above-mentioned steps and/or sub-steps.

The invention relates also to the computer program for measuring a difference of fill ratio of a tank, having instructions that, when executed by a processor, cause the processor to execute at least the following steps :. downloading S10 a time series of N+1 SAR images I, having pixels covering the tank and having been acquired by acquisitions made by a satellite from different positions P, S relatively to the tank, wherein N is a first integer higher than 2 and is prescribed, aligning S15 the SAR images h by co-registration, detecting S20 in each SAR images I, at least one point Q selected among a first point A corresponding to an upper reflector 40 located on a circumference 20 of the tank 1 and a second point B corresponding to a bottom 21 of the tank, selecting S30 K couples IP, IS of SAR images among the time series of SAR images, wherein K is a second integer and is prescribed, calculating S40 for each couple IP, IS of SAR images an interferogram AO,_nt - Os(As) , estimating S50 a topographic phase difference <X>to_Po,simu(tp, ts, M) of for each couple IP, Is of SAR images from a digital elevation model DEM, from satellite orbits during the acquisitions and from image timing information, calculating S60 a differential interferogram AODmt by subtracting the topographic phase difference <X>to_Po,simu(tp, ts, M) from the interferogram AO,_nt for each couple IP, IS of SAR images, estimating S70 the difference of fill ratio of the tank from the K differential interferograms AODint, <Pi taken at the at least one point Q. The computer program may have instructions that, when executed by a processor, cause the processor to execute the above-mentioned steps and/or sub-steps.

Of course, any step or sub-step may be provided without the others.

Of course, the aspects, embodiments, features, possibilities, steps, operations and examples of the disclosure mentioned above may be combined one with another or may be selected independently one from another.

Claims

1. A method for measuring a difference of fill ratio of a tank over time, comprising the following steps carried out by at least one calculator (CAL) : downloading (S10) a time series of N+1 SAR images (l,j having pixels covering the tank and having been acquired by acquisitions made by a satellite from different positions (P, S) relatively to the tank, wherein N is a first integer higher than 2 and is prescribed, aligning (S15) the SAR images (I,) by co-registration, detecting (S20) in each SAR images (l,j at least one point (Q) selected among a first point (A) corresponding to an upper reflector (40) located on a circumference (20) of the tank (1 ) and a second point (B) corresponding to a bottom (21 ) of the tank, selecting (S30) K couples (IP, IS) of SAR images among the time series of SAR images, wherein K is a second integer and is prescribed, calculating (S40) for each couple (IP, IS) of SAR images an interferogram (AOmt), estimating (S50) a topographic phase difference (C>topo,simu(tp, ts, M)) of for each couple (IP, Is) of SAR images from a digital elevation model (DEM), from satellite orbits during the acquisitions and from image timing information, calculating (S60) a differential interferogram (AODmt) by subtracting the topographic phase difference (C>topo,simu(tp, ts, M)) from the interferogram (AO,_nt) for each couple (IP, IS) of SAR images, estimating (S70) the difference of fill ratio of the tank from the K differential interferograms (AODmt, <Pi) taken at the at least one point (Q).

2. The method according to claim 1 , wherein the upper reflector (40) is an upper platform (4) located on the circumference (20) of the tank (1 ).

3. The method according to claim 1 or 2, wherein the at least one point (Q) is the first point (A) corresponding to the upper reflector (40) located on the circumference (20) of the tank (1 ) and the second point (B) corresponding to a bottom (21 ) of the tank, wherein estimating (S70) the difference of fill ratio of the tank is made from the K differential interferograms (A<I>Dint= <Pi)) calculated for the first point (A) corresponding to the upper reflector (40) located on the circumference (20) of the tank (1 ) and based on the K differential interferograms (A<DDi_nt= <Pi) calculated for the second point (B) corresponding to the bottom (21 ) of the tank.

4. The method according to claim any one of claims 1 to 3, wherein the differential interferogram (A<Wi_nt= <Pi) is calculated between the SAR image I, of the time series and the SAR image IM of the time series for i being a third integer going from 1 to N, wherein the SAR image IM is consecutive to the SAR image I, in the time series, the K couples (IP, IS) of SAR images are the N couples of SAR images l„ IM, wherein estimating (S70) the difference of fill ratio of the tank from the K differential interferograms (AODmt, <Pi) is made based on the differential interferograms (A<Wi_nt= <Pi) taken at the at least one point (Q) for each SAR image I, for i going from 1 to N.

5. The method according to claim 4, wherein the method comprises: choosing (S71 ) a search window (SW) of pixels surrounding a prescribed guard window (GW) containing the tank in each SAR image l„ calculating (S72) for each pixel Pc of the search window (SW) an amplitude dispersion DA(PC) equal to

wherein m,_A(P_c) is the mean of the amplitude values (c/Zj(P_c)) for the pixel Pc taken along the times series of N+1 SAR images,

^(Pc) ^1S the standard deviation of the amplitude values (c/Zj(P_c)) for the pixel Pc taken along the times series of N+1 SAR images, determining (S73) stable pixels (U) as the pixels Pc, which are situated in the search window (SW) and outside the prescribed guard window (GW) and whose amplitude dispersion D_rA(P_c) is lower than or equal to a prescribed amplitude dispersion threshold (T), calculating (S74) a reference phase (<p_ref,i) for the stable pixels (U) for each SAR image li, wherein estimating (S70) the difference of fill ratio of the tank is made based on (S75) a double-phase difference (6<pi(Q)) taken at the at least one point (Q), which is equal to the differential interferograms (A<Wi_nt= <Pi) taken at the at least one point (Q) to which the reference phase (<p_ref,i) is subtracted for each SAR image I, for i going from 1 to N.

6. The method according to claim 5, comprising calculating the reference phase <p_ref,i for the stable pixels (U) for each SAR image I, as being the differential interferogram (<p,j calculated for the pixels Pc whose amplitude dispersion is the lowest.

7. The method according to claim 5, comprising calculating the reference phase <p_ref,i for the stable pixels U for each SAR image I, as being an angle

wherein z. is an angle

V designates the stable pixels U,

6)t(U) is a weight for each stable pixel U,

<p, is the differential interferogram calculated between the SAR image h of the time series and the SAR image l,_+i of the time series for i going from 1 to N.

8. The method according to claim 7, wherein (^(LT) = 1.

9. The method according to claim 7, wherein (^(LT) is calculated depending on the amplitude value c/Zj(U) of the stable pixel U in SAR image I, and on the amplitude value c/Zj₊₁(U) of the stable pixel U in SAR image IM .

10. The method according to claim 9, wherein (L7) =

11. The method according to claim 9, wherein (^(LT) = c/Zj(U) c/Zj₊₁(U).

12. The method according to claim 9, wherein

13. The method according to claim 5, comprising calculating the reference phase <Pref,i for the stable pixels (U) for each SAR image I, as being an affine phase <p_re/,i(R) = ax_R + py_R + c for a point R of coordinates (x_K,y_K) in the SAR image l„ where a, ft and c are parameters fitted with

where

and c = Z-Y^a, ^ wherein z. is an angle

Nb is a number of stable pixels (U) in the search window (SW) and outside the prescribed guard window (GW).

14. The method according to any one of claims 5 to 13, wherein (S75) the double-phase difference 6<pi(Q) taken at the at least one point Q between the differential interferograms cpi(Q) taken at the at least one point Q and the reference phase <Pref,i for each SAR image I, for i going from 1 to N is equal to tyi(Q) = W(<Pi(Q) - <pref,t) where W(<p) = <p - is a wrap operator that transforms a phase <p into the ]-

n,n] interval.

15. The method according to claim 14, comprising calculating (S76) from each couple of SAR images l„ IM for point Q an orthogonal baseline which is a projection of a line (PS) joining positions (P, S) of at least one satellite (100) having acquired the SAR images l„ IM in a direction (D), which is orthogonal to a line of sight (PM) of point Q defined for SAR image I, and which departs from the position (P) for SAR image l„ wherein the projection of the line (PS) is made along the line of sight (PM), estimating (S77) by applying to the double-phase difference 6<pi(Q) a linear regression model, which follows 8(pi(Q) = K(Q) • Bi(Q) + v_t(Q), which depends on the orthogonal baseline and in which K(Q) is a first variable parameter and v,(Q) is a second variable parameter, and by varying K(Q) and v,(Q), the first variable parameter K(Q) being optimized into K(Q), which is equal to

K(Q) = argmin

MQ)

where p is the power of the L^p -norm minimized, calculating (S78) an estimated residual double-phase difference 8<pi__{topo res}(Q) equal to

calculating (S79) a compensated double-phase difference <5<Pi-_COmp(Q) by subtracting the estimated residual double-phase difference Sq>i-_topo,res(.Q) from the double-phase difference 6<pi(Q). wherein estimating (S70) the difference of fill ratio of the tank is made based on the compensated double-phase difference <5<Pi-_COmp(Q) •

16. The method according to claim 15, comprising calculating the difference ^u_{i noisy}(Q) of fill ratio of the tank as a function of the compensated double-phase difference <5<Pi-_COmp(Q) and of a prescribed or calculated slope

17. The method according to claim 15, wherein the at least one point (Q) is the first point (A) corresponding to the upper reflector (40) located on the circumference (20) of the tank (1 ) and the second point (B) corresponding to a bottom (21 ) of the tank, the K differential interferograms (AODmt) are calculated (S70) for the first point (A) corresponding to the upper reflector (40) located on the circumference (20) of the tank (1 ) and based on the K differential interferograms (AODmt) calculated for the second point (B) corresponding to the bottom (21 ) of the tank, the method comprising estimating (S82) the difference Auj .denoised of fill ratio of the tank as

1 v^-1

^^ui, denoised ~ 4^1 (Q) ' ^^ui, noisy (.Q')

wherein C is a set of points (A, B), constituted by the first point (A) corresponding to the upper reflector (40) located on the circumference (20) of the tank (1 ) and by the second point (B) corresponding to the bottom (21 ) of the tank, qji(Q) is a reliability weight, higher than 0 and lower than 1 , which is associated to each point (Q) of the set (C ) of points (A„ B),

and ?7(Q) is a prescribed or calculated slope ??(Q).

18. The method according to claim 17, wherein ipi

19. The method according to claim 17, wherein

Ipi G2) = SRCC.Q), wherein

NP is the number of pixels in the zone Z.

20. The method according to any one of claim 16 to 19, comprising estimating (S83), by applying to the compensated double-phase difference <5<p₇-comp(Q) having been calculated for another tank another linear regression model, which follows Auj = ?7(Q) • <5<p₇-comp(Q) + /((?), which depends on the compensated double-phase difference <5<p_7-Comp(Q) and in which ??(Q) is a third variable parameter and f(Q) is a fourth variable parameter, and by varying ??(Q), the third variable parameter ??(Q) being optimized into slope ??(Q) which is equal to

??(Q) = argmin p(Q)

I is a fifth integer going from 1 to L, wherein L is a sixth prescribed integer higher than or equal to 1 , calculating the difference of fill ratio of the tank as a function of the compensated double-phase difference <5<p_7-COmp(Q) and of the slope 77 (Q) .

21 . The method according to claim 20, wherein the slope ??(Q) is estimated based on the method carried out on the other tank having a floating roof.

22. The method according to claim 21 , wherein the slope ??(Q) estimated based on the method carried out on the other tank having a floating roof is used for estimating the difference Au; .denoised of fill ratio of the tank having a fixed roof.

23. A device for measuring a difference of fill ratio of a tank over time, comprising at least one calculator (CAL) configured to carry out the following steps: downloading (S10) a time series of N+1 SAR images (l,j having pixels covering the tank and having been acquired by acquisitions made by a satellite from different positions (P, S) relatively to the tank, wherein N is a first integer higher than 2 and is prescribed, aligning (S15) the SAR images (I,) by co-registration, detecting (S20) in each SAR images (l,j at least one point (Q) selected among a first point (A) corresponding to an upper reflector (40) located on a circumference (20) of the tank (1 ) and a second point (B) corresponding to a bottom (21 ) of the tank, selecting (S30) K couples (IP, IS) of SAR images among the time series of SAR images, wherein K is a second integer and is prescribed, calculating (S40) for each couple (IP, IS) of SAR images an interferogram (AOmt), estimating (S50) a topographic phase difference (C>topo,simu(tp, ts, M)) of for each couple (IP, Is) of SAR images from a digital elevation model (DEM), from satellite orbits during the acquisitions and from image timing information, calculating (S60) a differential interferogram (AODmt) by subtracting the topographic phase difference (C>topo,simu(tp, ts, M)) from the interferogram (AO,_nt) for each couple (IP, IS) of SAR images, estimating (S70) the difference of fill ratio of the tank from the K differential interferograms (AODmt, <Pi) taken at the at least one point (Q).

24. A non-transitory computer-readable storage medium having instructions that, when executed by a processor, cause the processor to execute steps for measuring a difference of fill ratio of a tank over time, wherein the steps comprise downloading (S10) a time series of N+1 SAR images (l,j having pixels covering the tank and having been acquired by acquisitions made by a satellite from different positions (P, S) relatively to the tank, wherein N is a first integer higher than 2 and is prescribed, aligning (S15) the SAR images (I,) by co-registration, detecting (S20) in each SAR images (l,j at least one point (Q) selected among a first point (A) corresponding to an upper reflector (40) located on a circumference (20) of the tank (1 ) and a second point (B) corresponding to a bottom (21 ) of the tank, selecting (S30) K couples (IP, IS) of SAR images among the time series of SAR images, wherein K is a second integer and is prescribed, calculating (S40) for each couple (IP, IS) of SAR images an interferogram (AOmt), estimating (S50) a topographic phase difference (C>topo,simu(tp, ts, M)) of for each couple (IP, Is) of SAR images from a digital elevation model (DEM), from satellite orbits during the acquisitions and from image timing information, calculating (S60) a differential interferogram (AODmt) by subtracting the topographic phase difference (C>topo,simu(tp, ts, M)) from the interferogram (AO,_nt) for each couple (IP, IS) of SAR images, estimating (S70) the difference of fill ratio of the tank from the K differential interferograms (AODmt, <Pi) taken at the at least one point (Q).