Antoine N André

Patrick Sandoz

Benjamin Mauzé

Maxime Jacquot

Guillaume J Laurent

Sensing One Nanometer over Ten Centimeters: A Micro-Encoded Target for Visual In-Plane Position Measurement

Keywords: Nano-positioning, Position encoding, Pseudo-random binary sequences, Phase measurements, Computer vision

At the nanoscale, accurate control and position sensing are recurrent issues conditioning advances of technologies and instruments in diverse fields of application. Recent advances were reported that combine optical microscopy and computer vision for non-contact and multi-axis position sensing of end-effectors with a high accuracy. This paper contributes to this domain by introducing a micro-encoded target aimed to high accurate ( , , ) measurements over an extended range. The micro-encoded target is designed to allow linear phase measurements, leading to high accuracy in position sensing, and unambiguous binary encoding leading to absolute position determination. In this way, ultimate performances and allowed measurement ranges are independent of each other and a current zone of observation can be registered with respect to the whole encoded area with an extra-large range-to-resolution ratio. Experiments reported in this paper demonstrate a precision below one nanometer in lateral position and of 4 • 10 -6 rad. in rotation with an encoded target of more than 11 × 11 cm 2 . The device demonstrated can be either inserted into experimental setups to allow accurate position control or used to calibrate in-plane displacements of one to three degrees of freedom actuators.

I. INTRODUCTION

Nanopositioning is basic to all areas of nanotechnology from manufacturing to metrology. Some applications like scanningprobe microscopy and nanoindentation only need very small movement controlled with a high resolution, while others such as wafer inspection and microrobotics require a high degree of accuracy over long ranges in several directions. As anticipated by NIST in 2008, emerging nanotechnology industries will require revolutionary metrology with higher resolution and accuracy than had previously been envisioned [START_REF]The programs of manufacturing engineering[END_REF].

At the macroscale, many instruments allow to measure position along several directions with high range-to-resolution ratios. For instance, coordinate measuring machines (CMM) are able to measure the geometry of physical objects by sensing discrete points along the X, Y, and Z axes, typically with micrometer precision over a meter. Laser trackers (LT) used to calibrate robots or to measure large objects up to tens of meters can achieve pose estimation with micrometers resolution by determining the positions of optical targets held against those objects. The range-to-resolution ratio of both instruments are generally better than 10 7 (see Fig. 1). Even at a larger scale, the global positioning system (GPS) has a relative high precision, pinpointing to within 30 centimeters over the earth's circumference leading to range-to-resolution ratio over 10 8 .

At the nanoscale, interferometers (IF) are considered and used as reference measuring systems. It's one of the only solutions to reach nanometer precisions over meter ranges but recent works involving optical microscopy and image processing have pushed the limits of multi-axis pose estimation at the millimeter, micrometer and nanometer scales [START_REF] Assa | Virtual visual servoing for multicamera pose estimation[END_REF]- [START_REF] Liu | Relative pose estimation for alignment of long cylindrical components based on microscopic vision[END_REF]. Although optical microscopy is a mature field, it continues to see exciting new developments. For instance, Zhao et al. proposed in 2019 an in-plane XYΘ pose estimation method that achieves a resolution of 9.35 nm by analyzing the intersection of regular sets of straight and circular lines micro-machined on a polar micro-structure [START_REF] Zhao | Nanoscale measurement with pattern recognition of an ultra-precision diamond machined polar microstructure[END_REF].

For the purpose of in-plane XYΘ pose estimation, visual measuring methods can be distinguished between feature-based and area-based methods. First methods utilize local features of the image, such as corners or edges, to retrieve displacements and can reach sub-micrometer resolution. Indeed, although several works have shown that local methods achieve resolution so far as one tenth of a pixel, the resolution in optical microscopy is limited by the diffraction of the light leading to ultimate resolution near 10 nm [START_REF] Liu | Locating end-effector tips in automated micromanipulation[END_REF]- [START_REF] Avci | Towards high-speed automated micromanipulation[END_REF]. Using a scanning electron microscope could be considered to improve the resolution of the image but at the detriment of the measurement range [START_REF] Shang | Vision-based nano robotic system for high-throughput non-embedded cell cutting[END_REF]. As can be seen on Fig. 1, range-to-resolution ratios of featured-based methods are below 10 5 . In addition, featured-based methods are sensitive to image sharpness and signal-to-noise ratio (SNR).

On the other hand, area-based methods consider a large area of the image and this results in averaging the local effects of sharpness and SNR. They are based on the existence of an image texture, either naturally present in the scene of interest or specifically inserted. Most of these methods are variations of the original correlation method and exploit the redundancy of spatial information to achieve resolution better than one thousandth of a pixel [START_REF] Li | Realtime in-plane displacements tracking of the precision positioning stage based on computer micro-vision[END_REF]- [START_REF] Kokorian | In-plane displacement detection with picometer accuracy on a conventional microscope[END_REF]. To increase the SNR and to

The authors are with FEMTO-ST Institute, Univ. Bourgogne Franche-Comté, UMR CNRS 6174, Besançon, France {antoine.andre, patrick.sandoz, benjamin.mauze, maxime.jacquot, guillaume.laurent}@femto-st.fr get calibrated measurements, some authors propose to use a periodic pattern attached to the object of interest. For instance, Sugiura implemented a Moiré fringe method to sense a micro probe position that reaches a resolution of 42 nm on a limited range [START_REF] Sugiura | On-chip method to measure mechanical characteristics of a single cell by using moiré fringe[END_REF]. Yamahata performed a differential phase measurement in the Fourier spectrum to characterize the displacement with a sub-nanometer resolution of a MEMS actuator [START_REF] Yamahata | Subnanometer translation of microelectromechanical systems measured by discrete fourier analysis of ccd images[END_REF].

However, phase-based methods are limited to displacements shorter than the pattern period. This issue of extending the allowed range of phase-based methods beyond the period has been addressed in several ways. The choice of an extended periodic frame makes the allowed displacement range independent of the field of observation of the camera but requires some kind of additional codification to remove ambiguities of an entire number of periods. A simple solution, inherited from interferometry, consists in using two fringe sets with different periods. In 2016, we proposed a 2D-interleaved twin-scale pattern and demonstrated a experimental resolution below 0.5 nm over an unambiguous range of 221 µm using direct phase measurements [START_REF] Guelpa | 2d visual micro-position measurement based on intertwined twin-scale patterns[END_REF]. Another common solution to enlarge the range is to combine a binary encoding to the periodicity. For that purpose, the use of pseudo-random binary sequences obtained by linear feedback shift register (LFSR) is very convenient [START_REF] Galeano | Position-referenced microscopy for live cell culture monitoring[END_REF]. The unique property of LFSR sequences is that every binary word of bits shares -1 bits with its immediate neighbors. Therefore, whatever the number of bits chosen, the LFSR sequence can be materialized on a single track position encoder. This technique was used by Stevenson and Jordan two decades ago with a Manchester code to design 2D position encoders [START_REF] Stevenson | Absolute position measurement using optical detection of coded patterns[END_REF]. The Manchester code was also used by Masa et al. in addition to regular gratings for direct phase measurement in their 2D absolute encoder [START_REF] Masa | Nanometric Resolution Absolute Position Encoder[END_REF]. Kim et al. use a similar method with a different bit-coding technique and achieved a 15 nm resolution from a smallest feature size and elementary period of 4 µm and 32 µm respectively [START_REF] Kim | Note: Nonlinearity error compensated absolute planar position measurement using a twodimensional phase-encoded binary grating[END_REF]. However, in Stevenson, Masa and Kim approaches, the shifts induced for differentiating 0 and 1 bit values result in non-linearities detrimental to resolution.

This paper proposes an absolute binary encoding method preserving the phase distribution of the periodic frame, and thus avoiding non-linearities, with a significant improvement of performances achieved. The proposed micro-encoded target and its associated processing is intended for XYΘ sensing over a range up to 11 cm with a sub-nanometer resolution.

Section II describes the design of the micro-encoded target and emphasizes on its phase preservation properties. Section III presents the image processing developed for in-plane position reconstruction and estimates the ultimate method performances through computer simulations. Section IV presents experimental results obtained and section V concludes on the performances demonstrated using the micro-encoded target.

II. MICRO-ENCODED TARGET DESIGN

Our initial aim was to demonstrate a 10 8 range-to-resolution ratio by achieving sub-nanometer precision combined with a 10 cm measurement range by means of a micro-encoded target. For this purpose and in regard to bibliography, the target has been designed to fulfill two complementary requirements: i) to present a periodic frame suitable for high-resolution phase measurements, ii) to allow the unambiguous identification of the consecutive periods to encode a range of 10 cm. A supplementary practical requirement is that the smallest features must be detected without requiring high magnification lenses that would reduce the available work-space because of shorter working distances.

The periodic frame designed is based on linear shift register sequences and has been thought in order to preserve phase continuity and linearity. Fig. 2 compares the code chosen with respect to the Manchester code, notably used by Stevenson and Jordan, and to the phase encoded scale of Kim et al [START_REF] Stevenson | Absolute position measurement using optical detection of coded patterns[END_REF], [START_REF] Kim | Note: Nonlinearity error compensated absolute planar position measurement using a twodimensional phase-encoded binary grating[END_REF]. The code that we used appears in Fig. 2.e and is based on a regular distribution of stripes that are grouped by sets of three periods. The binary value 1 is made of three present periods whereas, for the binary value 0 , the central period is missing. This choice maintains a high density in the frequency carrier without introducing phase non-linearities. Conversely, in Fig. 2.c and Fig. 2.d, the black areas are shifted in accordance to the binary value and this alters the local phase and introduces non-linearities. Following the proposed scheme, the encoded range depends on the number of bits and on the period size as shown in Table I. The final trade-off chosen is a period of 9 µm with a 12 bit encoding that leads to an encoded length of more than 11 cm and compatible with a 10× imaging lens.

The 2D code is obtained by reproducing the same code along the perpendicular direction. To avoid rotation ambiguities, we broke the /2 symmetry by removing on square at one corner. Then each intersection of a binary line with a binary column results in a coding cell of 9 spots as depicted in Fig. 3.a and 3.b. Furthermore, in the final LFSR sequence, the word of value 2 -1 is removed in order to have at least one zero in every set of n bits since the missing lines associated with the value 0 serves as embedded clock to discriminate coding cell frontiers.

In practice, the micro-encoded target has been realized in clean room with a commercial mask generator (Heidelberg DWL200) of which X,Y positions, displacements and alignment are continuously controlled by a set of three interferometers. The target obtained corresponds to the designed distribution of holes etched on a thin chromium layer deposited on a 12.7 × 12.7 cm 2 piece of glass with a thickness of 3 m. Fig. 4.a presents a partial view of the micro-encoded target realized as recorded experimentally with the setup presented in section IV. The size of the micro-encoded target can be reduced by dicing to fit various setup needs without performance alterations except the allowed range of measurement.

To summarize, the measurement principle consists in retrieving the coordinate ( , ) on the encoded target of the local image recorded with the camera as well as its orientation with respect to the lines and columns of the target. In regard to the phase measurements, the positions ( , ) correspond to absolute phases ( , ) respectively equal to 2 / and 2 / with the physical period of the encoded target.

III. POSE ESTIMATION AND COMPUTER CODE VALIDATION A. Image processing

Once an image of a small-sized zone of the micro-encoded target has been recorded, the aim of image processing is to determine with the best accuracy the localization of this zone within the whole encoded area of 11 × 11 cm 2 . The reconstruction of the ( , , ) coordinates relies on phase computations and digital decoding following several sequential steps as described in the following and in Fig. 4.

The first processing step consists in retrieving the phase related to the two sets of straight lines associated to the perpendicular directions of the target. This demultiplexing is based on band-pass filtering applied to the Fourier transform of the recorded image. The different spectral lobes visible in Fig. 4.b correspond to the spatial frequencies of the recorded image presented in Fig. 4.a. Since the initial image is real, the Fourier spectrum obtained is symmetrical with each lobe associated to a complex conjugate at the diagonal of the spectrum center. Then, a first maximum 1 is detected in the spectrum and a hyper-Gaussian filter centered on this particular frequency is applied. By selecting a single lobe, we break the symmetry and, after inverse Fourier transform, this results in a complex representation of the line associated to lobe 1 . The wrapped phase distribution is then directed given by the argument of the complex data provided by the inverse Fourier transform as depicted in Fig. 4. 1 for a sub-area. The same procedure is applied similarly to the frequency component 2 associated to the perpendicular direction (Fig. 4. 2 ). More details on this phase processing can be found in [START_REF] Galeano-Zea | Position encryption of extended surfaces for subpixel localization of small-sized fields of observation[END_REF].

The wrapped phases of Fig. 4. 1 and 4. 2 are used in different ways for complementary fine and coarse position determinations. The fine position determination consists in localizing the straight lines with respect to the pixel frame of the image independently from the LFSR absolute encoding. From the target design, we know that the wrapped phase of Fig. 4. corresponds to a continuous and linear distribution of phase associated to the set of parallel straight lines of each direction. Therefore, the wrapped phase can be unwrapped to obtain a phase plane (not shown, see [START_REF] Galeano-Zea | Position encryption of extended surfaces for subpixel localization of small-sized fields of observation[END_REF]). This operation corresponds to increment the order of the straight lines at each 2 jump to fit the continuous distribution of parallel lines. The phase planes obtained are fitted with a least square equation that allows the determination of the phase related to the central pixel of the image with a high resolution as well as the orientation of the lines represented by the red arrows in Fig. 4. 1 and 4. 2 .

At this stage, the order of the lines is not registered with respect to the whole target but is relative to the first image point considered in the unwrapping process. However, the sub-period phase obtained through this procedure is fully representative of the position of the lines with respect to the pixel frame of the camera. The equations of the phase planes retrieved can be expressed as given in Eq. 1 and 2.

1 ( , ) = 2 cos + 2 sin + / + 2 1 (1) 2 ( , ) = 2 cos ( + /2) + 2 sin ( + /2) + / + 2 2 (2)
with the line period in pixels. 1 , 2 stand for unknown entire numbers of lines/columns separating the unwrapped phase planes obtained at this stage from the expected absolute phase planes registered with respect to the encoded target origin. The correct association of the phase planes 1 and 2 with the target directions and , and with the absolute phases and , will be determined later from the position of the missing corner of the coding cells. The period in pixels is related to the physical distance between lines and columns by: = • Mag/PixelSize, where Mag is the optical magnification of the imaging system and PixelSize the physical size of pixels on the image sensor. In Eq. 1, the plane coefficients 2 cos / and 2 sin / correspond to the coordinates of the spectral lobe 1 in Fig. 4. while the lobe 2 is located at coordinates 2 cos ( + /2)/ and 2 sin ( + /2)/ . Since these phase plane coefficients differ only by the sine and cosine functions of , the orientation of the lines and columns within the image is directly obtained through an inverse tangent operation.

The aim of the coarse pose estimation step consists in determining the right line order through the decoding of the LFSR sequence present in the image. This is necessary to adjust the relative phase planes with the correct number of periods; i.e. to remove the 2 1 and 2 2 terms in Eq. 1 and 2, and thus register absolutely the zone observed within the entire target. For this purpose, among each set of three lines and columns, the one representative of the bit value as a function of its presence or absence has to be determined. This identification is based on the comparison of the relative intensities of consecutive lines and columns. The absent lines and columns associated to the 0 bit value lead to a minimum that allows the identification of the coding ones. Then the mean intensity of each coding line and column is compared to the background intensity to determine the associated binary value. In this process, the wrapped phase is used to localize the position of the white spots, of which phase is close to 0, with respect to the dark background of which phase is close to ± . This phase-guided local contrast determination makes the binary decoding very robust, notably against possible illumination non-uniformity across the field of observation. The missing spot aimed to remove /2 angular ambiguities is identified in the same way as shown in Fig. 4. 1 and 4. 2 . Finally, the unwrapped phase maps can be adjusted with the correct number of periods to be fully representative of the position of the observed zone in regard to the whole target as represented in Fig. 4. . The success of this absolute registration procedure requires the observation of at least one word of 12 bits; i.e. about 40 lines and columns of spots. However, these 40 lines and columns are not required to be present throughout the whole image since their binary value can be identified even when only several periods are visible in an image corner (Fig. 4. 1 ).

In this image processing, the knowledge that the physical period of the target is equal to 9 µm allows the conversion of phase data into actual positions and distances on the encoded target without any requirement for vision system calibration. The reliability of this procedure will be discussed in the result section.

A video illustrating the decoding process on a rotation example can be found as supplementary material. This process, coded in C++ has the advantage to be as close as theoretical results can give. Every function programmed comes with a test function that ensures the final position computed is equal (minus the numerical noise) to what theory predicts. Current ++ software offers a 5 Hz measurement rate on 1024 × 1024 pixel images and 30 Hz with 512 × 512 pixel images with a quantification depth of 12 bits. In terms of complexity, the process requires one direct 2D Fourier transform, two inverse ones and some additional ( ) computations for the decoding (with the total number of pixels). The method is easily shown to have the same ( log ) complexity as the Fourier transforms.

B. Estimation of process performances

Computer simulations were carried out in order to determine the level of errors induced by round-off of the numerical computations in the absence of experimental disturbances such as spurious environmental displacements and imaging noise. For that purpose, we generated synthesis images of the target for a set of 200 random ( , , ) positions and compared the reconstructed positions with the expected ones. We chose an arbitrary period of 20.123 pixels which is close to the one obtained experimentally from the actual encoded target with a physical period of 9 µm. This choice allows the conversion of numerical errors in pixels into actual distances representative of the uncertainty level. Typical residual errors observed are represented in Fig. 5. 1 , 1 , 1 for , , measures respectively. We note that the shape of the histograms of the errors values (Fig. 5. 2 , 2 , 2 ) is close to a Gaussian distribution of residual errors. Fig. 5. 1 and 2 show the location of the explored positions and angles with respect to the whole target and the ] -, ] interval respectively. The standard deviations of the errors observed are summarized in Table II and result from the propagation of round-off error in the discrete Fourier transforms and further processing. However, the extremely small values obtained, in the range of 10 -8 µm and of 10 -10 rad. demonstrate that, experimentally, the major obstacle to performances will result from setup imperfections and environmental disturbances.

IV. EXPERIMENTAL RESULTS

A. Large displacement measurements

The performances of the proposed method were demonstrated using different setups and, to avoid frame rate limitations due to software, images were saved on the hard disk drive and then post-processed. Fig. 6 presents a picture of the 6 degrees of freedom Cartesian robot used for large displacement demonstration. The robot is composed of 3 prismatic and 3 rotational joints with integrated sensors and home assembled from SmartAct nanopositioning stages (more details in [START_REF] Rauch | Smallest microhouse in the world, assembled on the facet of an optical fiber by origami and welded in the robotex nanofactory[END_REF]). The robot allows preliminary focus adjustments over the whole encoded target surface as well as 10 × 10 cm 2 lateral displacements. The target is imaged onto a CMOS camera (UI-3280CP-M with a 3.45 µm pixel size) through a 10× Mitutoyo lens (M plan Apo 10× N.A.0.28) and a imaging tube (Optem mini 1.5x 300 FL). With this hardware, a pixelic period of 19.5 pixels is expected and values retrieved by processing the recorded images software agree with this prediction. Fig. 7.a compares the ( , ) positions reconstructed with the vision system to the control positions sent to the robot and used as references to validate the measurements. Unfortunately the displacement range allowed by the robot is smaller than the 11 × 11 cm 2 surface of the encoded target. Nevertheless, these results demonstrate the ability of the method to perform absolute position measurements over a range of 10 × 10 cm 2 . Fig. 7.b shows the error observed between the reconstructed position and the control position sent to the robot. The line in red corresponds to a unique round-trip and its shape indicates that the error is strongly linked to the axis positions. These data could be used to calibrate the robot to improve its accuracy as it has been done in [START_REF] Tan | Accuracy quantification and improvement of serial micropositioning robots for in-plane motions[END_REF]. A video showing the position retrieving on a large scale X-Y translation example is provided as supplementary material.

B. Precision and linearity measurements

The serial assembly of several single-axis actuators results in an increased uncertainty in the final position actually reached by the target. Therefore, we explored angular performances as well as ultimate lateral performances on dedicated setups involving a single actuator of interest. These supplementary experiments were carried out in a metrology room with temperature and humidity regulations. In Fig. 8, a 2 × 2 cm 2 piece of the target is fixed onto a Nano-MET10 (Mad City Labs) piezo-electric linear actuator with a 10 µm range and an internal sensor resolution of 0.01 nm. A complementary Smaract SR-2013 rotation stage with a closed loop resolution of 4.4 • 10 -7 rad. was used for angular measurements as shown in supplementary material.

Fig. 9 presents typical noise figures obtained while the target is kept static whereas Fig. 10.a and 10.b present the linearity curves retrieved on a linear displacement of 9 µm and on a rotation of 2 rad. respectively. Table III summarizes the standard deviations of the reconstructed ( , , ) position with values on and smaller than one nanometer and of a few 10 -6 rad. for in-plane orientation. Fig. 10 presents the deviation of the linearity data from the fitted lines for displacements and angles respectively.

An excellent linearity is obtained in rotation with a coefficient of 0.999961, meaning that the linear accuracy error is of a few 10 -5 . The maximum distance to the fitting line is of 8.59 • 10 -4 rad. and corresponds to 0.00024% of the 2 range. Results are a bit less impressive in lateral position with a straight line coefficient of 1.011478 corresponding to a linear accuracy error of 1.15%. This linear accuracy error questions the actual value of the target period that may slightly differ from the expected value of 9 µm. This issue will be addressed in a near future by calibrating the period value as well as its variations across the encoded area by means of a traceable and certified metrological inspection instrument. The maximum distance to the mean straight line is of 39 nm in displacement and corresponds to a linearity of 0.4% of the 9 µm course performed in this test.

C. Resolution estimations

Resolution tests were also carried out and results obtained are presented in Figs. 11. For lateral resolution, we sent as input position to the piezoelectric actuator a square signal at a 5 Hz frequency with a 1 nm amplitude. Reconstructed positions are presented in Fig. 11.a. where the square shape of the signal is visible despite an important level of noise. This was confirmed by a Fourier transform applied to both the reconstructed positions and the data returned by the actuator internal sensor. Associated spectra are presented in Fig. 11.b where the expected spectral lobe is indeed present at the actual excitation frequency. The Fourier transform was applied to a 30 s sequence with some zones more noisy which explain the relatively small magnitude of this lobe in Fig. 11.b. As suggested by computer simulations, environmental disturbances, known to raise significant issues at the nanoscale, may constitute the primary cause of limitation of the resolution. Fig. 11.c presents the results obtained in rotation while an input control signal with an amplitude of 1.7 • 10 -5 rad. at a frequency of 2 Hz was applied to the rotation actuator. The expected angular variation is clearly visible and the signal-to-noise observed suggests that even smaller angles could be detected. The precision of 4 • 10 -6 rad. previously estimated suggests that an improved resolution in this range should be achievable.

Unfortunately, the actuators used is unable to performs reliable changes smaller than these values and we could not explore experimentally the ultimate resolution of the method.

D. Discussion of experimental performances

This section aims to discuss the lateral resolution obtained that is limited at about one nanometer with this experimental setup. Uncertainty in measurement can be influenced by many different factors. We identified several sources of uncertainty in measurement that drastically increase the level of measurement errors. We first explored environmental instabilities as a major source of disturbances (vibrations, thermal changes, acoustic noise, etc.). In this aim, the series of 5 megapixel images that led to Fig. 9 was sub-divided into two sets of 2.5 megapixel images made of the left and right halves of every image respectively. Both series of images obtained were thus exposed to the same environmental disturbances that can be removed by operating a differential measurement between the two halves. By applying this procedure, the standard deviation of the distance between the centers of the half-images dropped from 0.96 nm to 0.24 nm in X and from 0.77 nm to 0.19 nm in Y. This result indicates that environmental disturbances are responsible for more than 93% of the total variance of errors. Fig. 12 presents the errors of reconstructed positions observed through this differential measurement procedure in regard to initial results (blue). The second factor of uncertainty that we investigated is the noise level in recorded images. The encoded target is made of holes etched on a thin chromium layer deposited on a glass plate that are originally aimed to serve as photomask for UV photolithography. Chromium is opaque to UV light but partially transparent to the red light used for experiments and suited to the spectral response of a CMOS camera. Combined with light scattered by the backside of the glass plate and illuminated in reflection, this configuration results in a low contrast and a high level of noise in recorded images. The typical standard deviation of the background of images is about 110 gray levels; i.e. 3% of the 3560 dynamics of the recorded images. The introduction of a Gaussian noise with a magnitude of 3% of the dynamics in 100 synthesized images (kept at the same place as in case of Fig. 9) results in an error of 0.15 nm in X and 0.16 nm in Y which could explain up to 5% of the total experimental variance. The remaining sources of uncertainties such as the pattern fabrication, the optical distortions and the quantification are thus responsible for less than 2% of the variance of errors.

V. CONCLUSION AND PROSPECTS

This paper proposed a visual method for registering a micro-encoded target of large dimensions, i.e. 11 × 11 2 , that allows in-plane position measurements with a nanometer precision in displacement and of about 4 • 10 -6 rad. in orientation. The range-to-resolution ratio has been experimentally demonstrated to be at least 10 8 . Interferometry (IF) still provides a superior range-to-resolution ratio. However, IF is basically aimed to 1D linear measurements and multi-DoF measurements require cumbersome hardware. Furthermore, the measurement of angles with IF is indirect and restricted to short ranges. Therefore, the proposed method constitutes a relevant and convenient 3-DoF alternative to IF, especially for non-laboratory experiments, even if its ultimate range-to-resolution ratio remains smaller. As far as we know, there is no other 3-DoF measurement technique with such a resolution on such an extended measurement range.

The absolute binary encoding is based on LFSR sequences and allows unambiguous position recovery from a local image of less than 40 periods; i.e. 360 × 360 µm 2 or 0.001% of the encoded area. The technique chosen to discriminate between 0 and 1 binary values does not introduce phase distortions and this avoids phase non-linearities that would be detrimental to performances.

The encoded target can be scaled for a different dimensional range as long as the magnification of the vision system is adjusted accordingly. The encoded target can be inserted in experimental setups to allow visual position control, if necessary after size reduction to meet actual needs and workspace capabilities. It can also be used to characterize or calibrate the motion of one to three degrees of freedom actuators. The method is self-calibrating and uses the design period as size reference that is available in any image of the target. However, this period is subject to some unknown degree of uncertainty and a metrological certification of its actual value and variations will be carried out in a near future. Major prospects concern the extension of the method to out-of-plane position measurement, for instance by using an interferometric imaging system to convert z-axis motions into interferometric phase changes.

Fig. 1 .

 1 Fig. 1. State of the art comparison in term of range-to-resolution performances between standard measuring instruments and vision-based methods used in microscopy.

Fig. 2 .

 2 Fig. 2. Diverse methods for bit encoding: a) binary sequence; b) black and white code; c) Manchester code [19]; d) Kim et al.'s code [21]; e) phase preserving code proposed in this paper.

Fig. 3 .

 3 Fig. 3. Definition and instances of coding cells, (a) Coding cell configuration: Corners are independent of binary values with one corner switched to Off to remove /2 ambiguities. Central line and column are set On or Off in accordance with the current binary value. (b) Coding cells obtained for the 4 ( , ) binary combinations and for a single angular direction. Colored squares discriminate coding spots from non-coding ones.

Fig. 4 .

 4 Fig. 4. Pose estimation process, (a) acquired target image, (b) Modulus of the 2D discrete Fourier transform. Circles exhibit the two frequency peaks of the two directions of the target image. The application of a Gaussian filter and inverse Fourier transform, result in the wrapped phase maps ( 1 ) and ( 2 ) representative of the line position with respect to the image pixel frame. ( 1 ) presents the binary decoding necessary to identify the coding cells and the 2D LFSR sequences. ( 2 ) zooms on a single coding cell with colored squares discriminating corners from coding line and column. As a result, (e) shows the fine and absolute registration of the current image within the encoded target.

Fig. 5 .

 5 Fig. 5. Tested random positions ( 1 ), ( 1 ) , ( 1 ): Reconstructed ( , , ) positions versus expected ones for a set of 200 random positions distributed on the whole target and 2 rad. ( 2 ), ( 2 ), ( 2 ): corresponding histograms of deviations observed. ( 1 ) shows the tested random positions on the target and ( 2 ) shows the random angles used for the target simulation

Fig. 6 .

 6 Fig. 6. Experimental setup used for large displacement measurements.

Fig. 7 .

 7 Fig. 7. (a) Reconstructed positions from a 10×10 cm 2 displacement experiment. (b) Differences observed between reconstructed positions and input commands send to the robot.

Fig. 8 .

 8 Fig. 8. Experimental setup used for linear displacement characterization.

Fig. 9 .

 9 Fig. 9. Repetition of measurements with a static target.

Fig. 10 .

 10 Fig. 10. Linearity tests: ( 1 ) in translation = 1.01148 and ( 1 ) rotation = 0.999961 . Zooms in inserts ( 2 ) and ( 2 ) show the interval between most-distant lines. ( 3 ) and ( 3 ) show the distance between measured data and the regression fits.

Fig. 11 .

 11 Fig. 11. (a) 1 nm lateral pulse (5 Hz captured at 100 fps) sent during 1 second. (b) Fourier spectra of reconstructed and reference positions (30 s). (c) 1.7 • 10 -5 rad. pulse signal (2 Hz at 20 fps) over a 2 s duration.

Fig. 12 .

 12 Fig. 12. Relative position discrepancies observed between the left and right halves of images (in red) in regard to the initial results of Fig.9 (in blue): (a) Along X; (b) Along Y (c) Along .

TABLE I CONCEPTION

 I TABLE FOR TARGET SIZE CHOICE WITH A VARYING PERIOD SIZE ( : PHYSICAL PERIOD.)

	order for the LFSR generation	Code size	Resulting target size
			= 3µm	= 9µm	= 25µm
	3 bits	24	72 µm	216 µm	600 µm
	6 bits	192	1.15 mm	1.7 mm	4.8 mm
	8 bits	768	2.3 mm	6.9 mm	19.2 mm
	10 bits	3072	0.92 cm	2.76 cm	7.68 cm
	12 bits	12288	3.68 cm	11.06 cm	30.72 cm

TABLE III PRECISION

 III AND LINEARITY MEASUREMENTS (OVER EXCURSIONS OF 9 µM FOR AND AND OF 2 FOR )

	Axis	Precision (std)	Linearity	Accuracy error
		0.96 nm	0.4%	1.15%
		0.77 nm	-	-
		4.07 • 10 -6 rad	0.00024%	< 0.004%

ACKNOWLEDGMENT This work was supported by ROBOTEX (ANR-10-EQPX-44-01), by Cross-disciplinary Research (EIPHI) Graduate School (contract ANR-17-EURE-0002), the French Investissements d'Avenir program and the I-SITE Bourgogne Franche-Comté (BFC) project (contract ANR-15-IDEX-03). The encoded target was realized thanks to the RENATECH technological network and its FEMTO-ST facility MIMENTO. Authors acknowledge D. Raddenzati, L. Robert, M. Suarez and L. Gauthier-Manuel for technological and experimental assistance.


Olivier Maury
email: olivier.maury@ird.fr

Blaise Faugeras

Victor Restrepo

FASST : a fully age-size and space-time structured statistical model for the assessment of tuna populations

Keywords: Stock-assessment, Population dynamic, Recruitment, Tuna fisheries, Fishery statistics, Fishing effort, Size distribution, Tagging, Growth curves, Migrations, Otoliths, Catch-effort, Catchability, Fishing mortality, Gear selectivity, Stochastic models, Mathematical models, Simulation, Body Size

come   L'archive ouverte pluridisciplinaire

RÉSUMÉ

Le présent document décrit le modèle FASST. FASST est un modèle d'évaluation du stock totalement structuré par âge, taille, espace et temps. Il est exprimé comme un système déterministe d'équations différentielles partielles continues en âge, en taille et en temps. L'espace est représenté à l'aide de grandes zones hétérogènes interconnectées par des taux d'échange. Tous les paramètres utilisés pour représenter le recrutement, la croissance, la mortalité et les déplacements sont choisis de sorte à être dépendants de la taille et indépendants de l'âge. Par conséquent, aux fins de simplicité, la dimension par âge peut être simplifiée. Le modèle est intégré numériquement afin de simuler la dynamique du stock. Un cadre complexe d'estimation d'erreurs de processus et d'observation est utilisé afin d'estimer les paramètres du modèle et donc l'état passé du stock, dans un contexte bayésien. Le modèle utilise simultanément la plupart des données sur les pêcheries actuellement disponibles pour l'évaluation des stocks de thonidés, telles que la prise par flottille en nombre et en poids, les échantillons de fréquence de taille, les analyses des otolithes et les données de marquagerécupération. Il fournit une nouvelle façon d'intégrer cette information dans un cadre cohérent d'évaluation du stock.

RESUMEN

Este documento describe el modelo FASST. El FASST es un modelo de evaluación de stock estructurado por edad, talla, espacio y tiempo. Se expresa como un sistema determinista de ecuaciones diferenciales parciales continuas en edad, talla y tiempo. El espacio se representa utilizando grandes zonas independientes interconectadas por tasas de intercambio. Todas las parametrizaciones utilizadas para representar el reclutamiento, el crecimiento, la mortalidad y los movimientos se eligen para que sean dependientes de la talla e independientes de la edad. Por consiguiente, para mayor sencillez, la dimensión edad puede ser simplificada. El modelo se integra numéricamente para simular la dinámica del stock. Se utiliza un marco complejo de estimación de errores de proceso y observación para estimar los parámetros del modelo, y por tanto, el estatus pasado del stock en un contexto Bayesiano. El modelo utiliza simultáneamente la mayoría de los datos de pesquerías disponibles actualmente para la evaluación de stocks de túnidos, como la captura por flota en números y peso, las muestras

Introduction

The goal of this paper is to present the FASST model (Fully Age-Size and Space-Time structured statistical model). FASST is a new model for tuna stock assessment which extends to a complex size-structured model the MULTIFAN-CL methodology presented in Fournier et al. (1998). Primarily designed for representing Atlantic bigeye tuna fisheries, the model is designed to be generic enough to be used for other species/regions. Tuna populations and their associated fisheries exhibit very specific characteristics which have to be taken into account to model realistically their dynamics:

• Tuna fisheries are highly heterogeneous in space and time and such heterogeneity has a high functional importance in their functioning. The mixing rate of fish between different regions is generally not total.

Then, the stock biomass located in a given fishing area interacts more or less strongly with the biomass located in other areas or with an unavailable part of the population located outside the fishing area. Important migrations and movements of fishes occurs at various scales so that the different regions interact differently: Fish movements must be explicitly represented either with a continuous model (such as an advection-diffusion-reaction model) either with a bulk transfer rate compartment model;

• Recruitment is continuous in time but highly heterogeneous in space and time (Cayré 1986) so that singular annual cohorts are most often difficult to identify. It follows that either continuous recruitment either a discretization of it such as the concept of mensual micro-cohorts [START_REF] Fonteneau | Considérations sur l'état du stock de listao de l'Atlantique[END_REF]) should be used;

• Non uniform mortality over sizes (due to both size dependant natural mortality and selectivities) may lead to non Gaussian distribution of size at age. Hence, simply adding a gaussian size distribution to an age structured model could lead to biases on both growth and mortality estimates (Figure 1): An explicit age and size structured model should be tested.

• Growth is potentially variable in space (especially for skipjack) depending of the region considered (Bard, 1986) so that fishes of the same age may exhibit very different sizes depending of their various history: A spatialized approach taking explicitly into account the potential variability of growth in space should be used;

• Many different fleets with different selectivity patterns and heterogeneous and changing catchabilities (in general showing an increasing efficiency) due to technical progress or to changes in fishing strategies and tactics fish the same population so that the use of fishing effort is problematic. A process error structure can profitably be added to fishing effort/mortality relationship and catchability time series (Fournier et al. 1998;[START_REF] Maury | Multi-fleet non-equilibrium production models including stock surface to estimate catchability trends and fishery dynamics in a bayesian context. Application to the skipjack tuna's fishery (Katsuwonus pelamis) in the Atlantic Ocean[END_REF]Maunder and Watter 2003).

• The uncertainty of quantities needed for management needs to be estimated so that a statistical approach must be used [START_REF] Hilborn | Quantitative Fisheries Stock Assessment. Choice, Dynamics and Uncertainty[END_REF].

Conducting reliable stock assessments taking into account those problems requires the development of a complex spatially explicit model dealing with all the fleets present in the fishery and with the movements of fish.

Such a stock assessment tool should profitably integrate a deterministic modelling of major processes with a statistical structure for both observation and process components [START_REF] Fournier | MULTIFAN a likelihood-based method for estimating growth parameters and age composition from multiple length frequency data sets illustrated using data for southern bluefin tuna (Thunnus maccoyii)[END_REF]1998). Space may be either considered using continuous representations such as advection-diffusion-reaction models [START_REF] Bertignac | A spatial population dynamics simulation model of tropical tunas using a habitat index based on environmental parameters[END_REF][START_REF] Sibert | An advection-diffusion-reaction model for the estimation of fish movement parameters from tagging data, with application to skipjack tuna (Katsuwonus pelamis)[END_REF][START_REF] Maury | Spatial modeling of Atlantic yellowfin tuna population dynamics: a habitat based advection-diffusion-reaction approach with application to the local overfishing study[END_REF]) either using discrete compartment models with fish transfer rates between zones (Fournier et al. 1998). Statistical models provide confidence intervals for the estimated parameters and inferred management policies [START_REF] Hilborn | Quantitative Fisheries Stock Assessment. Choice, Dynamics and Uncertainty[END_REF]Fournier et al. 1998;Maunder and Watters 2003). The parameter estimation may be conducted using maximum likelihood methods in a bayesian framework to easily handle the uncertainty concerning reference points useful for management [START_REF] Mcallister | Bayesian stock assessment using catch-age data and the sampling/importance resampling algorithm[END_REF][START_REF] Punt | Fisheries stock assessment and decision analysis: the bayesian approach[END_REF].

FASST is an age-size and time-space structured synthetic model of tuna population dynamics designed to use simultaneously all the information available (catches, effort, size frequencies, tagging3 , otolith increments, length/weight relationships, ...) in a Bayesian framework to estimate the posterior probability distribution of the parameters (recruitment, growth, movement, natural mortality and catchability). Such statistical approach has already proven to be very powerful (Fournier et al.1990 and1998;Maunder and Watters 2003) and should be an important improvement for tuna stock assessment in ICCAT.

2 Model structure

Notations

Variables

N, the fish number C, the catches R, the recruitment f, the fishing effort q, the catchability F, the fishing mortality γ, the growth rate D, the dispersion parameter of size Q, the proportion of fish in a given data set lying in a given length interval matur indicating maturity Θ, the percentage of females in a size class W, the body weight 2.1.2 Subscripts t, the time in months a, the age in months l, the length in cm j, the zone j k, the fleet k α, indicates the length frequency data set g, indicates the tag group (all the fish tagged the same month) ο, concerns otolith

Parameters

A, the maximal age n, the total number of zones T, the fish transfer rates between zones M, the natural mortality P, the fishing power s, the selectivity L∞, K, a 0 , X, the generalized von Bertalanfy growth parameters β, χ, δ, the original von Bertalanfy growth parameters θ, m, the parameters for the length/weight relationship i, relates length at maturity to asymptotic length Ε, is the total number of length-frequency samples Ω, is the total number of otolith samples σ, the standard errors

Random variables

ε, the catchability random walk error η, the fishing mortality/ effort error μ, the selectivity random walk error ω, the natural mortality random walk error

Population dynamics

The dynamics of the fish population is described through the use of fish density functions n j (l, t, a) whith j being the indice for spatial zones, ( )

T t , 0 ∈ being the time dimension, ( ) A a , 0 ∈
being the age dimension, and

( ) L l , 0 ∈
the length dimension. The number of fish N of age between a 1 and a 2 , length between l 1 and l 2 at time t in zone i is given by the integral:

( ) ∫ ∫ 2 1 2 1 , , a a l l i dlda a t l n .
A compartment model which takes explicitly into account movement between homogeneous zones, growth in length (and potential growth spatial variability), natural and fishing mortality is used for population dynamics. It is based on the following system of partial differential equations:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ = = ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ ∂ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ + ∂ ∂ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + - = ∂ ∂ + ∂ ∂ ∂ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ + ∂ ∂ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + - = ∂ ∂ + ∂ ∂ ∑ ∑ ∑ ∑ - = → - = → = → = → a t l n t f t l q dt a t l dc a t l n t f t l q dt a t l dc l l a t l n t l D l a t l n t l a t l n t l T l M t l F a t l n t l T a a t l n t a t l n l l a t l n t l D l a t l n t l a t l n t l T l M t l F a t l n t l T a a t l n t a t l n n k n k n k n k k k n n n n n n i i n n n i i n i n n n i i n i i i , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , , 1 , 1 , 1
, 1

1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 1 M M γ γ (1.)
with T i→j , the transfer rate from zone i to zone j, q i the catchability coefficient in zone i, M the natural mortality rate and f i,k the fishing effort exerted by fleet k in zone i.

Simplifying the model

In the model, according to basic ecological theory, all the parameterizations used to represent recruitment, growth, mortality and movements are chosen to be size-dependant and age-independent. Consequently, for simplicity purposes, the age dimension can be simplified without loosing anything from the complexity of the full model dynamics. With this simplification, the number of fish N of length between l 1 and l 2 at time t in zone i is given by the integral:

( ) ∫ 2 1 , l l i dl t l n .
The simplified model is written as follows:

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ = = ⎪ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎧ ∂ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ + ∂ ∂ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + - = ∂ ∂ ∂ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ∂ ∂ ∂ + ∂ ∂ - ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + + - = ∂ ∂ ∑ ∑ ∑ ∑ - = → - = → = → = → t l n t f t l q dt t l dc t l n t f t l q dt t l dc l l t l n t l D l t l n t l t l n t l T l M t l F t l n t l T t t l n l l t l n t l D l t l n t l t l n t l T l M t l F t l n t l T t t l n n k n k n k n k k k n n n n n n i i n n n i i n i n n i i n i i i , , , , , , , , , , , , , , , , , , , , , , , , , , , , , 1 , 1 , 1
, 1

1 1 1 1 1 1 1 1 1 2 1 1 2 1 1 M M γ γ ( 2.)
with T i→j , the transfer rate from zone i to zone j, q i the catchability coefficient in zone i, M the natural mortality rate and f i,k the fishing effort exerted by fleet k in zone i.

Spatial structure of the model

The geographical structure of the model is based on the heterogeneity of the fishery and the data available (Figure 2). Considering fleets to be zone-specific enables the model equations to be fully nested if zones are grouped. Then, different level of resolution (from a model with no spatial structure to a model with a fully detailed spatial structure) can be compared by using usual statistical tests (AIC, PBF, likelihood ratios) to select the resolution according to the amount of information contained into the data.

Recruitment and boundary conditions

Homogeneous Neumann boundary conditions at l=0 and l=L express the fact that the length of individuals can not reach negative values or values larger than L:

( ) ( ) 0 , , , , 0 = ∂ ∂ = ∂ ∂ l a t L n l a t n j j
The initial size distribution is fixed:

( ) ( ) a l n a l n j j , , 0 , 0 =
We also need a boundary condition for a=0. In the model, the only source of individuals is the recruitment which corresponds to the entrance of age 0 (in the full model) or size 0 (in the simplified model) individuals into the system.

Numerical approximation

The systems (1) and ( 2) are integrated numerically on a monthly time basis. The method of the characteristics is used for the full model (1). Explicit approximations are used for the growth, movements, mortality and catch components with a time splitting method using the Lie formula [START_REF] Marchuk | Splitting and alternating direction methods[END_REF]. A 2 nd order MUSCL (Monotone Upstream-centered Schemes for Conservation Laws) explicit numerical scheme is used for the growth component. Such a MUSCL scheme has been chosen because it ensures the positiveness of the solutions, it is conservative, very stable and has a much lower numerical dispersivity than simpler difference schemes.

In the full model (1), a closed boundary is used for the age dimension at age a=A, the maximal age which accumulates old fishes (A is a "+" group).

The recruitment do correspond with the boundary condition in age 0 and length 0. It is assumed to vary in each zone. A Beverton and Holt type stock-recruitment relationship is used in each zone with a lognormal process error allowing for random fluctuations of recruitment around its mean: Bigeye are assumed to be recruited when their fork length reach 30cm. In the full model, between age 0 and the age corresponding to the recruitment, the growth and movement processes are applied without mortality in order to obtain a realistic distribution of size at age when the micro-cohorts enter the fishery.

( ) ∫ ∫ = - = Θ ∝ + Ψ = = A a l l j l l t
Both models (1) and ( 2) are integrated. The simplified model ( 2) is used during the parameter estimation process to decrease calculation time and memory requirements. The full model ( 1) is used only for simulations, when the parameters have been estimated with (2) and when both age and size distributions are useful for interpreting the results.

Parameterizations

Growth

We consider that, in a given zone, fish growth follows a generalized von Bertalanfy curve which allows for slower growth rates for young fish as they are frequently observed for tunas, especially for yellowfin tuna (Pauly and Moreau 1997;Gascuel et al. 1992). From a process point of view 4 , this growth curve is based on a lengthdependent growth rate, independent of age:

l l l t l m l l catabolism anabolism t l l t w t j t j m t j t j t l j m t j t j t j t j t j t j , , 2 , 2 , , , , , 2 , , , , 3 2
3 3 3 δ β δ β γ χ δ β χ χ - ≈ - = ∂ ∂ = ≈ ≈ - = - = ∂ ∂ = ∂ ∂ - - (4.)
The parameter m can be estimated externally with a length-weight relationship. In the general case where χ=2  and m=3, equation (3) reduces to a von Bertalanfy growth equation with According with that formulation, age a fish belonging to the same cohort will have a different growth rate if they have a different size. This will cause a non-gaussian distribution of size frequencies for a given age 5 , even if mortality applies homogeneously over sizes.

Different assumptions can be made about the dispersion rate over the length dimension. They will have implications on the standard deviations of the fish length at age distributions. The a priori more realistic assumption is a dispersion coefficient being linearly related to growth rate:

v u D t a l j t a l j + = , , , , , ,
γ which also allows for growth-independent dispersion rate when v=0.

Mortality

In the model, the natural mortality M is length-dependant [START_REF] Hampton | Natural mortality rates in tropical tunas: size really does matter[END_REF]. Its structure is assumed to be a random walk with respect to length:

( ) ω σ ω σ ω ω , 0 2 1 2 N e M M l l l - + =
4 From a processes point of view, growth results from the antagonism between anabolism proportional to a physiological surface and catabolism proportional to fish weight. 5 If a more standard analysis based on the usual generalized von Bertalanffy growth curve ( )

( ) ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ - = Χ - Χ - ∞ 1 0 1 a a K a e L l
is desired in the full model, the growth rate may be assumed to be age-dependent and written as follows:

( ) ( ) ( ) Χ Χ - - Χ - - Χ - ∞ - = ∂ ∂ = 1 , , , , , 0 , 0 , 1 a a K a a K t j t j t a j t j t j e e K L a l γ
With homogeneous mortality over sizes, this will lead to a gaussian distribution of size frequencies for a given age.

The fishing mortality F is also length-dependent and defined as the sum of the fishing mortality of the n fleets k : According with the separability assumption, the catchability for each fleet is split into a length component (varying slowly over time), the selectivity s, and two time-varying components, the fleet fishing power p and the fish disponibility d for a given fleet:

∑ = = n k t l j k t l j F F 1 * , , , * , , with t j k t l j k t l j k f q F , ,
∑ ∑ = = k t j k t j k t k t l k k t j k t l j k t l j k f d p s f q F , , , , , , , , , , , , * , , ,
To take into account a potential seasonality due to fish behavior (reproductive concentrations, …), the disponibility d may be allowed to have a sinusoidal structure :

( ) [ ] [ ] [ ] π φ τ ς φ τ π ς 2 ; 0 12 ; 0 1 ; 1 , 2 sin 2 * , , , , , , * , , , , ∈ ∈ - ∈ ∈ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + + = + + R d d d d d j k r j k j d j k j k j k r j k t j k
With the positive part of a real number noted () + .

Catchability may also be related to external factors such as environmental factors (thermocline depth, ...) with a deterministic function.

Stochasticity is added to the deterministic F * parameterization to take into account time variability of the key parameters due to non-explicit external processes. A random walk structure is assumed for selectivity with respect to length:

( ) μ σ μ σ μ μ , 0 2 , , , 1 , 2 , , N e s s t l k t l k t l k - + =
Selectivities by size may be allowed to vary in time on a slow time basis (e.g. every 5 years) for each fleet.

To account for potential fluctuations of fishing power, we assume a random walk structure for the annual fishing power time series for each fleet. Such trends in fishing power may be due to changes in targeting or due to technological progress for instance. A normal random walk structure for log(p t ) allows the catchability to vary slowly over time without a priori assumption on its trend (increase or decrease) (Fournier et al., 1998;[START_REF] Maury | Multi-fleet non-equilibrium production models including stock surface to estimate catchability trends and fishery dynamics in a bayesian context. Application to the skipjack tuna's fishery (Katsuwonus pelamis) in the Atlantic Ocean[END_REF]. For that purpose, we assume that the yearly variability of p has the following structure:

( ) ε σ ε σ ε ε , 0 2 , 1 , 2 , N e p p t k t k t k - + =
To address high-frequency variability of the catchability coefficient, a lognormal process-error structure is assumed for the fishing mortality. Then, the fishing mortality of fleet k at time t is written

2 * , , , , , , 2 , , . η σ η
- = t j k e F F t l j k t l j k
where the η k,j,t are robustified normally-distributed random variables with mean 0 (Fournier et al., 1998).

Movements

Fish movement rates are restricted to occur only between spatially adjacent zones. They are supposed to have two components. A diffusive one which does not vary in time for a given adjacent pair of zones and an advective one which has a truncated sinusoidal structure which enables for seasonal advective movements. To account for length-dependent variation of fish speed, movement rates are supposed to be proportional to a power of fish length:

( ) [ ] [ ] [ ] π φ φ π φ π 2 ; 0 12 ; 0 1 ; 1 , 2 sin 2 sin 2 * , , , , , * , , , , , , , * , , , , ∈ ∈
Λ - ∈ Υ ∈ ⎪ ⎪ ⎪ ⎩ ⎪ ⎪ ⎪ ⎨ ⎧ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + Λ + Υ + = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + Λ + Υ + = + → → - → → → → → + → → → → → R T T l T T T l T T T l j i r l j i y l j i T l j i l i j r l i j t l i j y l j i T l j i l j i r l j i t l j i
Because movement patterns are likely to be different for adults which may experience reproductive migration and stay more in equatorial reproductive waters than juveniles which exhibit home range movements, two sets of movement parameters will be estimated and tested: a set for juveniles whose length is smaller than the length at first maturity l matur and a set for adults whose length is greater. Length at first maturity is defined as being equals to a fixed fraction of the asymptotic length L ∞ :

[ [ 1 ; 0 1 ∈ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = = - ∞ i i iL l
x m matur δ β

Tagging

Tagged fish are modelled by using exactly the same population dynamics equations and parameters that for non tagged fish. Additional mortality, tag shedding and tag reporting rates can be added if some information is available about it.

Fitting the model in a Bayesian context

To estimate the parameters in a Bayesian context, we will use the method of the maximum of posterior distribution (Bard 1974) by maximizing the sum of the log-likelihood of the data plus the log of the prior density function. Then, given the data, the Bayesian posterior distribution function for the model parameters has 8 components (one for the log-likelihood of the catch by fleet estimates L C , one for the length frequency distributions of catches L l , one for the age-length relationship from otolith reading L o , one for the catchability statistical structure L q , one for the recruitment deviations L R ,one for the tag recovery estimates L tag , one for the natural mortality random walk over length L M , and one for priors and penalties L penal ). Then, the posterior distribution is equal to L:

penal tag R M q o l C L L L L L L L L L × × × × × × × =
The parameters of the model will be estimated by finding the values of the parameters which minimize the negative log of this posterior distribution function. This minimization will be performed with a quasi-Newton numerical function minimizer using exact derivatives with respect to the model parameters with the AD model builder software (ADMB © 1993-1996 by Otter Research Ltd). ADMB calculates the exact derivatives with a technique named automatic differentiation (Griewank and Corliss 1991).

Catches component

We assume that the log of the predicted catches are the expected values of a random variable with a normal distribution:

( ) ( ) ( ) ∏∏∏∏ = = = = - - ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ = n j m k l l t t C C C t l k j C k C t l k j t l k j k e C L 1 1 0 0 2 log log , , , max max 2 2 , , , , , , 2 1 σ π σ
With C ˆ, the observed catches and C, the predicted catches. Important additional information may be provided by fixing the variances σ by fleet which fix the weights of the corresponding likelihood components.

Length-frequencies component

The robust normal likelihood of [START_REF] Fournier | MULTIFAN a likelihood-based method for estimating growth parameters and age composition from multiple length frequency data sets illustrated using data for southern bluefin tuna (Thunnus maccoyii)[END_REF]1998) modified by [START_REF] Maunder | A-SCALA : an Age-Structured Statistical Catch-At-Length Analysis for assessing Tuna Stocks in the Eastern Pacific Ocean[END_REF] is used for the contribution of the length frequency distribution to the likelihood:

( ) ∏∏ = = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - - ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + = E l l N Q Q I l

Age-length relationship from otoliths component

With a simple assumption of gaussian distribution of length measurement errors (and no error on age measures) the otoliths component is written as follows:

( ) ( ) ( ) ∏ Ω = - - ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ = 1 2 ˆ2 2 2 1 o a l a l o o o o o e L σ π σ
In the case where parameters can be estimated with the full model, a multinomial likelihood function can be used to model the error between predicted von Bertalanffy length at age and length at age derived from otolith increments counting (preliminary trials indicate that this robust function gives the same results that a likelihood assuming normal error for age estimations but brings information about potential spatial variability of growth).

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ∏ ∑ Ω = = = 1 0 , , , , , , , , o A a
a o t o l o f o z o a o t o l o f o z o C C L

Catchability structure components

This component combines the log-normal structured random walks for selectivities and fishing power trends for each fleets and the effort/fishing mortality process error which has a robustified normal structure. This robustified normal distribution assumes a probability p for unlikely events (events which are more than e times the variance from the mean) and 1-p for the standard normal distribution [START_REF] Sibert | Skipjack movement and fisheries interaction in the western Pacific. In Status of interaction of Pacific tuna fisheries[END_REF]:

( ) ( ) ∏∏ ∏∏∏ ∏∏ = = - = = = - = = - ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + + ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ - × ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ × ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎣ ⎡ =

Priors and penalties

Priors are added to the likelihood to take into account potential external information. A prior assumption based on the literature (e.g. [START_REF] Hampton | Natural mortality rates in tropical tunas: size really does matter[END_REF] will be made for the parameters M l . Another based on fishbase (Froese and Pauly 1997) data will be developed for the parameter i relating length at first maturity to asymptotic length.

  the * indicates that the statistical part of the fishing mortality has not yet been introduced).

  S α the sample size. N I is the number of length intervals in the sample.

  mortality random walk componentThis component corresponds to the log-normal structured random walks for natural mortality over lengths: deviations components This component assumes a log-normal recruitment deviation from its mean in each zone and time strata: Observed numbers of tag returns are related to predicted numbers of tag returns by a Poisson likelihood function[START_REF] Sibert | An advection-diffusion-reaction model for the estimation of fish movement parameters from tagging data, with application to skipjack tuna (Katsuwonus pelamis)[END_REF]):

  Penalties on fishing mortality (cf. multifan-cl and A-scala) may be added (at least during the first phases of the estimation process) to keep results in a realistic range.

Figure 1 .Figure 2 .

 12 Figure1. Simulated size frequencies obtained by using a gaussian size structure added to an age structured model (dashed line) or by using a fully size and age structured model such the one presented here (continuous line). First column, no fishing, homogeneous natural mortality over sizes. Second column, no fishing for sizes smaller than 80cm and strong fishing for longer sizes, homogeneous natural mortality over sizes. First line, absolute scale; second line, log scale. When the mortality is not homogeneously distributed over sizes (second column), it appears that simply adding a normal size structure on the age structured model may cause strong biases on abundance, position of the mode and shape of the size distribution for each cohort.

The reliable parameterization of such complex spatially explicit statistical models has proven to very difficult without using tagging data[START_REF] Hilborn | Determination of fish movement patterns from tag recoveries using maximum likelihood estimators[END_REF] Fournier et al. 1998;[START_REF] Sibert | An advection-diffusion-reaction model for the estimation of fish movement parameters from tagging data, with application to skipjack tuna (Katsuwonus pelamis)[END_REF] which bring to the model an important amount of information concerning growth and length-specific natural mortality[START_REF] Hampton | Natural mortality rates in tropical tunas: size really does matter[END_REF], catchabilities and selectivities[START_REF] Anganuzzi | Estimation of size selectivity and movement rates from markrecovery data[END_REF]) and age-dependent movements. Such reliable parameterization should be profitably conducted by using simultaneously fishery data and tagging data and potential auxiliary information for integrating the fishing effort standardization into the stock assessment model.