Single Plant Fertilization using a Robotic Platform in an Organic Cropping Environment

The increasing demand of organic vegetables has driven conventional farmers to change their businesses in order to develop organic methods for cultivation and the use of alternative techniques to avoid damaging the soil and the quality of the products. In some cases, to enhance biodiversity and soil fertility, crops are established in a mixed pattern called ‘strip cropping’ where single or dual lines of a given species are alternated with a second compatible species, with the aim to enhance resilience, system sustainability, local nutrient recycling, and soil carbon storage. However, this husbandry of crops grown and mixed in a strip design pose new challenges regarding mechanisation, which in many cases can only be overcome by increasing human tasks. To counteract the additional labour of a multi-crop system, one of the main objectives of the ‘Sureveg’ CORE Organic Cofund ERA-Net European project is to evaluate the benefits of growing in alternate rows for the production of organic vegetables and includes the use of robots as a tool to facilitate the automation of the process, allowing the individual treatment of organic fertilization at plant level. Within the project framework a modular proof-of-concept version has been produced, combining several sensing technologies (3x LiDAR, plus a multispectral RGB-NIR camera) with actuation in the form of a robotic arm operating upside down. The present work describes a method to develop fertilization tasks with recycled organic waste in strip-cropping farms, based on detection of plant species (cabbage) using and liquid fertilizer application with a robot. For that, it is necessary to identify the crop (variable in the rows) and place the robot with respect to the plants to apply the product. Test fields were located at ETSIAAB - UPM (40°26’33.1"N, 3°43’41.9"W) where stripped crops were established and monitored along the growing season. In order to detect autonomously each single plant, point clouds of the three LiDAR units were combined, soil was removed applying a weighted threshold, and plants were identified using clustering and convolutional neural network methodologies. To trigger the actuation system, the decision on which plant had to be sprayed with the liquid fertilizer was taken according to two factors: 1) the estimated volume of every single plant, and 2) the multispectral indexes calculated using the RGB-NIR camera. The prototype is fully functional, and further test are needed to quantify its performance.


Introduction
The increase in food consumption worldwide needs to be tackled by producers and technicians (Ritson 2020;Pachapur et al., 2020).However, due to conventional agricultural practices, the environmental effect associated to the use of fertilisers and herbicides has grown (Srivasta et al., 2020).Fortunately, precision agriculture allows to produce higher quality products, through sustainable development (Loures et al., 2020), the use of technological tools, sensory systems (Poblete-Echeverría et al., 2020;Singh et al., 2020), and modern actuation systems (Cubero et al., 2020;Moysiadis et al., 2020;Fue et al., 2020;Hussain et al., 2020).The Sureveg project (CORE Organic Cofund, 2019) focuses on the application of diversified strip-cropping systems to intensive organic vegetable cultivation, the reuse of biodegradable waste, as well as development of 31 automated machinery for the management of strip-cropping systems.For this purpose, a proof of concept was proposed in the form of a manually operated robotic prototype, containing 3 lidar sensors, a multispectral camera, and a 5 degrees-of-freedom (DOF) manipulator.A nozzle connected to a tank of organic fertiliser is used as the end-effector of the actuator.

Materials and Methods
The mobile platform implemented to the proof of concept of this work consists of different subsystems: the mobile structure with wheels, the robotic arm, the sensory system, the actuation system for the application of fertiliser, and the storage tank of liquid fertiliser.

Mobile Structure
The mobile structure that supports the robot and sensory systems -actuation, was built with aluminium bars (Bosch Rexroth 45 x 45) and four wheels, according to previously established requirements in the framework of an organic stripped crop.The frame supports the different elements and does not rely on an autonomous traction system.Figure 1 shows the implemented structure.In Figure 1, the robot is replaced by a dummy volume shown hanging down from the middle of the structure.This location and orientation allow it to reach the ground to acquire data and apply treatments.
Figure 1.Structure implemented with aluminium profiles (45x45) and wheels.

Robotic Arm
The robotic arm acquired for the prototype (Robolink Igus CPR RL-DC-5 STEP RL-D-RBT-5532S-BF) has the following characteristics: a load capacity of 1.5 kg, 5 degrees of freedom, IP54 protection, and able to follow point to point trajectories.The robot has a reach of 790 millimetres with a precision of 1 mm.Additionally, the robot weighs around 20 kg, including the electromechanical components and excluding the control and power elements.Its load capacity is 2.5 kg.

Sensory System
The sensory system mainly relies on two types of sensors: a multi-spectral camera and 3 lidars.Both provide information on vegetational health status and development.The multispectral camera is used to obtain multispectral images and compute vegetation indices of plants.This camera was installed on the centre of the frame of the cart, to take pictures of the target plants from the desired perspective.Additionally, 3 lidars provide 3D models of the soil and plants from multiple perspectives.Besides information on crop architecture, the combined 3D world model facilitates the motion planning of the robotic arm and the application of treatments in the adequate locations.Figure 2 specifies the sensor locations on the mobile structure.
• 3 • Figure 2. Frontal view (right) and rear of the platform (left), the yellow circles show the 3 lidar sensors and the red circle indicates the multispectral camera.

Sensing Algorithms
The sensing algorithms can be divided into two sub-categories: multi-spectral sensing using the multispectral camera (aimed at searching for differences in plant health) and 3D sensing using the LiDAR (focused on detecting growth volume differences).At the beginning of the crop season, the newly transplanted crops provide very little altitude difference when compared to the coarseness of the soil, which is why at the beginning of the season, only the use of the multi-spectral algorithms is recommended.

Fertilization System
The implemented actuation system consists of a tank containing a liquid treatment (Figure 3, highlighted in red), a hose from the tank to the robot tooltip, and a nozzle to spray the liquid (Figure 3, Indicated in yellow).In this way, the robotic arm can apply the liquid treatment in the desired position and orientation, as shown in Figure 4.The drive of the nozzle is carried out by means of the electronic control of a hydraulic motor.Test fields were located at ETSIAAB -UPM (40°26'33.1"N,3°43'41.9"W)where stripped crops of cabbages and leaks were established and monitored along the growing season.

Multi-Spectral Sensing
The Parrot Sequoia camera produces each image in fivefold: green, red, red-edge, near-infrared, and (high resolution) RGB.The time-lapse setting in combination with a forward speed generates a series of partly overlapping images that were be stitched together using Hugin Panorama, as it is open-source and allows for manual manipulation of automatically extracted features such that the exact same mosaicking operations can subsequently be applied to the other spectral bands.Each of the spectral images is obtained through a different lens, as can be seen on the bottom of the sensor shown on the left in ¡Error!No se encuentra el origen de la referencia..As a result, the images themselves show a slight translation with respect to one another, something that was assessed and corrected for automatically using e.g.imgregtform in Matlab.With the corrected mosaicked images, a series of vegetation indexes were calculated for the entire field.The most optimal vegetation index can depend on the crop type and variant, and can be optimised regarding the specific characteristics of the crops and soil spectra present in a certain field.In this work, results using NDVI are presented.
• 4 • The vegetation index was used to separate plant matter from soil, as well as assess the health of the plant.In Figure 4, the distribution of the vegetation index values in the entire row is shown in the form of a histogram, clearly revealing the three types of areas present in the image.In the distribution of NDVI values three gaussian distributions can be recognised: a wide one around low (negative) values for soil; a narrow one around 0 describing the stark shadows, i.,e.areas where all spectral reflectance values were very small; and a wide one of positive NDVI values describing the plant matter.The relative heights of these three peaks depends on the number of plants with regard to the amount of visible soil, as well as the lighting conditions.
To facilitate the identification of the optimal cut-off value to define plant matter, the vegetation index values were passed through a filter, further separating the peaks in the histogram, revealing a minimum between the plant matter pixels and the remaining pixels that was not present in the original distribution, as shown in Figure 5.Using Otsu's method, the cutoff value that minimizes the variance in each of the pixel populations (plant matter vs. non plant matter) was identified.
Figure 5 The same NDVI map from Figure 4 after passing through a Gaussian filter, yielding the histogram shown in the middle.On the right the final result is included.
After applying the Otsu procedure, the range of the filtered image no longer corresponds to the [-1,1] range of the original NDVI index.The altered values reveal a minimum between the plant matter and the non-plant matter, as obtained through Otsu's method, indicated in the histogram with a red vertical line in Figure 5.For each resulting cluster in the new image, a range of characteristics can be considered: e.g.vegetation index mean value, vegetation index distribution, cluster area, or perimeter.In this proof of concept, the mean value within a cluster was used.This value is then compared to the mean NDVI values of the other clusters within the same row to identify fertilisation necessities on a single crop scale.

3D Sensing
The 3D data was gathered using multiple SICK A.G. lidars.Each lidar uses an infrared rotating laser pulse, where the time of flight between emission and reflection in each rotational increment is used to calculate the distance of the sensor to any object.The reflectance intensity of the beam was found to not provide any additional information in this setting and was therefore disregarded.From the cylindrical coordinates, each distance was converted to a Cartesian coordinate system, where the inclination of the sensor and the odometry data of the cart's wheels were combined to provide the dimension perpendicular to the plane of rotation of the laser.This results in a 3D point cloud per lidar, that all describe a single crop row.They point clouds provide complementary information inherent to their varying angles and installation heights, while the majority of the surfaces is picked up by all lidars.This overlap was used to merge the point clouds into a high-resolution single point cloud, describing all surfaces within a crop row, i.e. both plant matter and non-plant matter such as soil.
As three-dimensional information is the only information present in this type of point cloud, i.e. no colour or reflectance information is retained, the points needed to be identified as crop or non-crop purely based on their location and height.• 5 • Especially for smaller plants, the coarseness of the terrain yields height data that is comparable to the height data of lower hanging leaves.This means an identification based on only height data was not possible.
To overcome this, each point was assigned a value  based on the heights of all points in within a 150 mm radius, where the points that are closest are taken into account more heavily than those at the border of this sphere, i.e.,  is a weighted sum of height squared over distance.The number of points that are present in this vicinity also contributes to this value, as higher surfaces such as plant leaves are closer to the sensor and are therefore sampled at a higher frequency than any points further away.
The resulting values  for all points within a point cloud can be sorted in ascending order.The texture predictability and constant distance of the soil assures that the values of  of the soil points yield a nice linear distribution.This linearity wa used to identify the value of  for a given point cloud where the sorted distribution deviates from this expectation.The points with a value of  higher than this threshold were considered plant points, and lower values were considered to describe non-plant matter.This algorithm allows for low plant leaves and protruding surfaces like rocks to be identified correctly, despite the rock points possibly having higher altitude values than those of the low leaf, as shown in Figure 6.
Figure 6 The point cloud with colour denoting height data (left) alongside the same cloud coloured by the values of .
The resulting 3D point cloud containing only plant points can subsequently be cleaned, where any singular noise points were removed, yielding only clusters of a certain size and up.For each of these clusters, a range of characteristics were calculated, e.g.covered soil surface area, cluster volume, cluster surface to volume ratio, or distribution of  values.These were compared to the other clusters within the same row to identify anomalies.

Combination
The multi-spectral sensing algorithm is able to detect any individual crops with either delayed size in terms of covered soil or altered vegetation index values.The 3D sensing algorithm identifies delayed development in terms of volume or height as well.This means a combination of both sensors was preferred to determine the fertilisation necessity of individual plants.
Besides the prescription map, the data of the 3D algorithm was used as a reference world model for the actuation algorithms, which serves both to locate the cart in the field as well as to plan the movements of the robotic arm as discussed below.
Drawbacks and potential pitfalls of these algorithms are the detectability of newly transplanted crops, as they might be missed by both sensors and therefore wouldn't get any treatment.As the field is prepared before transplantation and all individuals will need fertiliser this is not assumed to pose a big problem.The distance between plants is another aspect that needs to be considered, as plants that reach overlap or physically touch each other cannot be separated by either of the proposed sensing algorithms.This should be taken into account when assessing the results once the plants reach a certain growth stage where this starts to occur.Finally, the algorithms now assume all plant matter to belong to desired crops, i.e. weeds could mistakenly be identified as individuals with a delayed development when compared to the plants of the desired crop in that row.Crop type identification can be added to these algorithms in the future, but was not yet considered in this proof of concept.

Geometrical Characterisation
The clusters extracted in the previous section served as a reference to establish a path for the movements of the robotic arm to apply the liquid treatment.The extracted parameters are shown in Figure 7: in the form of centres and defined edges.Colours indicate each identified group.For the processing of the clusters, the unsupervised learning process Kmeans has been used, where the clusters are grouped with labels "cluster 1-10".
• 6 •   In Figure 10, a box plot shows the declining localisation error found at each of the positions indicated in Figure 9. Initially, the average error was around 12 mm; which reduces to approximately 0.1% at the end of a row.This initial positioning error did not affect fertilisation, because a margin of 5 cm was added to the radius defined in the geometrical parameters extraction to encompass the plant and define the passage zone of the arm's trajectory, thus avoiding the fertiliser being applied directly on the plant and damaging it.As the platform progresses, new points are added to the L-PC, allowing for more key points and improved localisation, resulting in errors with an average of 5 mm.

Conclusions
A proof-of-concept prototype has been built, based on a mobile platform, a sensing system composed of lidar sensors, a multispectral camera and odometry, and an actuation system based on a robotic arm and a pumping system for liquid fertiliser.The system gathers information from the plants in the row about crop volume (lidar), crop health (multispectral images) and plant position (point cloud for navigation).Then a decision about fertilisation need (yes/no) is made and the actuation system applies liquid fertiliser to the soil around the plant avoiding contacting the crop.The prototype can be enhanced with more sensors, and it can be also adapted to perform other types of field work.

Figure 3 .
Figure 3. Robotic fertilization prototype in the cultivation row.Red: the tank with liquid treatment, yellow: the nozzle.

Figure 4
Figure 4 Process of calculation of the vegetation index from spectral data, and its correspondence with histogram information.

Figure 7 K
Figure 7 K-means classification and geometrical parameters extraction from a row crop.

Figure 8
Figure 8 Matching key points from the L-PC (upper-two captured partial sections) to the G-PC (lower).

Figure 9 Figure 9
Figure 9 shows the perception system (G-PC and robot visualisation), where the point cloud and the estimated position of the platform are shown for each instance.As the robotic platform advances from point 1 through 7, the L-PC accumulates earlier values, which continuously improves the location estimation.

Figure 10
Figure 10 Box diagram of mean localisation errors at each of the positions indicated in Figure 9.