The Fukushima nuclear fallout. On March 11, 2011, after an earthquake and the tsunamis had devastated the northeast coast of Japan, a multiple core meltdown happened at the Fukushima Daiichi Nuclear Power plant, leading to a substantial release of radioactive material into the environment. Due to a general lack of detailed radiation measurement, SafeCast, a citizen science group based in Japan, started building instruments for geo-localized measurements of radioactivity. The ideas of project was to reconstruct a map of radioactivity from these measurements. Let us first consider the following simple, one-dimensional problem. Suppose that you're interested in a continuous time function, f of t, that is only known through some discreet set of measurements. The known samples here on the right side. In addition, these measurements might not be uniformly sampled. That is, they do not lie on the regular grid. We would like nonetheless, to be able to say something about the value of the function in areas where no data is available. If it possible to make a few assumptions about the function itself, it is impossible to make an educated guess about the value of the function at locations where the signal is unknown. This process is referred to as interpolation. In our example, we will consider the particular case of band limited functions. We can see here a band limited interpolation from the known samples. As you can see, the dotted line, or the continuous time function, goes through all the red samples so it does a perfect match at the known values. And it interpolates in between. Let us consider a continuous time, t-periodic, and band limited function f of t. Since it is band limited, this function can be parameterized by its Fourier expansion coefficients up to a finite order m. This means there are two N plus 1 free coefficients that are conveniently represented in vector form by vector c. So here we have the expression of the periodic function, which is the expression of the Fourier's series that we study in this class. And the vector is given here by 2 N plus 1 coefficients. We assume now that we know capital end pairs, ti, zi, where zi is a value of the function at time ti. We do not make any assumptions about the locations of the ti's. Given N large enough with respect to the maximum free coefficient order M, we can make a least squares fit of the data to find the free expansion coefficient that best approximates data collected. We need at least M equal 2M plus 1 measurements. But in general, we would like a much larger number because of possible measurement noise. The least square fit can be rewritten matrix form and solved efficiently. But let's just skip the details here. For a full demonstration, see the corresponding IPython notebook. Once we have recovered the best guess of the Fourier expansion coefficients, it is possible to resample the function f of t at any location we like. Often, the function is resampled on a regular grid. Let us give a first example. The function we chose here satisfies perfectly the assumptions of periodicity and band limitedness. It just the sum of two signs of frequency, one hertz and four hertz. The true function is represented by the thick black line. We pick at random 20 samples, represented here by red dots. The model order is N is equal to 6. The interpolation given by our least squares fit is a dotted red line. Since all the assumptions are satisfied, and because there is no noise, the interpolation recovers perfectly the original signal. You might be wondering what happens if the original signal does not satisfy some of the assumptions, as we might very well be the case in practice. For example, assume the single is not band-limited. We take here the example of signal, continuous everywhere, but not differentiable at zero. Such a signal is not band-limited. However, it's Fourier coefficients decay proportionally to the square of the frequency. That is quite fast. The non-band limitedness of the function can be somewhat compensated by picking a larger order for our model and having more samples. We choose here the ordinal to be m is equal to 10, and take n is equal to 20 samples, picked, again, at random. We can observe that our interpolation is still very close to the true function. The approximation is, however, not so tight around the sharp angle at t is equal to 0. Let us now move on to the real world example of Fukushima. We are interested in measuring radiation levels all over Japan, and beyond after the tsunami of 2011 and the nuclear accident. Instruments were loaned to volunteers in the field and measurements were carried out by fixing the sensors on cars. The sensors would then take a measurement of the radiation every five seconds. Therefore, sampling is generally limited to areas accessible by car. A subset of the measurements is shown here on the map of the Fukushima prefecture area. Each colored dot is a measurement with yellow being low intensity and red high intensity. You can observe some fairly high measurement levels on the path of the plume, starting around the plant, situated at the far right of the mountain and going northeast. This is a typical example where the location of measurements is constrained by available infrastructure. Since the measurements are done with a car, we can only measure on roads, as plainly visible on the map. It is thus necessary to apply interpolation to get an estimate of the intensity of radiation in areas where no samples are available, such as forests, fields, and so on. The interpolation is done using a model with m is equal to 441 Fourier coefficients and a number of samples that is close to a million. The function so obtained is then resampled on the regular grid and displayed here on the same map as before. A decent estimate of the total radiation field is obtained in this way. A few caveats are to mention. First, radiation is notoriously non-bandlimited, and we would probably need the larger order model for better approximation. Second, large measurement noise is present and also degrades the quality of the estimation. A few tricks were used in this reconstruction to improve the quality of the estimation. The code and detailed explanations are available in the companion IPython Notebook. Check it out.
