System and method of carbon emission mapping

Abstract

A system and method for carbon emission mapping. The method comprises obtaining a plurality of first interferograms based on a plurality of radar images of a peatland; estimating an ionospheric phase contribution using an ionospheric correction on the plurality of first interferograms; using the ionospheric phase contribution, applying the ionospheric correction to a plurality of second interferograms; performing Small BAseline Subset (SBAS) network inversion on the plurality of second interferograms to obtain time series data; applying a closure phase bias correction to the time series data; applying a tropospheric noise correction to the time series data to obtain ground displacement data; and generating a carbon emission map of the peatland based on the ground displacement data.

Description

SYSTEM AND METHOD OF CARBON EMISSION MAPPING CROSS-REFERENCE TO RELATED APPLICATION

[0001] This application claims the benefit of priority to the Singapore application no.10202500656X filed 14 March, 2025, the contents of which are hereby incorporated by reference in their entirety for all purposes.TECHNICAL FIELD

[0002] This application relates generally to the field of carbon emission mapping, and more particularly, to a system for carbon emission mapping and a method of carbon emission mapping using radar images of a peatland.BACKGROUND

[0003] Peatlands, despite covering only 3-4% of global land area, store around 33% of global soil carbon, making them vital in the carbon budget. However, degradation of peatlands can turn them into significant carbon emitters. Peatland restoration and conservation are thus important for meeting carbon emission targets. However, the uptake of restoration and conservation projects in voluntary carbon markets has been impeded, largely due to the high cost and scarcity of robust carbon emission measurements, as field instruments such as gas flux towers and gas chambers are cumbersome and expensive to deploy at scale.SUMMARY

[0004] According to an aspect, disclosed herein is a method of carbon emission mapping. The method comprises obtaining a plurality of first interferograms based on a plurality of radar images of a peatland, wherein the plurality of radar images corresponds to radar images of thepeatland captured over a plurality of time instances; estimating an ionospheric phase contribution using an ionospheric correction on the plurality of first interferograms, wherein the plurality of first interferograms corresponds to a first bandwidth configuration; using the ionospheric phase contribution, applying the ionospheric correction to a plurality of second interferograms, wherein the plurality of second interferograms corresponds to a second bandwidth configuration; performing Small BAseline Subset (SBAS) network inversion on the plurality of second interferograms to obtain time series data; applying a closure phase bias correction to the time series data; applying a tropospheric noise correction to the time series data to obtain ground displacement data; and generating a carbon emission map of the peatland based on the ground displacement data.

[0005] According to another aspect, disclosed herein is a system. The system comprises: memory storing instructions; and a processor coupled to the memory and configured to process the stored instructions to implement: a module configured to perform the method of carbon emission mapping.BRIEF DESCRIPTION OF THE DRAWINGS

[0006] Various embodiments of the present disclosure are described below with reference to the following drawings:FIG. 1 is a system diagram for a system for carbon emission mapping according to embodiments of the present disclosure.FIG. 2 is a flowchart illustrating a method of carbon emission mapping according to various embodiments.FIG. 3A illustrates an exemplary interferogram.FIG. 3B illustrates the definition of the term bandwidth.FIG. 4 is a flowchart illustrating a method of ionospheric correction according to various embodiments.FIG. 5 illustrates an exemplary outcome of ionospheric correction, with the masking of decorrelated areas improving the estimated ionospheric phase and eliminating artefacts over decorrelated areas in the corrected interferogram - an example for an interferogram between 14 February 2020 and 13 March 2020, spanning the entire ALOS-2 ScanSAR image of 350 x 350 km.FIG. 6 is a flowchart illustrating a method of closure phase bias correction according to various embodiments.FIGs. 7A and 7B show a comparison of ground measured and InSAR peat motions using different number of interferogram bandwidths to reduce closure phase bias.FIG. 8 is a flowchart illustrating a method of tropospheric noise correction according to various embodiments.FIGs. 9A and 9B show a comparison of ground measured and InSAR peat motions with and without the spatially correlated tropospheric noise (SCTN) correction.FIGs. 10A and 10B show field-measured peat motion in drained peatlands of Palangkaraya, Central Kalimantan, Indonesia according to an exemplary implementation. FIG. 10A shows land use, drainage canals, peat camera sites, and InSAR coverage (inset). FIG. 10B shows vertical displacement measured by peat cameras and 30-day mean rainfall recorded at the nearby Tjilik Riwut Airport weather station between 2018 and 2024. Positive and negative values indicate uplift and subsidence respectively throughout. Dry periods with <5 mm 30-day mean rainfall, typically between May and September, are highlighted.FIGs. 11 A and 1 IB show field measurements of peat motion and groundwater level change. Peat motion and groundwater level changes measured by peat cameras at 8 sites in Palangkaraya, Central Kalimantan, Indonesia show high correlation. Dry periods with <5 mm 30-day mean rainfall, typically between May and September, are highlighted.FIGs. 12A and 12B show a comparison of InSAR-derived peat motions against peat camera field measurements. InSAR captured both long-term subsidence and episodic fluctuations in peat motion at several field validation sites. InSAR displacements were projected from the line-of-sight to vertical direction. Error bars indicate the InSAR phase standard deviation, <T which was derived from coherence that was propagated through the inversion matrix used for time series analysis. The peat camera at Misik plantation likely detected stronger, localized subsidence compared to its surroundings during the 2019 El Nino - a downshifted InSAR time series shows better fit with peat camera displacements after the 2019 El Nino.FIGs. 13A and 13B show local-scale heterogeneity of peat motions within the proximity of InSAR pixel centres. FIG. 13A shows peat motion measured (1) by the peat camera, (2) manually at the subsidence pole at the peat camera, (3) by InSAR at the InSAR pixel centre closest to the peat camera, and (4) manually at nearby subsidence poles. Note that manual measurements made at the subsidence poles are subject to human errors and inconsistencies in method of measurement over time. FIG. 13B shows locations of the various measurement types. At Misik plantation, subsidence at the surrounding subsidence poles were all less than at the peat camera after the 2019 El Nino. The InSAR measurements are more similar to the surrounding subsidence poles than the peat camera, possibly as InSAR reflects the landscapeaverage within the pixel rather than localised subsidence at the peat camera. At Scrubland south, peat motions at the surrounding subsidence poles vary from the peat camera. The InSAR measurements lie within the range of subsidence measured in the field, also possibly as InSAR reflects the landscape-average rather than specific.FIGs. 14A and 14B show ionospheric phases estimated for each image acquisition. The phase is assumed to be zero for the first acquisition. White polygons represent the location of Palangkaraya city. FIG. 14A shows wrapped ionospheric phases with spatial complexities that are not strictly linear or quadratic. FIG. 14B shows unwrapped ionospheric phases, where thelargest offset caused by ionospheric noise between any two locations is up to 6.7 m and averaging 2.3 ± 1.5 m across the entire ALOS-2 ScanSAR image of 350 km x 350 km.FIGs. 15A to 15D show the impact of processing parameters on InSAR-derived peat motions. FIG. 15A shows the displacement time series at 8 sites using different interferogram bandwidths (BW), with ionospheric, tropospheric, and MintPy closure phase bias corrections applied. The MintPy correction is not applicable for bandwidth- 1 as it requires a mix of temporal baselines. Bandwidth- 1 results thus have no closure phase bias mitigated, neither through the MintPy correction nor the use of both short and long temporal baselines; while bandwidth-20 results have the most bias mitigated through the MintPy correction and the use of largest bandwidth. FIGs. 15B to 15D show the RMSE, r, and long-term rate difference between peat cameras and InSAR depending on the interferogram bandwidth, and whether closure phase bias (clo. ph.) and tropospheric (tropo.) corrections (corr.) were applied, and where ionospheric corrections were applied in all variations. All InSAR displacements have been projected from the line-of-sight to vertical direction.FIGs. 16A to 16E show large-scale trends in InSAR-derived peat motions and their correlation with different environmental factors. FIGs. 16A and 16B show areas that subsided more during the 2019 El Nino did not necessarily subside more in the long term. Long-term rates showed weak correlation with land use, proximity from the peat edge, and drainage canal density. FIG.16C shows land use, where centre lines refer to the mean. FIG. 16D shows proximity from the peat edge. FIG. 16E shows drainage canal density. All InSAR motions were projected from the line-of-sight to vertical direction and are relative to a stable reference point in the city centre of Palangkaraya. Areas with high mean InSAR phase standard deviation, <T > 0.7 cm were excluded from the analysis.FIGs. 17A to 17F show the relationship between long-term rate and maximum subsidence by land uses. Areas that subsided more during the 2019 El Nino dry season, based on the maximumsubsidence detected between 15 May 2019 and 16 November 2019, do not necessarily subside more in the long-term - a comparison of pixels within each land use.FIG. 18 shows spatial and temporal baselines of ALOS-2 ScanSAR image acquisitions used. The triangle represents the reference image acquired on 19 June 2020, dots represent the 91 other images, and lines represent interferograms formed between image pairs that were used in the interferogram network inversion.FIG. 19 shows InSAR phase standard deviations, <T by land uses. Forests have higher mean <T compared to other land uses. However, the mean <T still varies greatly within each land use. Numbers in brackets indicate the number of pixels present for each land use within the analysed peat area in Central Kalimantan, Indonesia.FIG. 20 is a schematic diagram of a processor system.DETAILED DESCRIPTION

[0007] The following detailed description is made with reference to the accompanying drawings, showing details and embodiments of the present disclosure for the purposes of illustration. Features that are described in the context of an embodiment may correspondingly be applicable to the same or similar features in the other embodiments, even if not explicitly described in these other embodiments. Additions and / or combinations and / or alternatives as described for a feature in the context of an embodiment may correspondingly be applicable to the same or similar feature in the other embodiments.

[0008] In the context of various embodiments, the articles “a”, “an” and “the” as used with regard to a feature or element include a reference to one or more of the features or elements.

[0009] In the context of various embodiments, the term “about” or “approximately” as applied to a numeric value encompasses the exact value and a reasonable variance as generally understood in the relevant technical field, e.g., within 10% of the specified value.

[0010] As used herein, the term “and / or” includes any and all combinations of one or more of the associated listed items.

[0011] As used herein, terms “concurrently”, “simultaneously”, “at the same time”, or the like, may refer to events or actions that coincide or overlap within a period of time, regardless of whether the events start at the same time instant, and regardless of whether the events end at the same time instant.

[0012] As used herein, the term “interferogram” refers to a data or a data product derived from a phase comparison of two microwave radar images acquired by a Synthetic Aperture Radar (SAR) satellite. In Interferometric Synthetic Aperture Radar (InSAR), the interferogram encodes phase differences resulting from variations in the radar signal path length between image acquisitions, which may be caused by ground surface motion, atmospheric effects, topography, orbital variations, noise, or combinations of thereof. The interferogram may be processed to quantify ground surface motion, including subsidence or uplift in peatlands.

[0013] Conventional approaches towards carbon emission measurements on peatlands often involve the use of field instruments such as gas flux towers and gas chambers, both of which are cumbersome and expensive to deploy at scale. The application of InSAR is one approach in measuring peat surface motion, which may act as a high-resolution proxy for peat carbon flux. However, such approaches are often technically challenging, and especially so in the tropics due to denser vegetation over peat, in comparison to temperate and boreal regions.

[0014] The present disclosure discloses an InS AR-based method for measuring tropical peat surface motion as a high-resolution proxy for peat carbon flux, utilizing L-band SAR data. The proposed method utilizes a unique combination of InSAR processing techniques to achieve higher measurement accuracies than demonstrated by existing InSAR approaches in tropical peatlands, specifically for vegetated areas including plantations, cropland, and scrubland. Since the oxidation of peat and release of gaseous or dissolved carbon is accompanied by the loweringof the peat surface, the peat motions are proxies that can be converted to the peat carbon flux. In avoidance of doubt, the proposed method may also be used in temperate and boreal regions.

[0015] The proposed method utilizes peat motions as proxies for the measurement of peat carbon flux, and may comprise various advantages, such as:

[0016] 1. Cost-effective: The proposed method enables remote monitoring over large areas with continuous spatial coverage and significantly reduces dependency on existing methods involving labour-intensive measurements on the ground.

[0017] 2. Continuous monitoring: The proposed method enables frequent and regular remote monitoring of target areas, independent of weather conditions and time of day.

[0018] 3. Enhanced credibility: The proposed method may be rigorously validated against direct on-the-ground measurements and boost credibility, which is essential in applications such as carbon credit assessment. In contrast, conventional methods using InSAR for tropical peatlands lacked proper validation or failed to directly correlate with actual peat movement.

[0019] 4. Improved accuracy of long-term tropical peat motion measurement: The proposed method may achieve accuracies below 1 cm / yr, as opposed to unknown accuracies or estimated accuracies of several cm / yr in conventional methods using InSAR for tropical peatlands. Subcentimetre accuracy is significant as a 1 cm / yr peat motion translates to substantial carbon emissions of around 16.15 tonnes CCh / ha / yr.

[0020] 5. Intra-annual peat motion measurement capability: The proposed method comprises the ability to measure shorter-term, seasonal fluctuations in peat motion, going beyond the traditional approach of reporting only long-term, annual averages. This enables monitoring of dynamic changes, capturing seasonal variability and shifts in peat management and restoration. This deeper insight supports more adaptive strategies for peatland management and carbon projects.

[0021] Carbon emission estimates from direct measurements or approximations based on peat motion remain sparse in Southeast Asian tropical peatlands. There has been no robust validation of L-band InSAR peat motion measurements in the region. In view of the above challenges, disclosed herein is an InSAR processing method and / or approach in carbon emission mapping. The proposed method provides carbon emission estimations of high spatial and temporal resolution and the capability to scale up measurements wherever the L-band InSAR images are available. This offers project developers in the voluntary carbon market a cost-effective and reliable tool for monitoring, reporting, and verifying peat carbon at a large scale, reducing their reliance on sparse, heavily interpolated data, and boosting the credibility of carbon credits and their appeal to investors.

[0022] According to an aspect of the disclosure, referring to FIG. 1, disclosed herein is a carbon emission mapping system 100. The system 100 may be an InSAR processing-based system for carbon emission mapping of a peatland. The system 100 may be configured to perform a method of carbon emission mapping by determining peat carbon flux based on a plurality of radar images of the peatland. According to various embodiments, the system 100 may determine peat motions derived from L-band InSAR, and thereafter use the peat motions to approximate peat carbon flux by multiplying the peat motion, carbon concentration, and dry bulk density of the peat.

[0023] According to various embodiments, the proposed system 100 involves a combination of modules, with the operation of these modules using specific parameters to achieve accuracy in peat motions and fine-tuned to increase computational efficiency to enable scaling-up of the processing.

[0024] In various embodiments, the system 100 may comprise an interferogram module 200. The interferogram module 200 may determine 210 (or obtain) a plurality of interferograms based on a plurality of radar images 82 of a peatland. The interferogram module 200 maydetermine a multitude of plurality of interferograms, such as a plurality of first interferograms with (or corresponding to) a first bandwidth configuration, and a plurality of second interferograms with (or corresponding to) a second bandwidth configuration. In various embodiments, the second bandwidth configuration may be different from the first bandwidth configuration. In various embodiments, the second bandwidth configuration may be larger than the first bandwidth configuration.

[0025] The bandwidth may be defined as the maximum number of subsequently acquired radar images that may be paired with each radar image to form a network of interferograms. FIG. 3B illustrates the definition of bandwidth. Referring to FIG. 3B, if there are 50 radar images in total, a bandwidth-3 interferogram network may comprise unique pairs of radar images including: 1-2, 1-3, 1-4; 2-3, 2-4, 2-5; ...; 47-48, 47-49, 47-50; 48-49, 48-50; 49-50. It may be noted that each image is paired with 3 other images acquired after it. For the last few images (e.g. image 48 which forms 48-49 and 48-50 unique pairs), if there are less than 3 images acquired after it, there may be less than 3 unique pairs.

[0026] In various embodiments, the plurality of radar images 82 may correspond to radar images of the peatland captured over a plurality of time instances. In other words, the system 100 may perform the carbon emission mapping based on radar images 82 captured over a time duration or captured over multiple time instances.

[0027] The plurality of radar images 82 of the peatland may be obtained from one or more radar modules 80. Alternatively, the plurality of radar images 82 of the peatland may be obtained from a database comprising multiple radar images.

[0028] According to various embodiments, the system 100 may further comprise: an ionospheric correction module 300 for performing a method of ionospheric correction 310; a Small BAseline Subset (SBAS) network inversion module 400 for performing a method of network inversion 410; a closure phase bias correction module 500 for performing a method ofclosure phase bias correction 510; a tropospheric noise correction module 600 for performing a method of tropospheric noise correction 610; and a carbon emission mapping module 700 for performing a method of generating a carbon emission map 710.

[0029] Referring to FIG. 2, according to various embodiments, the system 100 may sequentially apply each of the method of ionospheric correction 310, the method of network inversion 410; the method of closure phase bias correction 510, the method of tropospheric noise correction 610; and the method of generating the carbon emission map 710 in a series or cascading workflow. Alternatively, the system 100 may selectively apply one or more of: the method of ionospheric correction 310, the method of closure phase bias correction 510, and the method of tropospheric noise correction 610 based on the artefacts / errors as present in the actual interferograms as obtained. FIG. 3A illustrates an exemplary interferogram comprising artefacts.

[0030] In various embodiments, the ionospheric correction may be performed at the interferogram level (i.e. for interferograms networks). Thereafter, a network inversion (such as SBAS network inversion) operates on the interferogram network to yield a time series. Thereafter, the closure phase bias correction and tropospheric noise correction may be applied to the time series to obtain the ground displacement. Carbon emissions may be computed upon stabilization of ground displacement.

[0031] Ionospheric correction module 300 and method 310

[0032] FIG. 4 illustrates a method of ionospheric correction 310 performed by the ionospheric correction module 300. The method of ionospheric correction 310 may be performed using the InSAR Scientific Computing Environment (ISCE) software, which performs range split-spectrum ionospheric correction.

[0033] Referring to FIG 4., in various embodiments, the ionospheric correction module 300 may be configured to improve the range split-spectrum ionospheric correction results. Aplurality of first interferograms may be formed from a plurality of radar images using a first bandwidth configuration. The plurality of first interferograms may be multi-looked to improve signal-to-noise ratio. In some embodiments, decorrelated pixels may be masked prior to ionospheric phase estimation, and spatial interpolation may be applied to reconstruct continuous ionospheric phases after ionospheric phase estimation. Ionospheric phase contributions associated with each radar image may be determined through network inversion. A plurality of second interferograms may then be formed from the same radar images using a second, larger bandwidth configuration. In some embodiments, the second interferograms may be multi-looked using fewer looks than the first interferograms, and the ionospheric phase contributions may be applied to the second interferograms.

[0034] In various embodiments, the method of ionospheric correction 310 may comprise: in 320, forming a plurality of first interferograms from a plurality of radar images using a first bandwidth configuration. In some embodiments, the first bandwidth configuration may correspond to a bandwidth-4 interferogram network.

[0035] In various embodiments, the method of ionospheric correction 310 may comprise: in 330, multi-looking each of the plurality of first interferograms to improve signal-to-noise ratio while maintaining sufficient spatial resolution to estimate long-wavelength ionospheric phase variations. For example, multi-looking may be performed using approximately 80 looks in range and 448 looks in azimuth, although other multi-look factors may also be used.

[0036] In various embodiments, the method of ionospheric correction 310 may comprise: in 340, identifying and removing decorrelated pixels in each of the plurality of first interferograms prior to ionospheric phase estimation. Decorrelated areas may occur over rivers or other water bodies that may not be captured in external datasets. Removal of decorrelated pixels may comprise applying a respective mask to the decorrelated areas of the interferograms. The mask may be generated using a mapping algorithm which converts delineated pixels on aninterferogram into bounding box coordinates compatible with processing performed by the ISCE alosStack module. An example of an improved ionospheric correction due to the removal of decorrelated areas is shown in FIG. 5.

[0037] In various embodiments, the method of ionospheric correction 310 may comprise: in 350, ionospheric phase estimation may be performed for each of the plurality of first interferograms using a split-spectrum ionospheric correction approach to obtain ionospheric phase estimates. Each interferogram may be separated into two sub-bands, and a respective ionospheric phase delay may be estimated based on frequency-dependent phase differences between the sub-bands.

[0038] In various embodiments, the method of ionospheric correction 310 may comprise: in 360, applying spatial interpolation to fdl gaps corresponding to the removed decorrelated pixels, thereby producing spatially continuous ionospheric phase estimates for each of the plurality of first interferograms.

[0039] In various embodiments, the method of ionospheric correction 310 may comprise: in 370, based on the ionospheric phase estimates of the plurality of first interferograms, obtaining an ionospheric phase contribution associated with each radar image using network inversion.

[0040] In various embodiments, the method of ionospheric correction 310 may comprise: in 380, a plurality of second interferograms may be formed from the plurality of radar images using a second bandwidth configuration that is larger than the first bandwidth configuration. In some embodiments, the second bandwidth configuration may correspond to a bandwidth-20 interferogram network. The ionospheric phase contribution associated with each radar image may be paired according to the second bandwidth configuration to derive ionospheric phase corrections for the plurality of second interferograms, and the ionospheric phase corrections may be applied to the second interferograms.

[0041] In various embodiments, the method of ionospheric correction 310 may comprise: in 390, each of the plurality of second interferograms may be multi-looked to improve signal-to-noise ratio. The number of multi-looks used for the plurality of second interferograms may be smaller than the number of multi-looks used for the plurality of first interferograms. For example, multi-looking may be performed using approximately 17 looks in range and 56 looks in azimuth for the plurality of second interferograms, although other multi-look factors may also be used.

[0042] According to various embodiments, the method of ionospheric correction 310 may comprise: estimating an ionospheric phase contribution using an ionospheric correction on the plurality of first interferograms; and using the ionospheric phase contribution, applying the ionospheric correction to a plurality of second interferograms. The method of ionospheric correction 310 may further comprise: removing decorrelated pixels in each of the plurality of first interferograms; and interpolating gaps corresponding to the decorrelated pixels in each of the plurality of first interferograms, using the ionospheric phase estimation. In various embodiments, removing decorrelated pixels in each of the plurality of first interferograms comprises: applying a respective mask to the decorrelated pixels in each of the plurality of first interferograms. The method of ionospheric correction 310 may further comprise: applying a split-spectrum ionospheric correction to each of the plurality of first interferograms, wherein each of the plurality of first interferograms is split into two sub-bands, and wherein a respective ionospheric phase delay is estimated based on frequency-dependent phase differences between the two sub-bands.

[0043] Network inversion module 400 and method 410

[0044] In various embodiments, the Small BAseline Subset (SBAS) network inversion module 400 may perform a method of network inversion 410. The method of network inversion 410 may comprise: performing Small BAseline Subset (SBAS) network inversion on theplurality of second interferograms to obtain a time series data. In various embodiments, the plurality of second interferograms are weighted.

[0045] Closure phase bias correction module 500 and method 510

[0046] FIG. 6 illustrates a method of closure phase bias correction 510 performed by the closure phase bias correction module 500. The method of closure phase bias correction 510 may be implemented in the ISCE software which creates a stack of interferograms, operated with specification of the ideal interferogram network configuration to reduce the effect of the closure phase bias.

[0047] In exemplary implementations, existing module in the Miami InSAR time-series software in Python (MintPy) software may be used to compute the amount of closure phase bias, operated with specification of the ideal number of connections and assumptions which may be determined based on characteristics of the InSAR pixel. The closure phase bias correction module 500 may be automated in the determination of the number of connections and assumptions for each InSAR pixel and apply a pixel-by-pixel closure phase bias correction as opposed to a uniform correction.

[0048] Closure phase bias refers to the discrepancy amongst triplets of phases (c|>), where <t>AC ~ 4>AB ~ 4>BC 0 and A, B, and C are different time epochs. To minimize the effect of closure phase biases, several parameters may be used: (i) a redundant number of interferograms from a bandwidth-20 network, weighted by nJ; for the SBAS network inversion, and (ii) a correction method in MintPy.

[0049] An infinitely large bandwidth may be the most ideal to minimize the effects of closure phase bias as demonstrated by several studies based on C-band Sentinel- 1 images. However, this is not practical for L-band ALOS-2 images which may decorrelate faster with increasing bandwidth due to the longer temporal baselines, commonly of 14 or 28 days, as compared to 6 or 12 days for Sentinel- 1. Large bandwidths are also not practical for most areasin the tropical peatlands where coherence is already low and large bandwidth interferograms would be mostly decorrelated. The inclusion of decorrelated interferograms will then have insignificant impact on the SBAS network inversion since the interferograms are weighted by

[0050] In view of the above considerations, it was found that by increasing the number of bandwidths used up to 20, convergence of the SBAS network inversion may be observed (see FIGs. 7A and 7B). This is computationally efficient and maximizes the performance of the closure phase bias correction in MintPy. As shown in FIGs. 7 A and 7B, the closure phase bias is continually reduced with increasing number of bandwidths used, especially for the scrubland and forest sites. Although the correction method had minimal impact on the accuracy when bandwidth-20 was used, the correction may remain beneficial especially if small bandwidths are used (FIGs. 7A and 7B). The correction may also be beneficial or even critical over scrubland and forests which experience the largest bias, followed by cropland, then plantations (FIGs. 7A and 7B). Finally, the quality of the estimated closure phase depends on the coherence and cannot be reliable estimated if the data is decorrelated. Hence, different levels of closure phase bias correction may be applied to each pixel to avoid introducing unwrapping errors during the correction - using smaller bandwidths to estimate the bias for pixels which decorrelate faster, and larger bandwidths for those which decorrelate slower. To determine the appropriate level of correction for each pixel, an algorithm may be applied to identify instances of phase unwrapping errors at different bandwidths, which indicate decorrelated data. The algorithm may then select, for each pixel, the maximum bandwidth at which no unwrapping errors are observed, thereby defining the appropriate level of correction. In some embodiments, smaller bandwidths may be used for pixels that decorrelate quickly, and larger bandwidths for pixels that remain coherent over longer temporal separations.

[0051] Referring to FIG. 6, in various embodiments, the method of closure phase bias correction 510, may comprise: in 520, applying a closure phase bias correction to the time series data. Applying the closure phase bias correction to the time series data in 520, may comprise: in 530, determining a plurality of closure phase biases based on the time series data; and in 540, for each pixel in the time series data, applying a pixel-by-pixel closure phase bias correction using a selected one of the plurality of closure phase biases to obtain the corrected time series data. In other words, different levels of closure phase biases may be determined for each pixel in the time series data.

[0052] In various embodiments, applying the closure phase bias correction to the time series data in 520, may further comprise: in 550, determining instances of unwrapping error for the corresponding pixel in the time series data; and in 560, determining the selected one of the plurality of closure phase biases corresponding to no unwrapping error.

[0053] In various embodiments, applying the closure phase bias correction to the time series data in 520, may further comprise: in 570, determining the selected one of the plurality of closure phase biases corresponding to a largest temporal separation. In other words, the closure phase biases for each pixel may be determined based on the largest temporal separation where no unwrapping error is present.

[0054] Tropospheric noise correction module 600 and method 610

[0055] FIG. 8 illustrates a tropospheric noise correction 610 performed by the tropospheric noise correction module 600. The tropospheric noise correction module 600 may be implemented to reduce anomalies in the time series data using the self-developed algorithm. The tropospheric noise correction method may be a spatially correlated tropospheric noise correction method.

[0056] The short-wavelength tropospheric noise in the InSAR time series may be corrected using an empirical approach. Tropospheric noise may have magnitudes up to tens ofcentimetres. It may be spatially correlated over varying distances of few to tens of kilometres, but temporally uncorrelated over a few hours or less. Shorter-wavelength noise of few kilometres is mostly attributed to the wet tropospheric component that is driven by turbulent water vapor changes, while longer-wavelength noise of tens of kilometres is attributed to the hydrostatic component that varies more smoothly with pressure, temperature, and elevation-related water vapor content.

[0057] Conventional techniques based on Global Atmospheric Models (GAMs) and height correlation models are generally not effective for estimating these tropospheric effects, likely because of the low resolution of the atmospheric variables in GAMs which are not suited for removing turbulent tropospheric noise. In tropical peatlands which are degraded and have mostly flat terrains, the height correlation-based method is not meaningful as it operates as a function of the thickness of the troposphere that the InSAR microwave propagates through.

[0058] Instead, the short-wavelength tropospheric noise may be reduced at a pixel scale, i.e. with the short-wavelength tropospheric noise of each pixel reduced utilizing the stacked time series of neighbouring non-deforming pixels. On a basis that tropospheric noise is spatially correlated, the signals at non-deforming pixels )non—def close to each other may contain predominantly tropospheric noise that are similar across pixels, and to a smaller extent other random noise that may be dissimilar across pixels:< <g non-def > . non-def , , non-defV—'rtropo "T" 'rnoise GG^targetrefers(0the phase of a target pixel in which the tropospheric noise is to be removed. Therefore, stacking the time series of these neighbouring non-deforming pixels, such as by taking the mean, may amplify the spatially correlated tropospheric noise while cancelling outother random noise. This may result in an approximate of the spatially correlated tropospheric noise around a target pixel that can then be removed from the target pixel itself:>

[0059] The noise correction may be performed for each pixel independently. In an example, the non-deforming neighbouring pixels may be defined as those with a velocity of less than 1 mm / year over the whole time series, an average spatial coherence above 0.3, and a maximum distance of 10 km away from the target pixel. The maximum distance may be selected to be within typical spatial correlation lengths of wet tropospheric noise of not more than tens of kilometres. Finally, only target pixels with at least 5 non-deforming neighbouring pixels are corrected for the spatially correlated tropospheric noise to avoid the inclusion of other random noise. This correction reduces the noise in the time series as shown in FIGs. 9A and 9B, where there are fewer number of anomalous values for several different test sites.

[0060] Referring to FIG. 8, in various embodiments, the method of tropospheric noise correction 610, may comprise: in 620: applying a tropospheric noise correction to the time series data to obtain ground displacement data. In various embodiments, applying the tropospheric noise correction to the time series data 620, comprises: in 630, for each target pixel in each of the time series data, correcting the target pixel with a respective spatially correlated tropospheric noise to obtain the ground displacement data.

[0061] In various embodiments, applying the tropospheric noise correction to the time series data 620, further comprises: in 640, for each target pixel in each of the time series data, determining the respective spatially correlated tropospheric noise based on a plurality of neighbour pixels around the target pixel. In other words, the system 100 may determine spatially correlated tropospheric noise based on neighbouring pixels.

[0062] In various embodiments, applying the tropospheric noise correction to the time series data 620, further comprises: in 650, for each target pixel in each of the time series data, determining the respective spatially correlated tropospheric noise around the target pixel based on a mean of stacked time series of the plurality of neighbour pixels of the target pixel.

[0063] In various embodiments, applying the tropospheric noise correction to the time series data 620, further comprises: in 660, determining at least one target pixel from the time series data, wherein the at least one target pixel comprises a plurality of non-deforming neighbour pixels, wherein a count of the plurality of non-deforming neighbour pixels is equal to or above a minimum neighbour pixel threshold. As such, the tropospheric noise correction 610 may be performed based on a minimum number of non-deforming neighbour pixels. In exemplary embodiments, the minimum neighbour pixel threshold is 5.

[0064] In various embodiments, each of the plurality of non-deforming neighbour pixels is characterized by: a velocity under a velocity threshold; an average spatial coherence above a spatial coherence threshold; and a distance away from the target pixel under a maximum distance threshold.

[0065] Carbon emission mapping module and method

[0066] In various embodiments, the carbon emission mapping module 700 may perform a method of generating a carbon emission map 710 of the peatland based on the ground displacement data. The carbon emission map 710 may be used to estimate or determine the respective carbon emission over a duration.

[0067] Exemplary implementation of the proposed method and system

[0068] This section describes an exemplary implementation of the proposed system and method, wherein an L-band InSAR system was demonstrated for mapping tropical peat motion The exemplary implementation was validated both spatially and temporally against groundtruth peat motion data. Using an InSAR framework with improved noise corrections and high-rate field measurements, the implementation demonstrated a route towards cost-effective, large-scale peat motion monitoring as a proxy for tropical peatland carbon balance, which is an essential prerequisite for addressing this globally important source of carbon emissions.

[0069] Peat motion from field measurements

[0070] Tropical peatlands of Palangkaraya, Central Kalimantan, Indonesia, was selected for the validation of InSAR-derived peat motion, where high-resolution field measurements at 8 sites were collected (FIG. 10A). All sites are within or adjacent to the former Mega Rice Project area, where over 1 million hectares of peat swamp forest were deforested and over 4,000 km of canals were constructed to drain peat for unsuccessful rice cultivation between 1996 and 1998. Excessive drainage and subsequent partial land abandonment led to dry, scrub-covered peat that is vulnerable to frequent fires, especially during dry seasons, notably causing major fires in 1997 and 2015. The peatlands mainly comprised smallholder cropland, smallholder and industrial plantations, scrubland, and secondary forests. Drainage continued through the formerly constructed network of canals and smaller canals subsequently dug (FIG. 10A) to support agricultural activities. This had resulted in rapid, ongoing subsidence driven by compaction, consolidation, and oxidation of the peat, with oxidation also leading to high CO2 emissions.

[0071] The field measurements showed both high levels of episodic fluctuation in peat surface elevation, linked to dry-wet cycles, and long-term subsidence at drained sites (FIG.10B). The 8 sites comprised different land uses including: oil palm plantations, cropland, fern-dominated scrubland, and secondary forests. Measurements were made using peat cameras between 1 November 2018 and 24 March 2024. These time-lapse cameras moved with the peat surface and captured images of a stable reference point, providing sub -millimetre precision of peat motion at a two-hourly frequency. Short periods of data loss affected all sites in late 2021 and early 2022, however, the datum point remained unchanged between measurement periods.Episodic subsidence was observed at most sites during the dry seasons from 2019 to 2023, typically occurring between May and September when minimal rainfall induces lowered groundwater levels (FIGs. 11 A and 1 IB). The dry seasons, with <5 mm 30-day mean rainfall recorded at the nearby Tjilik Riwut Airport weather station (2.22°S, 113.95°E, 10 to 65 km away from sites), are highlighted in FIG. 10B. UPR forest did not subside in most years as it was mostly flooded most of the time based on groundwater levels (FIGs. 11A and 11B) and field observations. The largest episodic subsidence, up to 8 cm, occurred in 2019 and 2023 during dry seasons. These subsidence events coincided with positive phases of the National Oceanic and Atmospheric Administration Oceanic Nino Index and were thus likely intensified by the interannual El Nino, which causes drought-like conditions every few years in the region. Both events were followed by several centimetres of uplift, coinciding with the onset of rain. However, rebounded peat surface elevations were lower than before the dry seasons despite similar groundwater levels, suggesting permanent peat loss due to oxidation, or irreversible compaction and consolidation due to drying.

[0072] Validation of InSAR peat motions

[0073] Using ALOS-2 ScanSAR imageries acquired between 12 October 2018 and9 August 2024, displacement time series were derived. The displacement time series showed peat motion at each acquisition date and its evolution over time, and estimated the corresponding long-term rates, which represent the average displacement per year over the full observation period. ALOS-2 is a L-band radar satellite capable of providing repeat observations every 14 days in ScanSAR mode, with 350 km x 350 km coverage per image. Using an InSAR time series analysis approach, spatially continuous measurements at 180 m x 180 m spatial resolution across Central Kalimantan were produced (FIG. 10A).

[0074] The L-band InSAR results captured episodic fluctuations observed in field measurements at majority of the sites, including plantations, cropland, and scrubland. Largesubsidence and rebound during and after both El Nino dry seasons, and greater stability during non-El Nino periods were similarly observed in InSAR (Error! Reference source not found.).The average root mean squared error (RMSE) and Pearson’s correlation coefficient (r) between epoch-to-epoch displacements of the peat cameras and InSAR were 1.70 + 0.31 cm and 0.75 + 0.09 across these 6 sites, and 1.75 + 0.34 cm and 0.70 + 0.14 across all sites, including forest sites (Table la and b). The lowest RMSE was 1.27 cm at Hampangen plantation, and the highest r was 0.88 at Misik plantation (Table la and b). The sites where InSAR performed less optimally were Kalampangan cropland (high RMSE of 2.21 cm, but high r of 0.83) and UPR forest (low r of 0.40, but low RMSE of 1.55 cm). In both cases, InSAR captured the long-term trend, showing subsidence at the cropland site and no change at the forest site, but overestimated the amount of short-term motion compared to observed values (Error! Reference source not found.). At the second forest site at KHDTK, a long data gap in the peat camera data precluded a full assessment of short-term variations.

[0075] The L-band InSAR results closely matched long-term, multiannual rates of peat motion measured in the field. Both the peat cameras and InSAR showed rates around -1 to 0 cm yr'1across all sites over 5.5 years (Table 1c and d, negative values indicate subsidence throughout). The smallest difference between peat camera and InSAR rates was 0.06 cm yr'1at UPR forest, and the largest was -0.59 cm yr'1at RePeat scrubland (Table le). The RMSE of the rates across the 8 sites was 0.36 cm yr'1, with no overall tendency of InSAR to underestimate or overestimate the rates.

[0076] No previous tropical peatland study had directly compared InSAR with in-situ field-measured peat motion. For example, previous studies only compared overall InSAR and field-measured peat motion ranges without site-to-site validation, or compared InSAR indirectly against peat motion inferred from in-situ groundwater levels. The results therefore demonstrateaccuracy and represent the most direct one-to-one field validation of InSAR-derived peat motion in the region thus far.Table 1: InSAR accuracy against peat camera field measurementsa, bRMSE and r were computed between epoch-to-epoch displacements of InSAR and peat cameras.c, , e £)ifferencesinthe long-term, multiannual rates were computed from the peat camera velocity minus the InSAR velocity.fLower mean InSAR phase standard deviation,indicates better ease of tracking motion for an InSAR pixel.InSAR motions were projected from the line-of-sight to vertical direction for all calculations.

[0077] Deviations between InSAR and field measurements are expected given their differing spatial and temporal resolutions. InSAR detects average peat motions within 180 m x 180 m pixels, while field measurements are point -based. Deviations thus arise from local-scale heterogeneity, especially in heavily managed areas where varying management or proximity to ditches caused varying subsidence within an image pixel. At Misik plantation, nearby manual subsidence pole measurements showed less subsidence than at the peat camera post-2019 El Nino and align better with InSAR (Error! Reference source not found.). InSAR may thus accurately reflect the landscape-average subsidence, rather than faster subsidence at the specific camera location during the 2019 El Nino. Excluding the 2019 El Nino, InSAR showed high consistency in RMSE and r with the peat camera (Error! Reference source not found., Table 1). Next, ALOS-2’s 14- or 28-day revisit interval was better suited for detecting longer-term trends over months, such as from El Nino dry seasons and persistent drainage effects. Shorter fluctuations were less detectable, such as three rapid uplift episodes at Kalampangan cropland in 2019 which corresponded to irrigation timings (FIGs. 12A and 12B).

[0078] The higher InSAR accuracy estimations given above were achieved through robust noise reduction, which was unaccounted for in previous works. The ionospheric noise was upto 6.7 m, averaging 2.3 ± 1.5 m between any two points across the 350 km x 350 km ALOS-2 image (FIG. 14B). A range split-spectrum correction was utilized which better captured the spatial complexity of ionospheric noise (FIG. 14 A) as opposed to ad-hoc use of linear or quadratic detrending in previous works. The closure phase bias caused up to 30 cm of cumulative false subsidence over 5.5 years, especially over scrubland and forests (FIG. 15A). A larger temporal network of interferograms were used in comparison to previous studies and an additional correction was performed to reduce the bias. Tropospheric noise was up to 15 cm at any given time (Error! Reference source not found.). Corrections using global atmospheric models in previous studies led to worse or no change in the RMSE and r of rates and displacements compared to without correction, while the proposed correction based on spatially correlated noise characteristics improved the performance across all sites (Table ).Table 2: Comparison of InSAR accuracy using different tropospheric correctionsaDifferences in the long-term, multiannual rates were computed from the peat camera velocity minus the InSAR velocity.b> cRMSE and r were computed between epoch-to-epoch displacements of InSAR and peat cameras. InSAR motions were projected from the line-of-sight to vertical direction for all calculations.

[0079] The sub-centimetre accuracy in L-band InSAR long-term rates with corrections applied as outlined above could significantly improve peat carbon emission estimations. Awidely used method by Couwenberg & Hooijer converts subsidence rates to CO2 emissions by relating subsidence to changes in anoxic peat thickness. Based on this approach, the long-term rate accuracy of up to 0.06 cm yr'1translates to 0.97 t CO2 ha'1yr'1of uncertainty in CO2 emissions, assuming if all subsidence is due to peat oxidation. This assumption is appropriate for long-drained peatlands where oxidation dominates subsidence, such as the Mega Rice Project area analysed here, but does not apply to recently drained peatlands where compaction and consolidation dominate subsidence. The sub-centimetre accuracy, at least two orders of magnitude smaller than major InSAR noises unaccounted for in prior works, thus enhanced reliability in estimating CO2 emissions from peat motion. Additionally, this level of accuracy was achieved cost-effectively and over a large area across sites, overcoming practical limitations of ground-based measurements, such as subsidence poles, flux towers, and chambers. Current carbon emission reporting frameworks rely heavily on point -based ground measurements, like in the case of Tier 1 emission factors in the 2013 IPCC Wetlands. Achieving broad coverage with these methods is costly and labour-intensive. Subsidence poles, in particular, require numerous readings to derive an accurate measurement that accounts for human error and variability due to microtopography within a single site. In contrast, the InSAR accuracy demonstrated here may be upscaled with ease and thus offers a practical solution for site- to regional-scale carbon emission estimates.

[0080] Variability and drivers of peat motion

[0081] High variability in long-term rates of subsidence was found within similar land uses and drainage intensities; this underpinned the value of InSAR for refining carbon emissions, since no single environmental factor strongly predicts peat motion. Extending the analysis to all peatlands of Central Kalimantan, it was observed that there is a large overlap and similar mean long-term rates across land uses (FIG. 16C). Correlations of long-term rates with distance from the peat edge (r = 0.16, FIG. 16D) and with drainage canal density (r = -0.20, FIG. 16E)were weak. A possible interpretation is that greater drainage density leads to higher subsidence rates as expected, but that areas close to the peat edge tend to have thinner peat and therefore subside less, or remain stable if most peat had been lost. The weak correlations may reflect a combination of high local-scale variations in land-management, InSAR measurement uncertainties, land use classification errors, or unquantified impacts of water management. For example, it was suggested that differing canal depths and flow control structures likely influenced peat motions more than drainage density in the former Mega Rice Project. As rewetting efforts gain traction, variations in water management could further accentuate peat motion variability. If water management is indeed influencing the InSAR results, this would imply that the method has potential to monitor changes in peat condition due to interventions such as peat re -wetting, though further analysis using data from areas with known water table changes is needed.

[0082] Over shorter timescales, the InSAR displacements provide insights on the peat’s response to disturbances and recovery. During the 2019 El Nino, all land uses experienced subsidence (Error! Reference source not found.), but several swamp, forest, and scrubland areas fully rebounded after the dry season, with no long-term subsidence (Error! Reference source not found.). In contrast, cropland and plantations showed partial rebound, implying irreversible subsidence and long-term peat loss. Within each land use, areas with less subsidence during El Nino do not necessarily experience less long-term subsidence (Error! Reference source not found. A to 17F). These could again be due to localised differences in drainage intensity and water management, or because peat in less degraded sites retains greater elasticity to respond to dry-wet cycles. This elasticity helps natural peatlands stay wetter during droughts and was observed in intact Sumatran forests, while drained plantations showed less short-term variability but faster long-term subsidence. As with high-latitude peatlands, InSAR’ sability to track seasonal and short-term motion in tropical peatlands may serve as a diagnostic tool for assessing peat condition and associated CO2 emissions.

[0083] The high spatial and temporal variabilities in peat motions found above underscore the benefit of direct, location-specific monitoring of heterogeneous tropical peat landscapes through L-band InSAR. Conventional approaches rely on sparse measurements, such as subsidence poles, and use proxies, such as land use or drainage intensities shown above, to extrapolate and estimate CO2 emissions over large areas. These proxies cannot capture fine-scale spatial or temporal variations in peat motion that InSAR reveals and may lead to under-or over-estimation of CO2 emissions. In contrast, the InSAR framework as disclosed herein may provide continuous coverage and direct measurements across heterogenous peatlands, reducing reliance on proxies, generalised assumptions, or static measurements, and enable more accurate and comprehensive emission estimates.

[0084] In conclusion, through the first direct validation against field-measured peat motion in tropical peatlands it was shown that L-band InSAR accurately captures episodic and interannual peat motion trends. The proposed framework targets degraded parts of peatlands, where CO2 emissions are highest, and monitoring is most needed. The results here will benefit peatland MRV, filling gaps in existing approaches that lack spatial and temporal coverage. InSAR also provides a source- and location-specific proxy for belowground peat degradation, complementing broader atmospheric CO2 measurements such as from flux towers and satellite spectrometers. Accuracy and scalability from InSAR are important, especially as Southeast Asia sees a rise in carbon markets such as IDX Carbon (Indonesia) and Climate Impact X (Singapore), where robust peatland MRV is essential for transparency. Because these markets build upon widely established or recognised protocols for carbon emission estimations such as those that have been used for IPCC reporting, integrating InSAR measurements into these protocols would enhance the robustness of emission estimates and improve the integrity of bothscientific reporting and carbon market accounting. With peat management and restoration recognized as key to meeting global CO2 emissions reduction targets, the InSAR framework here offers a cost-effective and scalable solution for supporting these critical efforts.

[0085] Methods

[0086] InSAR analysis

[0087] 89 PALSAR-2 ALOS-2 ScanS AR Ll.l images acquired between 12 October 2018 and 9 August 2024, along Path 28, Frame 3650 were used. To derive displacements, the alosStack module in the InSAR Scientific Computing Environment (ISCE) v2.6.3 software was used to resample images to a common grid using the 19 June 2020 image as the reference, interfere images to form interferograms, perform ionospheric corrections, and multi-look, filter and unwrap interferograms. The interferograms were multi-looked by 17 and 56 in range and azimuth directions, respectively, resulting in pixel sizes of around 180 m x 180 m after geocoding.

[0088] The Miami InSAR time-series software in Python (MintPy) vl.5.3 software was used to perform Small BAseline Subset (SBAS) network inversion to solve for stepwise displacements between images, correct for solid Earth tides (SET) and closure phase biases, and geocode the time series result. Pixels from all interferograms were included and weighted by the inverse of phase covariance for the SBAS network inversion, which ensures that noisy pixels are downweighted. For post -processing analyses relating to land uses, distances from peat edge, and drainage canal densities, pixels with a high mean phase standard deviation > 0.7 cm) were excluded to improve robustness of the analyses. SET refers to the deformation of the Earth’s crust caused by the gravitational attraction of the moon, Sun, and other celestial bodies. To remove these non-peat displacements, MintPy’s SET correction was applied, which uses PySolid to estimate the tidal deformation at each location and image acquisition time based on Love numbers and celestial ephemerides. Finally, the geocoded time series for spatiallycorrelated tropospheric noise were corrected using an independent algorithm, described in a later section. The topographic residual correction in MintPy was not applied as the ALOS-2 orbits are well controlled with typical perpendicular baselines of ALOS-2 data between 100 m and 500 m relative to the reference image (Error! Reference source not found. 18). These result in small height ambiguities where the topographic residual correction is not critical.

[0089] The reference point of the InSAR analysis was collocated with the CPKY Global Navigation Satellite System (GNSS) station located in the urban area of Palangkaraya (-2.220°, 113.926°). The station is stable, with a vertical velocity of <1 mm / year from 2018 to 2024, based on a GNSS time series processed. To ensure a fair comparison between peat cameras and InSAR, the line-of-sight InSAR displacements and rates to the vertical were projected by dividing the cosine of pixel-wise local incidence angles. The parameters which were significant for noise reduction in the derived displacements and rates are detailed in the following sections.

[0090] Ionospheric noise correction

[0091] The ionospheric phase was estimated using the range split-spectrum approach in the alosStack module. Interferograms were heavily multi-looked by 80 and 448 in range and azimuth, respectively, to improve the signal-to-noise ratio, while still maintaining sufficient detail to estimate the long-wavelength ionospheric phase. The correction was most robust over areas which suffer from minimal decorrelation. Thus, only bandwidth-4 interferograms were used which maintain reasonable coherence over most pixels to solve for the ionospheric contribution of each individual image. The bandwidth was defined as the maximum temporal separation used to pair images to form the network of interferograms. In a bandwidth-4 interferogram network, for example, each image is paired with the next 4 images acquired after itself to form 4 different interferograms. The bandwidth is thus agnostic to the number of days between images. Decorrelated areas were removed when estimating the ionospheric phase, andthe spatial interpolation of the estimated ionospheric phase were used to fill gaps over decorrelated areas. These decorrelated areas often occur over rivers and water bodies that may not be captured in the water body dataset from the Shuttle Radar Topography Mission. An example of an improved ionospheric correction due to the removal of decorrelated areas is shown in FIG. 5. FIGs. 14A and 14B show the offsets in displacement caused by ionospheric noise for each date when an image was acquired.

[0092] Closure phase bias correction

[0093] Closure phase bias refers to the discrepancy amongst triplets of phases ( ), where <t>AC ~ 4>AB ~ 4>BC 0 and A, B, and C are different time epochs. The bias is most prominent in a small bandwidth interferogram network likely because of short-lived signals such as vegetation and soil moisture changes that may cause inconsistency within phase triplets. To minimize the effect of closure phase biases, we used (i) a redundant number of interferograms from a bandwidth-20 network, weighted by nJ; for the SBAS network inversion, and (ii) a correction method in MintPy (for example, the method in Zheng, Y., Fattahi, H., Agram, P., Simons, M. & Rosen, P. A. On Closure Phase and Systematic Bias in Multilooked SAR Interferometry. IEEE Trans. Geosci. Remote Sens. 60, 1-11 (2022)). An infinitely large bandwidth is most ideal to minimize the effects of closure phase bias. However, this is not practical for tropical peatlands even with L-band ALOS-2 images, which decorrelate faster with increasing bandwidth due to the longer temporal baselines, commonly of 14, 28 days, or more, as compared to 6 or 12 days for Sentinel- 1. Large bandwidths are also not practical for most areas in the tropical peatlands where coherence is already low and large bandwidth interferograms would be mostly decorrelated. The inclusion of decorrelated interferograms will then have insignificant impact on the SBAS network inversion since the interferograms are weighted by nJ;. The number of bandwidths used were increased up to 20, where convergence of the SBAS network inversion was observed (FIGs. 15A to 15D). This is computationallyefficient and also maximizes the performance of the closure phase bias correction of Zheng et al. As shown in FIGs. 15 A to 15D, the closure phase bias is continually reduced with increasing number of bandwidths used, especially for the initial smaller bandwidths. The correction is also most critical over scrubland and forests which experience the largest bias, followed by cropland, then plantations (FIGs. 15A to 15D). Finally, the quality of the estimated closure phase depends on the coherence and cannot be reliably estimated if the data is decorrelated. Different levels of closure phase bias correction were applied for each pixel to avoid introducing unwrapping errors during the correction - using smaller bandwidths to estimate the bias for pixels which decorrelate faster, and larger bandwidths for those which decorrelate slower.

[0094] Spatially correlated tropospheric noise correction

[0095] Tropospheric noise may have magnitudes up to tens of centimetres. It is spatially correlated over varying distances of few to tens of kilometres, but temporally uncorrelated over a few hours or less. The wet (turbulent) tropospheric component, driven by differential spatial variations in water vapour, contributes noise across a broad range of spatial scales from a few to tens of kilometres, as turbulent processes cascade energy from large atmospheric motions to smaller scales. The dry (hydrostatic) component varies more smoothly with pressure and temperature, and can also be modulated by elevation-dependent water vapour. In the study area, spatially correlated tropospheric noise is significant as observed from anomalous values of up to 15 cm in the InS AR time series that occur at the same time epoch across several InS AR pixels that are located close to each other (FIGs. 9A and 9B). For example, Hampangen plantation and UPR forest are just a few kilometres apart, while the other 6 sites form a separate, closely located group. Each group shares the same peaks in their time series (FIGs. 9A and 9B).

[0096] It was found that corrections based on global atmospheric models (GAMs) such as PyAPS and GACOS worsened the InSAR accuracy, while corrections based on height correlation were not effective for reducing these tropospheric effects (Table ). This is likelybecause of the low resolution of the atmospheric variables in GAMs which are not suited for removing turbulent tropospheric noise. Since the terrain around the Palangkaraya peatlands is mostly flat, the elevation-dependent tropospheric noise is not significant, as the thickness of the troposphere that the InSAR microwave propagates through is almost constant. Common scene stacking corrections which rely on temporal averaging of numerous scenes to reduce noise were also inapplicable due to the low frequency of revisits of ALOS-2.

[0097] Instead, the spatially correlated tropospheric noise at each pixel was reduced by utilizing the stacked time series of neighbouring non-deforming pixels. Based on the assumption that tropospheric noise is spatially correlated, the signals at non-deforming pixels ( non-de / ^ciose(0 cac,other would contain predominantly tropospheric noise that are similar across pixels, and to a smaller extent other random noise that may be dissimilar across pixels:> < < < >"phase of a target pixel which we want to remove the tropospheric noise from. Thus, stacking the time series of these neighbouring non-deforming pixels, by taking the mean, would amplify the spatially correlated tropospheric noise while cancelling out other random noise. This results in an approximate of the spatially correlated tropospheric noise around a target pixel that can then be removed from the target pixel itself:

[0099] This noise correction was applied for each pixel independently. Non-deforming neighbouring pixels were defined as those with a velocity of less than 1 mm / year over the whole time series, an average spatial coherence above 0.3, and a maximum distance of 10 km awayfrom the target pixel. The maximum distance was selected to be within typical spatial correlation lengths of wet tropospheric noise of not more than tens of kilometres. Adaptive distances for each target pixel based on spatial correlation metrics may be alternatively applied. However, in this instance, it was found that a maximum 10 km distance is sufficient for removing large-magnitude tropospheric signals in this study area (Error! Reference source not found. 9A and 9B). Finally, only target pixels with at least 5 non-deforming neighbouring pixels were corrected for the spatially correlated tropospheric noise to avoid the inclusion of other random noise. The correction improved the average RMSE and r between InSAR and peat camera epoch-to-epoch displacements, while the RMSE of long-term rates was mostly unaffected (Table 2). The reduced noise in the time series is also shown by the fewer number of anomalous values for all sites (Error! Reference source not found. 9A and 9BError! Reference source not found.).

[0100] Phase variance and standard deviation

[0101] The phase variance, nJ; or standard deviation,were used to indicate the quality of InSAR measurements (Table 1 column f, FIGs. 7A and 7B). nJ; was estimated from the coherence, y for each InSAR pixel following Equation 7, where y refers to the similarity between two SAR images acquired on different dates:where y is higher, nJ;are lower, and the InSAR phase is more similar between images, facilitating more reliable displacement tracking,of each date in the displacement time series (Error! Reference source not found. 9A and 9B) was estimated via linear propagation of <T of each interferogram used during the SB AS network inversion. The meanin Table 1 column f was then derived from the mean of allof each date over the entire displacement time series.

[0102] InSAR captured the episodic fluctuations in peat motion better at peat camera sites when <T is lower. The meanranged from 0.37 cm to 0.73 cm across non-forested sites where InSAR performed better and exceeded 0.80 cm at forest sites where InSAR performed worse (Table 1 column f, FIGs. 7 A and 7B) due to several expected causes of coherence loss: Firstly, forest sites have dense, multi-layered vegetation structures that typically cause strong volume scattering and thus coherence loss. Secondly, vegetation growth causes temporal decorrelation while also increasing volume scattering. The South scrubland (meanof 0.73 cm) for example, was initially covered in short vegetation in late 2018, but later consisted of trees and taller vegetation in later years due to forest regrowth. Thirdly, longer than usual time gaps between images further decrease coherence at any site, such as from mid 2023 to early 2024 where the area was imaged every 42 to 84 days instead of 14 to 28 days,is thus expected to be higher at forest sites and to varying degrees, scrubland, compared to the plantation and cropland sites. Similar variations inalso extend beyond the peat camera sites across the analysed peat area in Central Kalimantan, where the meanwas often higher over forests than other land uses (Error! Reference source not found. 19). However,still varies widely within each land use (Error! Reference source not found. 19) and is a more reliable indicator of the InSAR measurement quality than the generic land use alone.

[0103] Conversion of peat motions to carbon emissions

[0104] The method of Couwenberg, J. & Hooijer, A. Towards robust subsidence-based soil carbon emission factors for peat soils in south-east Asia, with special reference to oil palm plantations. Mires Peat 12, (2013), which has been applied in several tropical peatland studies, was used to approximate the amount of carbon emissions (t CO2 ha'1yr'1) from peat motion rates as follows:CO emission —(8) subsidence rate x dry bulk density x carbon concentration x 3.67

[0105] Long-term InSAR subsidence rates was used, i.e. carbon concentration of 55% and dry bulk density of 0.08 g cm'3which are mean values across Southeast Asia peatlands, and a C to CO2 conversion factor of 3.67. The equation assumes that all subsidence is due to peat oxidation but may be scaled if a portion of the total subsidence should be attributed to compaction and consolidation. A higher percentage of subsidence due to oxidation would be more reflective of peat that was drained longer before, where most compaction and consolidation have already occurred during initial stages of drainage. For simplicity, all measured subsidence was assumed to reflect oxidation, which is appropriate for long-drained areas such as the study area but not for recently drained areas where compaction and consolidation will dominate.

[0106] Analysis of environmental factors

[0107] 2023 maps for land uses obtained from the Ministry of Environment and Forestry of Indonesia were referred to. For simplicity, land uses were aggregated into broad categories of swamp, forest (primary and secondary dry land, mangrove and swamp forests, and plantation forests), open shrub / scrub (open land, shrubs, and swamp shrubs), smallholder agriculture (dry land and mixed dry land farming, and rice fields), plantation, and urban area (port, transmigration, and settlement areas). Actual land uses at the peat cameras are based on field observations and may differ from the land use maps.

[0108] The boundary of peatlands in Central Kalimantan within the ALOS-2 image and distances of each InSAR pixel to the peat edge were determined from peat maps as of 2019 from the Indonesian Soil Research Institute.

[0109] To compute drainage canal densities, published drainage canal maps derived from optical imagery acquired by PlanetScope CubeSat sensors between July and September 2017 were used and processed using a convolutional neural network. The maps were converted from a binary raster format to a line vector shapefile format, then used the shapefile to sum the lengthof drainage canals at a resolution of 1 km to obtain drainage canal densities in units of km per km2- following a similar convention used in the original publication. Finally, the lower-resolution (1 km) raster of drainage canal densities were resampled to the same grid as the higher-resolution (180 m) raster of InSAR peat motion rates to compute the pixel-to-pixel correlation between the two.

[0110] The proposed system for carbon emission mapping may be implemented by a processor system 900 as illustrated in the schematic block diagram of FIG. 20. Components of the processing system 900 may be provided within one or more computing device to carry out the functions of the modules or any other modules. One skilled in the art will recognize that the exact configuration or arrangement illustrated in FIG. 20 is provided by way of example only, e.g., each processing system provided may be different and the exact configuration of processing system 900 may vary.

[0111] In embodiments of the present disclosure, the processing system 900 may include a controller 901 and user interface 902. User interface 902 is configured to enable manual interactions between a user and the computing module as required. For this purpose, the processing system 900 includes the input / output components required for the user to enter instructions to provide updates to each of the modules. A person skilled in the art will recognize that components of user interface 902 may vary from embodiment to embodiment but may typically include one or more input devices 935 such as but not limited to a touchscreen, a keyboard, a joystick, a mouse, a microphone, etc. The user interface 902 can also include a media player 940, which can be in the form of one or more playback devices, including but not limited to a display, a speaker, earphones, headsets, etc.

[0112] The controller 901 is configured to be in data communication with the user interface 902 via bus 915. The controller 901 includes memory 920 and processor 905 mounted on a circuit board to process instructions and data, e.g., to perform the method of the presentdisclosure. The controller 901 includes an operating system 906, an input / output (I / O) interface 930 for communicating with user interface 902, and a communications interface, e.g., a network card 950. The network card 950 may, for example, be configured to send data from the controller 901 via a wired or wireless network to other processing devices or to receive data via the wired or wireless network. Wireless networks that may be utilized by the network card 950 include, but are not limited to, Wireless-Fidelity (Wi-Fi), Bluetooth, Near Field Communication (NFC), cellular networks, satellite networks, telecommunication networks, Wide Area Networks (WAN), and etc.

[0113] Memory 920 and operating system 906 are in data communication with central processing unit (CPU) 905 via bus 910. The memory 920 may include both volatile and nonvolatile memory. The memory 920 may include more than one of each type of memory, e.g., Random Access Memory (RAM) 923, Read Only Memory (ROM) 925, and a mass storage device 945. The mass storage device 945 may include one or more solid-state drives (SSDs). One skilled in the art will recognize that the memory described above includes non-transitory computer-readable media and shall be taken to include all computer-readable media except for a transitory, propagating signal. Typically, instructions are stored as program code in the memory but can also be hardwired. Memory 920 may include a kernel and / or programming modules such as a software application that may be stored in either volatile or non-volatile memory.

[0114] Herein, the term “processor” is used to refer generically to any device or component that can process computer-readable instructions, including for example, a microprocessor, microcontroller, programmable logic device, or other computational device. That is, processor 905 may be provided by any suitable logic circuitry for receiving inputs, processing them in accordance with instructions stored in memory, and generating outputs (for example to the memory components or media player 940). In the present disclosure, processor 905 may be asingle core or multi-core processor with memory addressable space. In one example, processor 905 may be multi-core, comprising — for example — an 8 core CPU. In another example, it could be a cluster of CPU cores operating in parallel to accelerate computations.

[0115] Further, one skilled in the art will recognize that certain functional units in this description have been labelled as modules throughout the specification. The person skilled in the art will also recognize that a module may be implemented as circuits, logic chips or any sort of discrete component. Still further, one skilled in the art will also recognize that a module may be implemented in software which may then be executed by a variety of processor architectures. In embodiments of the disclosure, a module may also comprise computer instructions or executable code that may instruct a computer processor to carry out a sequence of events based on instructions received. In further embodiments, the module may comprise a combination of different types of modules or sub-modules. The choice of the implementation of the modules may be determined by a person skilled in the art and does not limit the scope of the claimed subject matter in any way.

[0116] All examples described herein, whether of apparatus, methods, materials, or products, are presented for the purpose of illustration and to aid understanding, and are not intended to be limiting or exhaustive. Modifications may be made by one of ordinary skill in the art without departing from the scope of the invention as claimed.

Claims

CLAIMS1. A method of carbon emission mapping, the method comprising:obtaining a plurality of first interferograms based on a plurality of radar images of a peatland, wherein the plurality of radar images corresponds to radar images of the peatland captured over a plurality of time instances; estimating an ionospheric phase contribution using an ionospheric correction on the plurality of first interferograms, wherein the plurality of first interferograms corresponds to a first bandwidth configuration; using the ionospheric phase contribution, applying the ionospheric correction to a plurality of second interferograms, wherein the plurality of second interferograms corresponds to a second bandwidth configuration, wherein the second bandwidth is larger than the first bandwidth;performing Small BAseline Subset (SBAS) network inversion on the plurality of second interferograms to obtain a time series data;applying a closure phase bias correction to the time series data;applying a tropospheric noise correction to the time series data to obtain a ground displacement data; andgenerating a carbon emission map of the peatland based on the ground displacement data.

2. The method as recited in claim 1, wherein estimating the ionospheric phase contribution comprises:removing decorrelated pixels in each of the plurality of first interferograms; and interpolating gaps corresponding to the decorrelated pixels in each of the plurality of first interferograms, using the ionospheric phase estimation.

3. The method as recited in claim 2, wherein removing decorrelated pixels in each of the plurality of first interferograms comprises: applying a respective mask to the decorrelated pixels in each of the plurality of first interferograms.

4. The method as recited in any one of claims 2 to 3, wherein estimating the ionospheric phase contribution further comprises:applying a split-spectrum ionospheric correction to each of the plurality of first interferograms, wherein each of the plurality of first interferograms is split into two sub-bands, and wherein a respective ionospheric phase delay is estimated based on frequency-dependent phase differences between the two sub-bands.

5. The method as recited in any one of the above claims, wherein each of the plurality of first interferograms corresponds to unique pairs of the plurality of radar images obtained based on the first bandwidth configuration.

6. The method as recited in claim 5,wherein each of the plurality of second interferograms corresponds to unique pairs of the plurality of radar images obtained based on the second bandwidth configuration.

7. The method as recited in any of the above claims, wherein applying the closure phase bias correction to the time series data further comprises:determining a plurality of closure phase biases based on the time series data;for each pixel in the time series data, applying a pixel-by-pixel closure phase bias correction using a selected one of the plurality of closure phase biases to obtain the corrected time series data.

8. The method as recited in any claim 7, wherein applying the closure phase bias correction to the time series data further comprises:determining instances of unwrapping error for the corresponding pixel in the time series data; anddetermining the selected one of the plurality of closure phase biases corresponding to no unwrapping error.

9. The method as recited in any claim 8, wherein applying the closure phase bias correction to the time series data further comprises:determining the selected one of the plurality of closure phase biases corresponding to a largest temporal separation.

10. The method as recited in any of the above claims, wherein applying the tropospheric noise correction to the time series data comprises:for each target pixel in each of the time series data, correcting the target pixel with a respective spatially correlated tropospheric noise to obtain the ground displacement data.

11. The method as recited in any claim 10, wherein applying the tropospheric noise correction to the time series data further comprises:for each target pixel in each of the time series data, determining the respective spatially correlated tropospheric noise based on a plurality of neighbour pixels around the target pixel.

12. The method as recited in any claim 11, wherein applying the tropospheric noise correction to the time series data further comprises:for each target pixel in each of the time series data, determining the respective spatially correlated tropospheric noise around the target pixel based on a mean of stacked time series of the plurality of neighbour pixels of the target pixel.

13. The method as recited in any one of claims 10 to 12, wherein applying the tropospheric noise correction to the time series data further comprises:determining at least one target pixel from the time series data, wherein the at least one target pixel comprises a plurality of non-deforming neighbour pixels, wherein a count of the plurality of non-deforming neighbour pixels is equal to or above a minimum neighbour pixel threshold.

14. The method as recited in any claim 13, wherein the minimum neighbour pixel threshold is 5.

15. The method as recited in any one of claims 13 and 14, wherein each of the plurality of non-deforming neighbour pixels is characterized by: a velocity under a velocity threshold; an average spatial coherence above a spatial coherence threshold; and a distance away from the target pixel under a maximum distance threshold.

16. The method as recited in any one of the above claims, further comprising:performing multi-look to each of the plurality of first interferograms.

17. The method as recited in claim 16, further comprising:performing multi-look to each of the plurality of second interferograms.

18. The method as recited in claim 17, wherein the number of multi-look performed to each of the plurality of first interferograms is larger than the number of multi-look performed to each of the plurality of second interferograms.

19. The method as recited in any one of the above claims, wherein the plurality of corrected interferograms are weighted.

20. A carbon emission mapping system, comprising:memory storing instructions; anda processor coupled to the memory and configured to process the stored instructions to implement a module configured to perform the method as recited in any one of the above claims.

Images