**Author details**

Ramasamy Velmurugan\* and Balasubramanian Subramanian Electroplating and Metal Finishing Division, Central Electrochemical Research Institute (CSIR), Tamilnadu, India

\*Address all correspondence to<EMAIL_ADDRESS>
© 2021 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*Physicochemical Approaches for Thin Film Energy Storage Devices through PVD Techniques DOI: http://dx.doi.org/10.5772/intechopen.99473*


**Author details**

Woojeong Kim1 , Kui Hyun Kang1 and Jung-Keun Suh2,3\*

\*Address all correspondence to<EMAIL_ADDRESS>

[11] Wakankar A, Chen Y, Gokarn Y, Jacobson FS. Analytical methods for physicochemical characterization of antibody drug conjugates. MAbs. 2011;**3**(2):161-172. DOI: 10.4161/

Characterization of Biopharmaceuticals Focusing on Antibody Therapeutics

http://dx.doi.org/10.5772/intechopen.79107

35

[12] Fan G, Wang Z, Hao M, Li J. Bispecific antibodies and their applications. Journal of

[13] De Nardis C, Hendriks LJA, Poirier E, Arvinte T, Gros P, Bakker ABH, de Kruif J. A new approach for generating bispecific antibodies based on a common light chain format and the stable architecture of human immunoglobulin G(1). The Journal of Biological

[14] Godar M, de Haard H, Blanchetot C, Rasser J. Therapeutic bispecific antibody formats: A patent applications review (1994-2017). Expert Opinion on Therapeutic Patents.

[15] Strohl WR. Current progress in innovative engineered antibodies. Protein & Cell.

[16] Johnson JA. Biologics and Biosimilars: Background and Key Issues. CRS Report.

[17] Crespi-Lofton J, Skelton JB. The growing role of biologics and biosimilars in the United States: Perspectives from the APhA Biologics and Biosimilars Stakeholder Conference. Journal of the American Pharmaceutical Association. 2017;**57**(5):e15-e27. DOI: 10.1016/j.

[18] Peterson J, Budlong H, Affeldt T, Skiermont K, Kyllo G, Heaton A. Biosimilar products in the modern U.S. health care and regulatory landscape. Journal of Managed Care &

[19] U.S. Food and Drug Administration. Biosimilar and Interchangeable Products [Internet]. Available from: https://www.fda.gov/Drugs/DevelopmentApprovalProcess/ HowDrugsareDevelopedandApproved/ApprovalApplications/TherapeuticBiologic-

[20] Song KE, Byeon J, Moon DB, Kim HH, Choi Y-J, Suh J-K. Structural identification of modified amino acids on the interface between EPO and its receptor from EPO BRP, human recombinant erythropoietin by LC/MS analysis. Molecules and Cells. 2014;**37**(11):819-

[21] Byeon J, Choi Y-J, Suh J-K. Structural characterization of modification on the interface between a ligand and its receptor for biopharmaceuticals. Receptors & Clinical

[22] Byeon J, Kang KH, Jung H-K, Suh J-K. Assessment for quantification of biopharmaceutical protein using a microvolume spectrometer on microfluidic slides. BioChip Journal.

Specialty Pharmacy. 2017;**23**(12):1255-1259. DOI: 10.18553/jmcp.2<IP_ADDRESS>5

Hematology & Oncology. 2015;**8**:130. DOI: 10.1186/s13045-015-0227-0

Chemistry. 2017;**292**(35):14706-14717. DOI: 10.1074/jbc.M117.793497

2017;R44620. Available from: https://fas.org/sgp/crs/misc/R44620.pdf

Applications/Biosimilars/ucm580419.htm [Accessed: 15-04-2018]

826. DOI: 10.14348/molcells.2014.0214

Investigation. 2015;**2**(2):1-5. DOI: 10.14800/rci.536

2017;**11**(1):21-29. DOI: 10.1007/s13206-016-1104-9

2018;**28**(3):251-276. DOI: 10.1080/13543776.2018.1428307

2018;**9**(1):86-120. DOI: 10.1007/s13238-017-0457-8

mabs.3.2.14960

japh.2017.05.014

3 Seoul Media Institute of Technology, Seoul, Republic of Korea

### **References**

[11] Wakankar A, Chen Y, Gokarn Y, Jacobson FS. Analytical methods for physicochemical characterization of antibody drug conjugates. MAbs. 2011;**3**(2):161-172. DOI: 10.4161/ mabs.3.2.14960

**Author details**

34 Biopharmaceuticals

Woojeong Kim1

**References**

10.1016/j.jpba.2016.05.028

DOI: 10.1038/nbt.3040

0929867324666170113112759

Pharmacy. 1989;**46**(11 Suppl 2):S9-S11

2006;**17**:638-642. DOI: 10.1016/j.copbio.2006.10.002

, Kui Hyun Kang1

\*Address all correspondence to<EMAIL_ADDRESS>
3 Seoul Media Institute of Technology, Seoul, Republic of Korea

and Jung-Keun Suh2,3\*

1 BIOnSYSTEMS, Research Center, Seongnam-si, Gyeonggi-do, Republic of Korea

2 Immersive Media Lab, Seoul Media Institute of Technology, Seoul, Republic of Korea

[1] Parr MK, Montacir O, Montacir H. Physicochemical characterization of biopharmaceuticals. Journal of Pharmaceutical and Biomedical Analysis. 2016;**130**:366-389. DOI:

[2] The MJ. Human insulin: DNA technology's first drug. American Journal of Hospital

[3] Walsh G. Biopharmaceutical benchmarks 2014. Nature Biotechnology. 2014;**32**:992-100.

[4] Sato AK, Viswanathan M, Kent RB, Wood CR. Therapeutic peptides: Technological advances driving peptides into development. Current Opinion in Biotechnology.

[5] Uhlig T, Kyprianou T, Martinelli FG, Oppici CA, Heiligers D, Hills D, Calvo XR, Verhaert P. The emergence of peptides in the pharmaceuticalbusiness: From exploration to exploitation. EuPA Open Proteomics. 2014;**4**:58-69. DOI: 10.1016/j.euprot.2014.05.003

[6] Jafari R, Zolbanin NM, Rafatpanah H, Majidi J, Kazemi T. Fc-fusion proteins in therapy: An updated view. Current Medicinal Chemistry. 2017;**24**(12):1228-1237. DOI: 10.2174/

[7] Rogstad S, Faustino A, Ruth A, Keire D, Boyne M, Park J.A retrospective evaluation of the use of mass spectrometry in FDA biologics license applications. Journal of the American

[8] Top 20 Drugs in the World 2017 [Internet]. Available from: https://igeahub.com/2017/08/

[9] Zolot RS, Basu S, Million RP. Antibody-drug conjugates. Nature Reviews. Drug

[10] Carter PJ, Lazar GA. Next generation antibody drugs: Pursuit of the 'high-hanging fruit. Nature Reviews. Drug Discovery. 2018;**17**(3):197-223. DOI: 10.1038/nrd.2017.227

Society for Mass Spectrometry. 2017;**28**:786-794. DOI: 10.1007/s13361-016-1531-9

08/top-20-drugs-in-the-world-2017/ [Accessed: 01-04-2018]

Discovery. 2013;**12**(4):259-260. DOI: 10.1038/nrd3980

[23] International Conference on Harmonisation of Technical Requirements for Registration of Pharmaceuticals for Human Use. ICH Q5E guideline, Comparability of biotechnological/biological products subject to changes in their manufacturing process. 2004

[35] Liu L. Antibody glycosylation and its impact on the pharmacokinetics and pharmacodynamics of monoclonal antibodies and Fc-fusion proteins. Journal of Pharmaceutical

Characterization of Biopharmaceuticals Focusing on Antibody Therapeutics

http://dx.doi.org/10.5772/intechopen.79107

37

[36] Zhang Z, Pan H, Chen X. Mass spectrometry for structural characterization of therapeutic antibodies. Mass Spectrometry Reviews. 2009;**28**(1):147-176. DOI: 10.1002/mas.20190

[37] Ayoub D, Jabs W, Resemann A, Evers W, Evans C, Main L, Baessmann C, Wagner-Rousset E, Suckau D, Beck A. Correct primary structure assessment and extensive glyco-profiling of cetuximab by a combination of intact, middle-up, middle-down and bottom-up ESI and MALDI mass spectrometry techniques. MAbs. 2013;**5**(5):699-710.

[38] Wang B, Gucinski AC, Keire DA, Buhse LF, Boyne MT 2nd. Structural comparison of two anti-CD20 monoclonal antibody drug products using middle-down mass spectrom-

[39] Fornelli L, Ayoub D, Aizikov K, Beck A, Tsybin YO. Middle-down analysis of monoclonal antibodies with electron transfer dissociation orbitrap fourier transform mass spectrometry. Analytical Chemistry. 2014;**86**(6):3005-3012. DOI: 10.1021/ac4036857 [40] Hartman PA, Stodola JD, Harbour GC, Hoogerheide JG. Reversed-phase high-performance liquid chromatography peptide mapping of bovine somatotropin. Journal of

[41] Recny MA, Scoble HA, Kim Y. Structural characterization of natural human urinary and recombinant DNA-derived erythropoietin. Identification of des-arginine 166 erythropoietin. The Journal of Biological Chemistry. 1987;**262**(35):17156-17163. Available from:

[42] Dougherty JJ Jr, Snyder LM, Sinclair RL, Robins RH. High-performance tryptic mapping of recombinant bovine somatotropin. Analytical Biochemistry. 1990;**190**(1):7-20. DOI:

[43] Schär M, Börnsen KO, Gassmann E. Fast protein sequence determination with matrixassisted laser desorption and ionization mass spectrometry. Rapid Communications in

[44] Bongers J, Cummings JJ, Ebert MB, Federici MM, Gledhill L, Gulati D, Hilliard GM, Jones BH, Lee KR, Mozdzanowski J, Naimoli M, Burman S. Validation of a peptide mapping method for a therapeutic monoclonal antibody: What could we possibly learn about a method we have run 100 times? Journal of Pharmaceutical and Biomedical Analysis.

[45] Hernandez P, Müller M, Appel RD. Automated protein identification by tandem mass spectrometry: Issues and strategies. Mass Spectrometry Reviews. 2006;**25**(2):235-254.

[46] Zhang Y, Fonslow BR, Shan B, Baek MC, Yates JR 3rd. Protein analysis by shotgun/bottom-up proteomics. Chemical Reviews. 2013;**113**(4):2343-2494. DOI: 10.1021/cr3003533

Mass Spectrometry. 1991;**5**(7):319-326. DOI: 10.1002/rcm.1290050705

2000;**21**(6):1099-1128. DOI: 10.1016/S0731-7085(99)00181-8

Chromatography. 1986;**360**(2):385-395. DOI: 10.1016/S0021-9673(00)91687-X

http://www.jbc.org/content/262/35/17156.full.pdf

10.1016/0003-2697(90)90126-T

DOI: 10.1002/mas.20068

etry. The Analyst. 2013;**138**(10):3058-3065. DOI: 10.1039/c3an36524g

Sciences. 2015;**104**(6):1866-1884. DOI: 10.1002/jps.24444

DOI: 10.4161/mabs.25423

[35] Liu L. Antibody glycosylation and its impact on the pharmacokinetics and pharmacodynamics of monoclonal antibodies and Fc-fusion proteins. Journal of Pharmaceutical Sciences. 2015;**104**(6):1866-1884. DOI: 10.1002/jps.24444

[23] International Conference on Harmonisation of Technical Requirements for Registration of Pharmaceuticals for Human Use. ICH Q5E guideline, Comparability of biotechnological/biological products subject to changes in their manufacturing process. 2004 [24] International Conference on Harmonisation of Technical Requirements for Registration of Pharmaceuticals for Human Use. ICH Q6B guideline, Specifications: Test Procedures

[25] International Conference on Harmonisation of Technical Requirements for Registration of Pharmaceuticals for Human Use. ICH Q8 R2 guideline, Pharmaceutical development. 2009. Available from: https://www.ich.org/fileadmin/Public\_Web\_Site/ICH\_Products/

[26] Li Y, Champion MM, Sun L, Champion PA, Wojcik R, Dovichi NJ. Capillary zone electrophoresis-electrospray ionization-tandem mass spectrometry as an alternative proteomics platform to ultraperformance liquid chromatography-electrospray ionization-tandem mass spectrometry for samples of intermediate complexity. Analytical

[27] Srebalus Barnes CA, Lim A. Applications of mass spectrometry for the structural characterization of recombinant protein pharmaceuticals. Mass Spectrometry Reviews.

[28] Bondarenko PV, Second TP, Zabrouskov V, Makarov AA, Zhang Z. Mass measurement and top-down HPLC/MS analysis of intact monoclonal antibodies on a hybrid linear quadrupole ion trap-Orbitrap mass spectrometer. Journal of the American Society for

[29] Tran BQ, Barton C, Feng J, Sandjong A, Yoon SH, Awasthi S, Liang T, Khan MM, Kilgour DPA, Goodlett DR, Goo YA. Comprehensive glycosylation profiling of IgG and IgGfusion proteins by top-down MS with multiple fragmentation techniques. Journal of

[30] He L, Anderson LC, Barnidge DR, Murray DL, Hendrickson CL, Marshall AG.Analysis of monoclonal antibodies in human serum as a model for clinical monoclonal Gammopathy by use of 21 tesla FT-ICR top-down and middle-down MS/MS. Journal of the American Society for Mass Spectrometry. 2017;**28**(5):827-838. DOI: 10.1007/s13361-017-1602-6 [31] Thompson NJ, Rosati S, Heck AJ. Performing native mass spectrometry analysis on therapeutic antibodies. Methods. 2014;**65**(1):11-17. DOI: 10.1016/j.ymeth.2013.05.003 [32] Beck A, Debaene F, Diemer H, Wagner-Rousset E, Colas O, Van Dorsselaer A, Cianférani S. Cutting-edge mass spectrometry characterization of originator, biosimilar and biobetter antibodies. Journal of Mass Spectrometry. 2015;**50**(2):285-297. DOI: 10.1002/jms.3554

[33] Siuti N, Kelleher NL. Decoding protein modifications using top-down mass spectrom-

[34] Macht M. Top-down characterization of biopharmaceuticals. Trends in Analytical

etry. Nature Methods. 2007;**4**(10):817-821. DOI: 10.1038/nmeth1097

Chemistry. 2013;**48**:62-71. DOI: 10.1016/j.trac.2013.01.016

Mass Spectrometry. 2009;**20**(8):1415-1424. DOI: 10.1016/j.jasms.2009.03.020

and Acceptance Criteria for Biotechnological/Biological Products. 1999

Guidelines/Quality/ Q8\_R1/Step4/Q8\_R2\_Guideline.pdf

Chemistry. 2012;**84**(3):1617-1622. DOI: 10.1021/ac202899p

Proteomics. 2016;**134**:93-101. DOI: 10.1016/j.jprot.2015.10.021

2007;**26**(3):370-388. DOI: 10.1002/mas.20129

36 Biopharmaceuticals

[47] Byeon J, Yim Y-R, Kim H-H, Suh J-K. Structural identification of a non-glycosylated variant at Ser126 for O-glycosylation site from EPO BRP, human recombinant erythropoietin by LC/MS analysis. Molecules and Cells. 2015;**38**(6):496-505. DOI: 10.14348/ molcells.2015.2256

[59] Motoyama A, Venable JD, Ruse CI, Yates JR 3rd. Automated ultra-high-pressure multidimensional protein identification technology (UHP-MudPIT) for improved peptide identification of proteomic samples. Analytical Chemistry. 2006;**78**(14):5109-5118. DOI:

Characterization of Biopharmaceuticals Focusing on Antibody Therapeutics

http://dx.doi.org/10.5772/intechopen.79107

39

[60] Omamogho JO, Glennon JD. Comparison between the efficiencies of sub-2 μm C18 particles packed in narrow bore columns. Analytical Chemistry. 2011;**83**(5):1547-1556. DOI:

[61] Water TH. Recent innovations in UHPLC columns and instrumentation. Trends in

[62] González-Ruiz V, Olives AI, Martín MA. Core-shell particles lead the way to renewing high-performance liquid chromatography. Trends in Analytical Chemistry. 2015;**64**:17-

[63] Kochling J, Wu W, Hua Y, Guan Q, Castaneda-Merced J. A platform analytical quality by design (AQbD) approach for multiple UHPLC-UV and UHPLC-MS methods development for protein analysis. Journal of Pharmaceutical and Biomedical Analysis.

[64] Ghosh PK. Introduction to Protein Mass Spectrometry. 1st ed. London: Academic Press; 2015. Available from: https://www.sciencedirect.com/science/book/9780128021026 [65] Domon B, Aebersold R. Mass spectrometry and protein analysis. Science. 2006;**312**(5771):

[66] Yu L, Remmele RL Jr, He B. Identification of N-terminal modification for recombinant monoclonal antibody light chain using partial reduction and quadrupole time-of-flight mass spectrometry. Rapid Communications in Mass Spectrometry. 2006;**20**(24):3674-

[67] Wang D, Wynne C, Gu F, Becker C, Zhao J, Mueller HM, Li H, Shameem M, Liu YH. Characterization of drug-product-related impurities and variants of a therapeutic monoclonal antibody by higher energy C-trap dissociation mass spectrometry.

[68] Colgrave ML, Craik DJ. Thermal, chemical, and enzymatic stability of the cyclotide kalata B1: The importance of the cyclic cystine knot. Biochemistry. 2004;**43**(20):5965-

[69] Fass D. Disulfide bonding in protein biophysics. Annual Review of Biophysics. 2012;**41**:

[70] Goyder MS, Rebeaud F, Pfeifer ME, Kálmán F. Strategies in mass spectrometry for the assignment of Cys-Cys disulfide connectivities in proteins. Expert Review of Proteomics.

[71] Wiesner J, Resemann A, Evans C, Suckau D, Jabs W. Advanced mass spectrometry workflows for analyzing disulfide bonds in biologics. Expert Review of Proteomics.

Analytical Chemistry. 2015;**87**(2):914-921. DOI: 10.1021/ac503158g

63-79. DOI: 10.1146/annurev-biophys-050511-102321

2013;**10**(5):489-501. DOI: 10.1586/14789450.2013.837663

2015;**12**(2):115-123. DOI: 10.1586/14789450.2015.1018896

Analytical Chemistry. 2014;**63**:14-20. DOI: 10.1016/j.trac.2014.07.016

10.1021/ac060354u

10.1021/ac102139a

28. DOI: 10.1016/j.trac.2014.08.008

212-217. DOI: 10.1126/science.1124619

3680. DOI: 10.1002/rcm.2790

5975. DOI: 10.1021/bi049711q

2016;**125**:130-139. DOI: 10.1016/j.jpba.2016.03.031

[59] Motoyama A, Venable JD, Ruse CI, Yates JR 3rd. Automated ultra-high-pressure multidimensional protein identification technology (UHP-MudPIT) for improved peptide identification of proteomic samples. Analytical Chemistry. 2006;**78**(14):5109-5118. DOI: 10.1021/ac060354u

[47] Byeon J, Yim Y-R, Kim H-H, Suh J-K. Structural identification of a non-glycosylated variant at Ser126 for O-glycosylation site from EPO BRP, human recombinant erythropoietin by LC/MS analysis. Molecules and Cells. 2015;**38**(6):496-505. DOI: 10.14348/

[48] Wang Y, Li X, Liu YH, Richardson D, Li H, Shameem M, Yang X. Simultaneous monitoring of oxidation, deamidation, isomerization, and glycosylation of monoclonal antibodies by liquid chromatography-mass spectrometry method with ultrafast tryptic

[49] Pace CN, Vajdos F, Fee L, Grimsley G, Gray T. How to measure and predict the molar absorption coefficient of a protein. Protein Science. 1995;**4**(11):2411-2423. DOI: 10.1002/

[50] Proc JL, Kuzyk MA, Hardie DB, Yang J, Smith DS, Jackson AM, Parker CE, Borchers CH. A quantitative study of the effects of chaotropic agents, surfactants, and solvents on the digestion efficiency of human plasma proteins by trypsin. Journal of Proteome

[51] Yu YQ, Gilar M, Lee PJ, Bouvier ES, Gebler JC. Enzyme-friendly, mass spectrometrycompatible surfactant for in-solution enzymatic digestion of proteins. Analytical

[52] Han JC, Han GY. A procedure for quantitative determination of tris(2-carboxyethyl) phosphine, an odorless reducing agent more stable and effective than dithiothreitol.

[53] Svoboda M, Meister W, Vetter W. A method for counting disulfide bridges in small proteins by reduction with mercaptoethanol and electrospray mass spectrometry. Journal of

[54] Switzar L, Giera M, Niessen WM. Protein digestion: An overview of the available techniques and recent developments. Journal of Proteome Research. 2013;**12**(3):1067-1077.

[55] Tsiatsiani L, Heck AJ. Proteomics beyond trypsin. The FEBS Journal. 2015;**282**(14):2612-

[56] Bunkenborg J, Espadas G, Molina H. Cutting edge proteomics: Benchmarking of six commercial trypsins. Journal of Proteome Research. 2013;**12**(8):3631-3641. DOI: 10.1021/

[57] Giansanti P, Tsiatsiani L, Low TY, Heck AJ. Six alternative proteases for mass spectrometry-based proteomics beyond trypsin. Nature Protocols. 2016;**11**(5):993-1006. DOI:

[58] Plumb R, Castro-Perez J, Granger J, Beattie I, Joncour K, Wright A. Ultra-performance liquid chromatography coupled to quadrupole-orthogonal time-of-flight mass spectrometry. Rapid Communications in Mass Spectrometry. 2004;**18**(19):2331-2337. DOI:

Analytical Biochemistry. 1994;**220**(1):5-10. DOI: 10.1006/abio.1994.1290

Mass Spectrometry. 1995;**30**:1562-1566. DOI: 10.1002/jms.1190301107

digestion. MAbs. 2016;**8**(8):1477-1486. DOI: 10.1080/19420862.2016.1226715

Research. 2010;**9**(10):5422-5437. DOI: 10.1021/pr100656u

Chemistry. 2003;**75**(21):6023-6028. DOI: 10.1021/ac0346196

molcells.2015.2256

38 Biopharmaceuticals

pro.5560041120

DOI: 10.1021/pr301201x

10.1038/nprot.2016.057

10.1002/rcm.1627

pr4001465

2626. DOI: 10.1111/febs.13287

[72] Gorman JJ, Wallis TP, Pitt JJ. Protein disulfide bond determination by mass spectrometry. Mass Spectrometry Reviews. 2002;**21**(3):183-216. DOI: 10.1002/mas.10025

[83] Berkowitz SA, Engen JR, Mazzeo JR, Jones GB. Analytical tools for characterizing biopharmaceuticals and the implications for biosimilars. Nature Reviews. Drug Discovery.

Characterization of Biopharmaceuticals Focusing on Antibody Therapeutics

http://dx.doi.org/10.5772/intechopen.79107

41

[84] Diepold K, Bomans K, Wiedmann M, Zimmermann B, Petzold A, Schlothauer T, Mueller R, Moritz B, Stracke JO, Mølhøj M, Reusch D, Bulau P. Simultaneous assessment of Asp isomerization and Asn deamidation in recombinant antibodies by LC-MS following incubation at elevated temperatures. PLoS One. 2012;**7**(1):e30295. DOI: 10.1371/ journal.

[85] Beck A, Wagner-Rousset E, Ayoub D, Van Dorsselaer A, Sanglier-Cianferani S. Characterization of therapeutic antibodies and related products. Analytical Chemistry. 2013;**85**:715-

[86] Fekete S, Guillarme D, Sandra P, Sandra K. Chromatographic, electrophoretic, and mass spectrometric methods for the analytical characterization of protein biopharmaceuticals.

[87] Rogers RS, Nightlinger NS, Livingston B, Campbell P, Bailey R, Balland A. Development of a quantitative mass spectrometry multiattribute method for characterization, quality control testing and disposition of biologics. MAbs. 2015;**7**:881-890. DOI: 10.1080/19420862

[88] Wang T, Chu L, Li W, Lawson K, Apostol I, Eris T. Application of a quantitative LC–MS multiattribute method for monitoring site-specific glycan heterogeneity on a monoclonal antibody containing two N-linked glycosylation sites. Analytical Chemistry.

[89] Huhn C, Selman MH, Ruhaak LR, Deelder AM, Wuhrer M. IgG glycosylation analysis.

[90] Abès R, Teillaud JL. Impact of glycosylation on effector functions of therapeutic

[91] Sibéril S, Dutertre CA, Fridman WH, Teillaud JL, Fcgamma R. The key to optimize therapeutic antibodies? Critical Reviews in Oncology/Hematology. 2007;**62**(1):26-33. DOI:

[92] Takahashi M, Kuroki Y, Ohtsubo K, Taniguchi N. Core fucose and bisecting GlcNAc, the direct modifiers of the N-glycan core: Their functions and target proteins. Carbohydrate

[93] Abès R, Dutertre CA, Agnelli L, Teillaud JL. Activating and inhibitory Fcgamma receptors in immunotherapy: Being the actor or being the target. Expert Review of Clinical

[94] Zhang Q, Li Z, Wang Y, Zheng Q, Li J. Mass spectrometry for protein sialoglycosylation.

[95] Tharmalingam T, Adamczyk B, Doherty MA, Royle L, Rudd PM. Strategies for the profiling, characterisation and detailed structural analysis of N-linked oligosaccharides.

Glycoconjugate Journal. 2013;**30**(2):137-146. DOI: 10.1007/s10719-012-9443-9

IgG. Pharmaceuticals (Basel). 2010;**3**(1):146-157. DOI: 10.3390/ph3010146

Research. 2009;**344**(12):1387-1390. DOI: 10.1016/j.carres.2009.04.031

Immunology. 2009;**5**(6):735-747. DOI: 10.1586/eci.09.57

Mass Spectrometry Reviews. 2017. DOI: 10.1002/mas.21555

Analytical Chemistry. 2016;**88**:480-507. DOI: 10.1021/acs.analchem.5b04561

2017;**89**:3562-3567. DOI: 10.1021/acs.analchem.6b04856

Proteomics. 2009;**9**(4):882-913. DOI: 10.1002/pmic.200800715

2012;**11**:527-540. DOI: 10.1038/nrd3746

pone.0030295

736. DOI: 10.1021/ac3032355

10.1016/ j.critrevonc.2006.12.003

[83] Berkowitz SA, Engen JR, Mazzeo JR, Jones GB. Analytical tools for characterizing biopharmaceuticals and the implications for biosimilars. Nature Reviews. Drug Discovery. 2012;**11**:527-540. DOI: 10.1038/nrd3746

[72] Gorman JJ, Wallis TP, Pitt JJ. Protein disulfide bond determination by mass spectrometry. Mass Spectrometry Reviews. 2002;**21**(3):183-216. DOI: 10.1002/mas.10025

[73] Lakbub JC, Shipman JT, Desaire H. Recent mass spectrometry-based techniques and considerations for disulfide bond characterization in proteins. Analytical and Bioanalytical

[74] Liu H, May K. Disulfide bond structures of IgG molecules: Structural variations, chemical modifications and possible impacts to stability and biological function. MAbs.

[75] Wu SL, Jiang H, Hancock WS, Karger BL. Identification of the unpaired cysteine status and complete mapping of the 17 disulfides of recombinant tissue plasminogen activator using LC-MS with electron transfer dissociation/collision induced dissociation.

[76] Wu SL, Jiang H, Lu Q, Dai S, Hancock WS, Karger BL. Mass spectrometric determination of disulfide linkages in recombinant therapeutic proteins using online LC-MS with electron-transfer dissociation. Analytical Chemistry. 2009;**81**(1):112-122. DOI: 10.1021/

[77] Xie H, Gilar M, Gebler JC. Characterization of protein impurities and site-specific modifications using peptide mapping with liquid chromatography and data independent acquisition mass spectrometry. Analytical Chemistry. 2009;**81**(14):5699-5708. DOI:

[78] Xie H, Chakraborty AB, Chen W. Towards fast mapping of protein disulfide bonds: An integrated workflow for automatic assignment of disulfide pairing. In: Proceedings of Practical Applications of Mass Spectrometry in the Biotechnology Industry; 8-10 September 2010; Marina Del Ray. CA: Mass Spec; 2010. p. 107. Available from: http://

[79] Walsh CT, Garneau-Tsodikova S, Gatto GJ Jr. Protein posttranslational modifications: The chemistry of proteome diversifications. Angewandte Chemie (International Ed. in

[80] Li W, Kerwin JL, Schiel J, Formolo T, Davis D, Mahan A, Benchaar SA. Structural elucidation of post-translational modifications in monoclonal antibodies. In: Schiel JE, Davis DL, Borisov OV, editors. State-of-the-art and Emerging Technologies for Therapeutic Monoclonal Antibody Characterization Volume 2. Biopharmaceutical Characterization: The NIST mAb Case Study. Washington, DC: American Chemical Society; 2015. pp. 119-

[81] Walsh G. Post-translational modifications of protein biopharmaceuticals. Drug Discovery

[82] Brinks V, Hawe A, Basmeleh AH, Joachin-Rodriguez L, Haselberg R, Somsen GW, Jis5oot W, Schellekens H. Quality of original and biosimilar epoetin products. Pharmaceutical

Chemistry. 2018;**410**(10):2467-2484. DOI: 10.1007/s00216-017-0772-1

Analytical Chemistry. 2010;**82**(12):5296-5303. DOI: 10.1021/ac100766r

www.waters.com/webassets/cms/ library/docs/720003751en.pdf

English). 2005;**44**(45):7342-7372. DOI: 10.1002/anie.200501023

Today. 2010;**15**:773-780. DOI: 10.1016/j.drudis.2010.06.009

Research. 2011;**28**(2):386-393. DOI: 10.1007/s11095-010-0288-2

183. DOI: 10.1021/bk-2015-1201.ch003

2012;**4**(1):17-23. DOI: 10.4161/mabs.4.1.18347

ac801560k

40 Biopharmaceuticals

10.1021/ac900468j

[96] Lauber MA, Yu YQ, Brousmiche DW, Hua Z, Koza SM, Magnelli P, Guthrie E, Taron CH, Fountain KJ. Rapid preparation of released N-glycans for HILIC analysis using a labeling reagent that facilitates sensitive fluorescence and ESI-MS detection. Analytical Chemistry. 2015;**87**(10):5401-5409. DOI: 10.1021/acs.analchem.5b00758

[108] Wakankar AA, Liu J, Vandervelde D, Wang YJ, Shire SJ, Borchardt RT. The effect of cosolutes on the isomerization of aspartic acid residues and conformational stability in a monoclonal antibody. Journal of Pharmaceutical Sciences. 2007;**96**(7):1708-1718. DOI:

Characterization of Biopharmaceuticals Focusing on Antibody Therapeutics

http://dx.doi.org/10.5772/intechopen.79107

43

[109] Rehder DS, Chelius D, McAuley A, Dillon TM, Xiao G, Crouse-Zeineddini J, Vardanyan L, Perico N, Mukku V, Brems DN, Matsumura M, Bondarenko PV. Isomerization of a single aspartyl residue of anti-epidermal growth factor receptor immunoglobulin gamma2 antibody highlights the role avidity plays in antibody activity. Biochemistry.

[110] Yan Y, Wei H, Fu Y, Jusuf S, Zeng M, Ludwig R, Krystek SR Jr, Chen G, Tao L, Das TK. Isomerization and oxidation in the complementarity-determining regions of a monoclonal antibody: A study of the modification-structure-function correlations by hydrogen-deuterium exchange mass spectrometry. Analytical Chemistry.

[111] Harris RJ, Kabakoff B, Macchi FD, Shen FJ, Kwong M, Andya JD, Shire SJ, Bjork N, Totpal K, Chen AB. Identification of multiple sources of charge heterogeneity in a recombinant antibody. Journal of Chromatography. B, Biomedical Sciences and

[112] Chu GC, Chelius D, Xiao G, Khor HK, Coulibaly S, Bondarenko PV. Accumulation of succinimide in a recombinant monoclonal antibody in mildly acidic buffers under elevated temperatures. Pharmaceutical Research. 2007;**24**(6):1145-1156. DOI: 10.1007/

[113] Kroon DJ, Baldwin-Ferro A, Lalan P. Identification of sites of degradation in a therapeutic monoclonal antibody by peptide mapping. Pharmaceutical Research.

[114] Chelius D, Rehder DS, Bondarenko PV. Identification and characterization of deamidation sites in the conserved regions of human immunoglobulin gamma antibodies.

[115] Dick LW Jr, Qiu D, Mahon D, Adamo M, Cheng KC. C-terminal lysine variants in fully human monoclonal antibodies: Investigation of test methods and possible causes. Biotechnology and Bioengineering. 2008;**100**(6):1132-1143. DOI: 10.1002/bit.21855 [116] Luo J, Zhang J, Ren D, Tsai WL, Li F, Amanullah A, Hudson T. Probing of C-terminal lysine variation in a recombinant monoclonal antibody production using Chinese hamster ovary cells with chemically defined media. Biotechnology and Bioengineering.

[117] Jiang G, Yu C, Yadav DB, Hu Z, Amurao A, Duenas E, Wong M, Iverson M, Zheng K, Lam X, Chen J, Vega R, Ulufatu S, Leddy C, Davis H, Shen A, Wong PY, Harris R, Wang YJ, Li D. Evaluation of heavy-chain C-terminal deletion on product quality and pharmacokinetics of monoclonal antibodies. Journal of Pharmaceutical Sciences.

Analytical Chemistry. 2005;**77**(18):6004-6011. DOI: 10.1021/ac050672d

2012;**109**(9):2306-2315. DOI: 10.1002/bit.24510

2016;**105**(7):2066-2072. DOI: 10.1016/j.xphs.2016.04.027

Applications. 2001;**752**(2):233-245. DOI: 10.1016/S0378-4347(00)00548-X

10.1002/jps.20823

s11095-007-9241-4

1992;**9**(11):1386-1393

2008;**47**(8):2518-2530. DOI: 10.1021/bi7018223

2016;**88**(4):2041-2050. DOI: 10.1021/acs.analchem.5b02800

[108] Wakankar AA, Liu J, Vandervelde D, Wang YJ, Shire SJ, Borchardt RT. The effect of cosolutes on the isomerization of aspartic acid residues and conformational stability in a monoclonal antibody. Journal of Pharmaceutical Sciences. 2007;**96**(7):1708-1718. DOI: 10.1002/jps.20823

[96] Lauber MA, Yu YQ, Brousmiche DW, Hua Z, Koza SM, Magnelli P, Guthrie E, Taron CH, Fountain KJ. Rapid preparation of released N-glycans for HILIC analysis using a labeling reagent that facilitates sensitive fluorescence and ESI-MS detection. Analytical

[97] Reusch D, Haberger M, Falck D, Peter B, Maier B, Gassner J, Hook M, Wagner K, Bonnington L, Bulau P, Wuhrer M. Comparison of methods for the analysis of therapeutic immunoglobulin G Fc-glycosylation profiles-part 2: Mass spectrometric meth-

[98] Terashima I, Koga A, Nagai H. Identification of deamidation and isomerization sites on pharmaceutical recombinant antibody using H(2)(18)O. Analytical Biochemistry.

[99] Liu H, Gaza-Bulseco G, Chumsae C. Glutamine deamidation of a recombinant monoclonal antibody. Rapid Communications in Mass Spectrometry. 2008;**22**(24):4081-4088.

[100] Yang H, Zubarev RA. Mass spectrometric analysis of asparagine deamidation and aspartate isomerization in polypeptides. Electrophoresis. 2010;**31**(11):1764-1772. DOI:

[101] Robinson AB, Scotchler JW, McKerrow JH. Rates of nonenzymatic deamidation of glutaminyl and asparaginyl residues in pentapeptides. Journal of the American Chemical

[102] Scotchler JW, Robinson AB. Deamidation of glutaminyl residues: Dependence on pH, temperature, and ionic strength. Analytical Biochemistry. 1974;**59**(1):319-322. DOI:

[103] Tyler-Cross R, Schirch V. Effects of amino acid sequence, buffers, and ionic strength on the rate and mechanism of deamidation of asparagine residues in small peptides. The Journal of Biological Chemistry. 1991 Nov 25;**266**(33):22549-22556. Available from:

[104] Geiger T, Clarke S. Deamidation, isomerization, and racemization at asparaginyl and aspartyl residues in peptides. Succinimide-linked reactions that contribute to protein degradation. The Journal of Biological Chemistry. 1987;**262**(2):785-794. Available from:

[105] Wakankar AA, Borchardt RT, Eigenbrot C, Shia S, Wang YJ, Shire SJ, Liu JL. Aspartate isomerization in the complementarity-determining regions of two closely related monoclonal antibodies. Biochemistry. 2007;**46**(6):1534-1544. DOI: 10.1021/bi061500t [106] Sreedhara A, Cordoba A, Zhu Q, Kwong J, Liu J. Characterization of the isomerization products of aspartate residues at two different sites in a monoclonal antibody.

Pharmaceutical Research. 2012;**29**(1):187-197. DOI: 10.1007/s11095-011-0534-2

[107] Liu H, Gaza-Bulseco G, Faldu D, Chumsae C, Sun J. Heterogeneity of monoclonal antibodies. Journal of Pharmaceutical Sciences. 2008;**97**(7):2426-2447. DOI: 10.1002/

Chemistry. 2015;**87**(10):5401-5409. DOI: 10.1021/acs.analchem.5b00758

ods. MAbs. 2015;**7**(4):732-742. DOI: 10.1080/19420862.2015.1045173

2007;**368**(1):49-60. DOI: 10.1016/j.ab.2007.05.012

http://www.jbc.org/content/266/33/22549.full.pdf

http://www.jbc.org/content/262/2/785.full.pdf

DOI: 10.1002/rcm.3831

42 Biopharmaceuticals

10.1002/elps.201000027

Society. 1973;**95**(24):8156-8159

10.1016/0003-2697(74)90040-2

jps.21180

[118] Antes B, Amon S, Rizzi A, Wiederkum S, Kainer M, Szolar O, Fido M, Kircheis R, Nechansky A. Analysis of lysine clipping of a humanized Lewis-Y specific IgG antibody and its relation to Fc-mediated effector function. Journal of Chromatography. B, Analytical Technologies in the Biomedical and Life Sciences. 2007;**852**(1-2):250-256. DOI: 10.1016/j.jchromb.2007.01.024

**Section 3**

**Development**

**Section 3**


*Article* **A Near Real-Time Method for Forest Change Detection Based on a Structural Time Series Model and the Kalman Filter**

## **Martin Puhm \*, Janik Deutscher, Manuela Hirschmugl, Andreas Wimmer, Ursula Schmitt and Mathias Schardt**

Joanneum Research, Institute for Information and Communication Technologies, Steyrergasse 17, 8010 Graz, Austria<EMAIL_ADDRESS>(J.D<EMAIL_ADDRESS>(M.H.);

<EMAIL_ADDRESS>(A.W<EMAIL_ADDRESS>(U.S<EMAIL_ADDRESS>(M.S.) **\*** Correspondence<EMAIL_ADDRESS>Tel.: +43-316-876-5107

Received: 16 July 2020; Accepted: 22 September 2020; Published: 24 September 2020

**Abstract:** The increasing availability of dense time series of earth observation data has incited a growing interest in time series analysis for vegetation monitoring and change detection. Vegetation monitoring algorithms need to deal with several time series characteristics such as seasonality, irregular sampling intervals, and signal artefacts. While common algorithms based on deterministic harmonic regression models account for intra-annual seasonality, inter-annual variations of the seasonal pattern related to shifts in vegetation phenology due to different temperature and rainfall are usually not accounted for. We propose a transition to stochastic modelling and present a near real-time change detection method that combines a structural time series model with the Kalman filter. The model continuously adapts to new observations and allows to better separate phenology-related deviations from vegetation anomalies or land cover changes. The method is tested in a forest change detection application aiming at the assessment of damages caused by storm events and insect calamities. Forest changes are detected based on the cumulative sum control chart (CUSUM) which is used to decide if new observations deviate from model-based forecasts. The performance is evaluated in two test sites, one in Malawi (dry tropical forest) and one in Austria (temperate deciduous, coniferous and mixed forests) based on Sentinel-2 time series. Both forest areas are characterized by a distinct, but temporally varying leaf-off season. The presented change detection method shows overall accuracies above 99%, users' accuracies of 76.8% to 88.6%, and producers' accuracies of 68.2% to 80.4% for the forest change stratum (minimum mapping unit: 0.1 ha). Results are based on visually interpreted points derived by stratified random sampling. A further analysis revealed that increasing the time series density by merging data from two Sentinel-2 orbits yields better forest change detection accuracies in comparison to using data from one orbit only. The resulting increase in users' accuracy amounts to 7.6%. The presented method is capable of near real-time processing and could be used for a variety of automated forest monitoring applications.

**Keywords:** state space models; forest disturbance mapping; near real-time monitoring; Sentinel-2; CUSUM

#### **1. Introduction**

Current Earth observation (EO) missions employing optical sensors such as Sentinel-2 acquire a vast volume of data: a new image every 5 days of almost every place on earth. By taking orbit overlaps into account, the time between consecutive images of the same region is reduced even further and the chance of acquiring cloud-free observations is further increased. Through high-quality georeferencing and atmospheric correction of the satellite images, it is possible to create consistent time series of measured reflectance values for any given spectral band. The vast availability of high-resolution optical data allows—for the first time—to also map small changes in near real-time. Here, "small" may apply to both spatial extent as well as spectral change magnitude. However, dense time series of high-resolution optical data have a number of characteristics that pose a challenge to change detection applications. In addition to noise effects remaining after atmospheric correction and uncertainties in the geometric registration, these challenging characteristics include:

Algorithms can be divided based on how they deal with the time series characteristics described above. The different approaches used to handle these characteristics strongly affect the algorithms' suitability for monitoring changes in near real-time. Approaches that do not account for seasonality form a first group of algorithms. A review of change detection studies using Landsat time series concludes that many older studies focused on mapping changes only at annual or biannual time scales based on series of cloud-free composite images, which always represent the same season [2]. In this context, both the Vegetation Change Tracker (VCT) [3] as well as the LandTrendr approach [4] represent widely used algorithms, but they are not designed for near real-time mapping.

A second group of algorithms explicitly accounts for seasonality by using regression models based on trigonometric functions to capture the intra-annual variations (variations within one year, i.e., seasonality) of the spectral signatures independently for each pixel. With this approach, also frequently referred to as harmonic regression, periods of stable land cover are modelled as a deterministic, continuous function of time. Irregular sampling intervals and data gaps are therefore not a problematic issue, but the deterministic nature of the model does not allow inter-annual variations of the seasonal pattern (variations between different years, e.g., shifts in seasonality). The seasonal model represents an average of different conditions occurring within a stable period, e.g., dry and wet years, late and early leave outbreak. A widely used algorithm belonging to this group is the Breaks for Additive Season and Trend (BFAST) algorithm [5] and its evolution BFAST Monitor [6]. While the latter is tailored to near real-time mapping of new changes, the original version is intended for the analysis of historic time series. Both versions have been used in a variety of studies and can be applied to detect both abrupt and gradual changes. Concerning the robustness of BFAST to invalid observations, it has been stated that occasional signal artefacts are well handled, but temporally aggregated occurrences such as several consecutively un-masked clouds can be a source of error [7]. Also, additional pre-processing to eliminate artefacts was applied [8–10]. The second widely used implementation of the harmonic regression approach is the Continuous Change Detection and Classification algorithm (CCDC), where the original concept [11] is extended to include more types of land cover besides forest, as well as a classification framework [12]. From the beginning, CCDC was designed to work with dense Landsat time series and can handle seasonality, irregularly spaced observations, and signal artefacts to some extent. Both abrupt and gradual changes can be detected. Some further updates to the algorithm include (i) a mechanism to automatically adjust the complexity of the time series model based on the number of available clear observations, as well as (ii) a different method to estimate the model parameters which reduces overfitting [13]. A third algorithm employing harmonic regression

utilizes residuals from the regression together with statistical quality control charts [14]. This approach comprises signal artefact detection with Shewhart X-bar charts. After the elimination of artefacts, both abrupt and gradual changes are indicated by exponentially weighted moving average (EWMA) charts in a near real-time manner. All of the described algorithms of the second group share certain basic concepts, but the individual implementations vary. They are designed to process large amounts of data in a highly automated way and therefore rely on data-driven statistical boundaries for detecting change, although the distinct nature and computation of these boundaries is quite different.

The third group of algorithms stands out by also taking inter-annual variations of the seasonal pattern into account. Structural time series models are set up in terms of components, such as trends and cycles, which have a direct interpretation. Their statistical treatment is based on the state space form and the Kalman filter and first described for time series analysis in econometrics [15]. Compared to the harmonic regression approach, the model is no longer deterministic. Kalman filtering denotes a versatile parameter estimation technique which yields optimal estimates in a statistical sense [16]. It is well established in many application fields and has been applied to numerous signal tracking problems [17]. The combination of structural time series models with the Kalman filter and the concept's suitability for remote sensing purposes has been investigated in a proof of concept study, where Landsat time series are used to detect storm damages in a small forest test site in Germany [18]. The Kalman filter has also been applied to time series of MODIS 8-day composites in order to detect insect-induced defoliation in near real-time at a forest test area in northern Sweden [19]. This study makes use of the CUSUM control chart [20] to indicate changes, but it does not combine structural time series models with the Kalman filter. Instead, the filter is used to derive a smoothed time series of the Normalized Difference Vegetation Index (NDVI) based on a global model trajectory.

With the advantages and limitations of existing algorithms in mind, this work combines the pixel-by-pixel modelling typical for existing harmonic regression algorithms with the Kalman filter's capability to dynamically adjust the model based on new observations. The main aim of this study is to present an innovative change detection algorithm for optical EO data which is based on a structural time series model and the Kalman filter. It is largely data-driven and designed especially for near real-time mapping in web- or cloud-based monitoring services. The algorithm presented in this paper accounts for seasonality and also allows inter-annual variations of the seasonal pattern, e.g., vegetation phenology. Furthermore, strategies for handling irregular sampling intervals and signal artefacts are presented. The algorithm is tested in a forest change detection application using time series of Sentinel-2 data (S-2). Forest disturbances are detected at two complex forest test sites in Austria and Malawi. The first test site is an alpine area in Austria characterized by frequent cloud cover, snow cover, strong topographic effects, and pronounced forest seasonality and phenology. The second test site is located in the dry tropical forests of Malawi, where forests show a strong and varying seasonality between dry and rainy seasons. In the Malawi test site, we also analyze and compare accuracy results for two different data scenarios: first, using S-2 images from only one orbit, and second, using all available S-2 imagery from two orbits. The aim of this analysis is to investigate if different viewing angle and inconsistent geo-location of the pixels resulting from the combination of two orbits decrease the overall change detection accuracies despite the boosted time series density.

#### **2. Materials and Methods**

#### *2.1. Test Sites, Data, and Pre-Processing*

For the forest change detection demonstration, we selected two test sites, one in Austria and one in Malawi. The location of the two sites is shown in Figure 1. The Austrian test site is located in the south-eastern part of the country. The test site is characterized by strong topography in the northern part (Alpine area), where coniferous forests dominate. The southern part of the test site is located in the foothills of the Alps with moderate topography and the forests are predominantly mixed forests composed of coniferous and deciduous trees. The annual temperature amplitude in the Alpine area is stronger than in the foothills; however, the Alpine coniferous forests have less pronounced phenological dynamics in comparison to the mixed and deciduous forests that dominate in the foothills. The eastern half of the study area is covered by two orbits, and the western half only by one orbit (see Figure 1). This is a typical data scenario encountered in practical applications.

**Figure 1.** Location of the test sites.

The Malawi test site is characterized by flat to slightly hilly topography. The land cover is very heterogeneous and is subject to a precipitation gradient from east to west. As a result, the area is characterized by a vegetation phenology gradient from east to west. The forested areas differ in tree-cover density and tree-type composition and therefore, show very different spectral behavior. Dry tropical forests and the surrounding land use classes are more difficult to classify and monitor than humid evergreen forests, as they show a typical phenological development from highly vital in the rainy season to dry and leafless in the dry season [21]. Understory and grassland fires beneath the forest canopy can further complicate forest classification. Two S-2 orbits (relative orbits 92 and 135) cover the Malawi study area. The size of the test site was clipped to the overlap area of two orbits to investigate the effect of separately processing data from one and two orbits.

Both test sites are located in areas that show distinct seasonal patterns due to phenology. Austria has the typical European summer growing season with leaf-off time in winter for deciduous species due to low temperatures. The forests of Malawi also show strong phenological variation as water scarcity during the dry season (typically May–October) causes leaf-fall. The typical temporal NDVI signatures of different forest types are shown in Figure 2. Each time series represents two years of NDVI observations for a single pixel corresponding to a specific forest type and test site. In all cases, data from two orbits is used and cloud/snow masking has been applied as described below. Aside from the different seasonal patterns, several other characteristic properties of the data can be observed:

**Figure 2.** Typical temporal NDVI (Normalized Difference Vegetation Index) signatures over two years for single pixels corresponding to specific forest types and test sites. Larger data gaps occur in the winter months due to snow cover in Austria (AUT) or the rainy season in Malawi (MWI). Outside of these periods, the observation density is high enough to capture the seasonal patterns.

The current data quality report issued by the Sentinel-2 Mission Performance Centre (S2 MPC) gives statistics for the multi-temporal geometric registration performance [22]. For about half of the images, the co-registration error is larger than 0.5 pixels at 10 m resolution. Within homogenous land-cover areas, this can be treated as an additional noise component. The high geometric uncertainty becomes a larger problem at the border regions between land-cover classes, especially if a given pixel jumps between forested and non-forested states.

All available Sentinel-2 scenes with a nominal cloud cover below 90% were downloaded for the test sites (Table 1) at Level-1C and atmospherically corrected to surface reflectance (SR) values using the Sen2Cor processor version 2.5.5 [23]. We then resample the 20 m bands to 10 m spatial resolution and stack all bands to a 10-band output image. We calculate a combined cloud, cloud shadow, and snow mask with the FMask algorithm [24,25]. This mask is slightly altered by morphological operations (erode, expand) to fill cloud holes and the masked pixels are then removed from the pre-processed S-2 imagery by assigning no-data values to them. We also perform a topographic correction based on a modified Minnaert correction [26] using the Shuttle Radar Topography Mission (SRTM) model at 30 m spatial resolution as digital elevation model (DEM).

**Table 1.** Earth ObservationData Information.

#### *2.2. Change Detection Method Using a Structural Time Series Model and the Kalman Filter*

The underlying assumption of the monitoring approach presented in this section is that the normal temporal trajectory of a given spectral band over the course of the year can be captured by a univariate structural time series model. Within the structural model, we can further distinguish between an observation sub-model and a dynamic sub-model. The observation sub-model on the one hand defines the relationship of the measurements to a set of state variables which cannot be observed directly. In a structural time series model, the state variables usually represent the series' additive decomposition into trend, seasonal, and long-term cyclical components. The dynamic sub-model on the other hand describes the expected temporal evolution of the state variables. By formulating the dynamic sub-model in continuous time, the problem of irregular sampling intervals is addressed. The Kalman filter is used to fit the model to the data and operates recursively from one point in time to the next. Each recursion may be divided into two steps. In the time-update step, the states' temporal evolution is predicted based on the dynamic sub-model. It is followed by the measurement update step, where the predicted state estimate is enhanced by incorporating newly available observations. Abrupt changes of the spectral signature, possibly linked to a forest disturbance, are indicated by statistically significant deviations between new observations and the Kalman filter predictions. The recursive operation of the filter further implies that prior knowledge about the initial state is required. Because the ability to distinguish anomalies from normal seasonal changes depends on the quality of the state estimates, a proper initialization is of high importance. Therefore, a robust least squares method is used to estimate the initial state from a historic time series covering at least one full year prior to the beginning of the monitoring period (see initialization time window in Table 1). The following sub-sections describe the implementation in more detail.

#### 2.2.1. Time Series Model

Structural time series models are mathematically formulated using the discrete-time state space representation [15]. This concept assumes that a linear, time-variant system can be described by a set of state variables. Because these variables can usually not be observed directly, an observation sub-model linking the system state to a set of measurements is required. In case of univariate time series, the measurement equation is

$$z\_k = \mathbf{h}\_k \mathbf{x}\_k + r\_{k\prime} \tag{1}$$

where **h***k* is a row vector, **x***k* denotes the state vector, and *zk* is a scalar observation made at time *tk*. In addition, *rk* represents serially uncorrelated observation noise with mean zero and variance *Rk*. We will refer to discrete points in time as epochs and define an epoch index *k* ∈ N to indicate a possible time-dependency of variables. The temporal evolution of the state vector is described by a dynamic sub-model using the transition equation:

$$\mathbf{x}\_{k} = \Phi\_{k}\mathbf{x}\_{k-1} + \mathbf{q}\_{k'} \tag{2}$$

where **Φ***k* denotes the transition matrix and **q***k* is a vector of serially uncorrelated process noise with mean zero and covariance matrix **Q***k*. Note that no assumptions regarding the distributions of the observation and process noises are made at this point, but they are supposed to be uncorrelated with each other in all epochs. It is further assumed that the initial state **x**0 is known with a level of uncertainty characterized by the state error covariance matrix **P**0.

The structural time series model implemented in this study enables a decomposition of the time series into additive trend, seasonal, and irregular components. We assume a constant trend and directly introduce the state variable μk, which represents the trend level at time *tk*. The seasonal component on the other hand is constructed by a sum of *P*-periodic cosine waves, with *P* denoting the fundamental duration of the seasonal cycle. Considering the nature of the time series investigated here, it is appropriate to measure time in days and thus set *P* to 365.25. For each cosine wave added to the seasonal component, two more variables are appended to the state vector. Taking a wave with frequency ω*i*, the state variable γ*i,k* representing the cosine waves' level at time *tk* is added along with the variable γ∗ *i*,*k* , whose interpretation is not particularly important. The structural model for the Malawi test site features a seasonal component with frequencies corresponding to one and two periods per year. In case of the Austrian test site, we use just the fundamental frequency because of the less complex seasonal pattern. The measurement Equation (1) of the Malawi model becomes

$$z\_k = \begin{bmatrix} 1 & 1 & 0 & 1 & 0 \end{bmatrix} \begin{bmatrix} \mu \\ \nu\_1 \\ \nu\_1^\* \\ \nu\_2^\* \\ \nu\_2^\* \\ \nu\_2^\* \end{bmatrix}\_k + r\_k \tag{3}$$

Note that *zk* represents the observed reflectance value at time *tk* and **h** is time-invariant. Using Δ*tk* = *tk* −*tk*−1 measured in days and ω*i* = 2π*i*/365.25, the transition Equation (2) of the Malawi model becomes

$$
\begin{bmatrix}
\mu\\\nu\_1\\\nu\_1'\\\nu\_2\\\nu\_k'\end{bmatrix}\_k = \begin{bmatrix}
1 & 0 & 0 & 0 & 0\\0 & \cos\left(\omega\_1 \Delta t\_k\right) & \sin\left(\omega\_1 \Delta t\_k\right) & 0 & 0\\0 & -\sin\left(\omega\_1 \Delta t\_k\right) & \cos\left(\omega\_1 \Delta t\_k\right) & 0 & 0\\0 & 0 & 0 & \cos\left(\omega\_2 \Delta t\_k\right) & \sin\left(\omega\_2 \Delta t\_k\right)\\0 & 0 & 0 & -\sin\left(\omega\_2 \Delta t\_k\right) & \cos\left(\omega\_2 \Delta t\_k\right)
\end{bmatrix}
\begin{bmatrix}
\mu\\\gamma\_1\\\gamma\_2'\\\gamma\_2'\\\gamma\_2'\\\gamma\_2'
\end{bmatrix}\_{k-1} + \begin{bmatrix}
q\_{\mu}\\q\_{\nu\_1}\\q\_{\nu\_2}\\q\_{\nu\_2}\\q\_{\nu\_2}\\\end{bmatrix} \tag{4}
$$

In order to fully specify the state space model, two more quantities need to be defined: *Rk* and **Q***k*. The first one will be addressed in the next sub-section. The covariance matrix **Q***k* quantifies the uncertainties within the dynamic sub-model and is defined depending on the nature of the process noise. Here, we assume that the state transition is affected by multivariate, continuous time white noise with constant variances *Qc* independently specified for the trend and seasonal component, that is:

$$\mathbf{Q}\_{k} = \Delta t\_{k} \begin{bmatrix} Q\_{c}^{\text{tred}} & 0 & 0 & 0 & 0 \\ 0 & Q\_{c}^{\text{eass}} & 0 & 0 & 0 \\ 0 & 0 & Q\_{c}^{\text{eass}} & 0 & 0 \\ 0 & 0 & 0 & Q\_{c}^{\text{eass}} & 0 \\ 0 & 0 & 0 & 0 & Q\_{c}^{\text{eass}} \end{bmatrix} \tag{5}$$

Note that all off-diagonal elements of **Q***k* equal zero and we therefore assume that the different process noise components are uncorrelated. Implementing the dynamic model in continuous time means that the time series model is specified for arbitrary values of Δ*tk* and irregular sampling intervals do not pose a problem. Finding appropriate values for *Qc* is one of the difficulties. Our experiments showed that setting individual values for each pixel and band proportional to the respective observation variance *R* works quite well. Since *R* is estimated from the data (see the next sub-section), the user has to specify two proportionality factors, one for the trend component and one for the seasonal component. We recommend to let *Q*seas *c* be larger than *Q*trend *c* in order to make the time series model more responsive to phenological variations.

#### 2.2.2. Initial State Estimation

In order to obtain estimates of the initial state vector **x**0, its associated covariance matrix **P**0, and the observation variance *R*, the state space model outlined above is transformed to a linear regression model with the measurement equation:

$$
\begin{bmatrix} z\_1 \\ z\_2 \\ \vdots \\ z\_m \end{bmatrix} = \begin{bmatrix} \mathbf{h}\Phi(t\_1, t\_0) \\ \mathbf{h}\Phi(t\_2, t\_0) \\ \vdots \\ \mathbf{h}\Phi(t\_m, t\_0) \end{bmatrix} \mathbf{x}\_0 + \begin{bmatrix} r\_1 \\ r\_2 \\ \vdots \\ r\_m \end{bmatrix} \tag{6}
$$

Equation (6) also illustrates the key difference between a structural time series model and a regression model. The latter does not distinguish between observation sub-model and dynamic submodel and has no notion of process noise. Standard least squares methods may be used to obtain estimates of **x**0 and **P**0 from a batch of *m* historic observations acquired before the monitoring period. We applied a robust parameter estimation approach known as iteratively reweighted least squares (IRLS). The technique belongs to the class of *M*-estimators [27] and the implementation follows [28]. We further use the mean squared error of the weighted least squares fit as an estimate for the observation variance *R*. If the number of valid historic observations is low, for example in areas with extremely high cloud probability, the observation variance can be severely underestimated. Existing studies applying harmonic regression also report this issue and we take up the suggestion to define a certain minimum value for *R* that is used if the estimate is lower [12,29].

#### 2.2.3. Kalman Filter

Once the time series model and the initial values **x**0 and **P**0 are defined, the discrete-time Kalman filter recursion can be used to obtain estimates for the state and its error covariance matrix in subsequent epochs. The time update step yields the predicted (a-priori) estimates **~ x***k* and **~ P***k* based on the dynamic model and the previous estimates at time *tk* − 1:

$$
\tilde{\mathbf{x}}\_k = \boldsymbol{\Phi}\_k \hat{\mathbf{x}}\_{k-1} \tag{7}
$$

$$
\tilde{\mathbf{P}}\_k = \boldsymbol{\Phi}\_k \hat{\mathbf{P}}\_{k-1} \boldsymbol{\Phi}\_k^\top + \mathbf{Q}\_k \tag{8}
$$

Then, the a-priori measurement residual *yk* and its variance *Ck* is computed according to Equations (9) and (10). The residual represents the difference of the prediction to the actual measurement and is referred to as innovation, since it contains new information currently not present in the predicted state [17].

$$y\_k = z\_k - \mathbf{h}\tilde{\mathbf{x}}\_k\tag{9}$$

$$\mathbf{C}\_{k} = \mathbf{h}\tilde{\mathbf{P}}\_{k}\mathbf{h}^{\mathrm{T}} + \mathbf{R} \tag{10}$$

In the measurement update step of each recursion, new information is merged with the predictions to obtain improved (a-posteriori) estimates **^ x***k* and **^ P***k*. The Kalman gain **k***k* (11) determines how much the newly acquired measurement will influence the a-posteriori estimates of the state and its error covariance and appears in both update Equations (12) and (13), where **I** represents an identity matrix.

$$\mathbf{k}\_k = \mathbf{\tilde{P}}\_k \mathbf{h}^T \mathbf{C}\_k^{-1} \tag{11}$$

$$
\hat{\mathbf{x}}\_k = \tilde{\mathbf{x}}\_k + \mathbf{k}\_k y\_k \tag{12}
$$

$$
\hat{\mathbf{P}}\_k = (\mathbf{I} - \mathbf{k}\_k \mathbf{h}\_k) \tilde{\mathbf{P}}\_k \tag{13}
$$

#### 2.2.4. Signal Artefact Handling and Change Detection

When the sequence of measurements processed by the filter may contain artefacts, an additional anomaly detection step should be included before the measurement update. The properties of the innovations can be exploited to detect anomalous measurements by means of a statistical test. Provided that the underlying model assumptions are valid, and the observation noise is Gaussian, the innovations should be normally distributed with mean zero and variance *Ck*. The test statistic *T*ˆ *k* given in (14) follows the χ2-distribution with a single degree of freedom. The hypotheses to be tested on a significance level α are stated in (15) and (16), respectively:

$$
\mathcal{T}\_k = \frac{y\_k^2}{\mathcal{C}\_k}, \quad \mathcal{T}\_k \sim \chi^2 \tag{14}
$$

$$H\_0: \ y\_k = 0 \text{ if } \hat{T}\_k \le \chi^2\_{1, \ 1-\alpha} \tag{15}$$

$$H\_1: \ y\_k \neq 0 \text{ otherwise} \tag{16}$$

Considering that anomalous observations will cause large innovations, the null hypothesis will be rejected. In that case, the measurement update is not carried out so that the state estimate remains unbiased.

The change detection part of the algorithm is based on the cumulative sum control chart (or CUSUM) [20]. The key quantity here is again the sequence of filter innovations, which should have zero mean if the Kalman filter model assumptions reflect the truth closely enough. We use the CUSUM control chart to monitor if the innovations deviate from mean zero. Of course, the presence of artefacts like un-masked clouds presents a problem in this regard, because a single innovation corresponding to an artefact may shift the mean significantly and hence trigger a (inaccurate) change signal. The previously described anomaly test cannot distinguish between an artefact and an abrupt land cover change. To work around this limitation, a new quantity we call edited innovation *y*ˇ*k* is introduced. It represents the original innovation divided by its standard deviation, but also limited in magnitude based on the significance level α specified for the anomaly test, that is:

$$\vec{y}\_k = \begin{cases} \min\left(\frac{y\_k}{\sqrt{\zeta\_k}}, \sqrt{\chi\_{1,1-\alpha}^2}\right) \text{if } y\_k \ge 0\\ \max\left(\frac{y\_k}{\sqrt{\zeta\_k}}, -\sqrt{\chi\_{1,1-\alpha}^2}\right) \text{if } y\_k < 0 \end{cases} . \tag{17}$$

Dividing by the standard deviation means that sequences of edited innovations should be close to having unit variance across different bands and pixels. The limit operation on the other hand ensures that the CUSUM control chart with *y*ˇ*k* as input is less sensitive to single statistical outliers. A temporal aggregation of edited innovations with the same sign is required to shift the mean of the sequence significantly. The CUSUM test statistic for a positive mean shift with respect to spectral band *i* is implemented as:

$$\begin{array}{ll} S\_i^+ (0) &=& 0\\ S\_i^+ (k) &=& \max(0, S\_i^+(k-1) + \nexists j\_i(k) - d) \end{array} \tag{18}$$

where *d* is a drift parameter which generally compensates small deviations and also ensures that effects of occasional signal artefacts on the test statistic fade away over time. A mean shift is signaled if the test statistic crosses a predefined threshold. Instead of evaluating all processed bands separately, we decided to aggregate the test statistics of several bands by simple summation and then use a global threshold. We look for anomalous reflectance increases in the red, red edge, and short wave infrared

(SWIR) bands, as these have proven to be sensitive to vegetation changes [30]. The implemented criterion for detecting a forest disturbance at time *tk* is given in Equation (19) using S-2 band numbers.

$$\sum\_{i} S\_i^+(k) > \text{change threshold where } i \in \{\text{B04, B05, B11, B12}\} \tag{19}$$

Both the change threshold in (19) and the drift parameter in (18) are user-defined tuning parameters. Appropriate values need to be determined empirically. Because of the limit operation described by (17), the maximum of a single edited innovation is a known constant and amounts to ~2.57 if α = 1%. This knowledge provides a helpful yardstick for setting both threshold and drift parameters. The statistical normalization applied to the edited innovations means that the same values can be used globally for all pixels. The flowchart depicted in Figure 3 illustrates how the methods discussed in the preceding sections are joined together in order to create a data-driven algorithm capable of detecting abrupt changes on the pixel level. A summary and some additional explanatory comments are given below.

**Figure 3.** Flow chart of the Kalman filter approach.

Figure 4 visualizes the method for a forest change pixel (10 by 10 m) in a deciduous forest in the Austrian test site. The first four sub-plots show the surface reflectance of Sentinel-2 bands 4 (red), 5 (red edge), 11 (SWIR1), and 12 (SWIR2) over time. The sub-plot at the bottom shows the value of the CUSUM test statistic defined in Section 2.2.4 over time. A phenological cycle that is typical for deciduous tree species can be observed. Higher reflectance values correspond to the leaf-off season during winter. The red line represents the Kalman fit and the grey area corresponds to the 90% confidence interval of the model forecast. While blue observations are considered in the model and measurement update step (compare Figure 3), the orange observations are flagged as anomalous and therefore ignored during the measurement update step. Please note that the level of significance for the anomaly test is α = 1%, hence blue observations may also appear outside of the plotted confidence interval. Occasional positive anomalies will cause a short-lived increase of the related CUSUM, but in most cases, it should not exceed the threshold (e.g., middle of 2016 or end of 2017 in Figure 4). Such anomalies occur if for example small cloud artefacts remain in the pre-processed data. A prolonged increase of the CUSUM is triggered by a persistent signal shift. At some point, the threshold is exceeded, and the pixel is flagged as changed, in this case, at the end of 2018. The vertical dashed line in Figure 4 marks the date on which the change is visible for the first time in the S-2 imagery. Note that anomalies with significant magnitude are first detected only in the SWIR bands, while the signal shift in the red and red edge band initially occurs within the 90% confidence interval; however, the aggregated CUSUM test statistic allows a timely detection of the change. After a threshold crossing, the test statistic is reset to zero. Repeated change alerts, as shown in this example, may thus occur. This could be used to increase the confidence about a detected change at the cost of delaying the detection; however, this aspect has not been investigated in detail in this work. For the resulting forest change maps at the two test sites, we first include all Sentinel-2 single pixel changes (10 by 10 m) that were detected by the described approach within the change detection time window. For the final forest change maps, a minimum mapping unit of 0.1 ha is applied to the forest change stratum, which relates to change areas represented by at least ten connected 10 m pixels. Detected forest change areas smaller than 10 pixels are removed from the final forest change maps.

**Figure 4.** Example of the time series model for a deciduous forest pixel at the Austrian test site, where a forest change occurred at the end of 2018. The plots show the surface reflectance of Sentinel-2 bands 4, 5, 11, and 12. The vertical dashed line marks the date on which the change is visible for the first time in the S-2 image time series. The plot at the bottom shows the cumulative sum control chart (CUSUM) used for detecting the change event (coordinates: X 559045 m, Y 5235095 m, in EPSG 32633).

#### *2.3. Evaluation Method*

In order to test the operational application of the presented forest change detection method, a full area-based validation of the resulting change maps is performed in both test sites based on stratified random sampling points that are located within forest areas outlined by benchmark forest masks (status 2016 for Malawi and summer 2018 for Austria). The overall number of samples and samples per stratum are based on recommendations for land cover accuracy estimation [31,32]. Points that were found to be already deforested before the beginning of the change detection time window were excluded from the analysis since these cases correspond to errors of the initial forest masks.

For Malawi, 849 reference sample points are used for statistical analysis, with 735 sample points belonging to stratum forest and 114 sample points belonging to stratum change. During blind interpretation, we flagged all sample points that were located at or very close (~10 m) to the border of a forest change patch. This allows us to treat these points as correctly classified in a subsequent plausibility analysis. This approach is usually termed "plausibility analysis", because it is deemed plausible that the border point can also belong to the other stratum. This is especially true for the small-scale 10 by 10 m change assessment we use in this study considering that some of the spectral bands of Sentinel-2 used to derive the changes only have a nominal spatial resolution of 20 m. For Malawi, assessments are based on the following validation approaches and input data options:

The plausibility analysis is only performed for the two-orbit input data. Thus, the comparison of the two input data options in Malawi is based only on the blind assessment approach.

For Austria, we interpreted a total of 1585 reference points, of which 21 points were removed as they were found to be non-forest already before the beginning of the change detection window. From the remaining 1564 points, 1212 belong to stratum forest and 352 to stratum change. The plausibility approach was carried out in the same manner as in Malawi. Please note that the Austrian test site is only partly covered by two orbits, which is a typical data scenario when working with Sentinel-2 data. Combined input data from two orbits is used where possible. For all assessments, we provide unbiased estimates of the mapped area proportions as well as the products', users', and producers' accuracies by applying Equations (1) and (6)–(8) from Reference [31].

#### **3. Results**

Table 2 gives a summary of key validation results for both test sites, different assessment approaches, and input data scenarios. The detailed results of the forest disturbance detection are shown in Tables 3–7, where the upper part presents the sample counts and the lower part presents the unbiased estimates of area proportion and accuracy measures. Overall accuracies are very high (96.4–99.3%). This is not surprising since the unchanged forest class accounts for 98.7% of the validation area at the Malawi test site and for 98.6% at the Austrian test site. It is also evidence for a low rate of false-positive change detections. For better comparison, Table 2 also lists the users' accuracies and producers' accuracies of the different validation approaches for the change class. The plausibility analysis increases the users' accuracy of the change class by 17.5% in Malawi and by 4.1% in Austria. Producers' accuracies show a strong increase of 31.1% in Malawi and 29.1% in Austria. Overall accuracies after plausibility analysis reach 99.3% at both test sites. Results show that combining data from two orbits leads to an increase in users' accuracy of 7.6% (Malawi—blind validation approach) compared to using data from only one orbit, while producers' accuracies remain the same. The detailed accuracy metrics for Malawi are shown in Tables 3–5, and for Austria, they are listed in Tables 6 and 7. For comparison between one and two orbits, please compare Tables 3 and 4. For comparison of blind and plausibility results, please compare Tables 4 and 5 for Malawi and Tables 6 and 7 for Austria.

**Table 2.** Summary of the accuracy measures for the forest change maps at both test sites for different assessment approaches (blind interpretation and plausibility analysis).

Confidence interval of accuracy measures: 95%.

#### **Table 4.** Error matrix for Malawi, blind approach—2 orbits.

Confidence interval of accuracy measures: 95%.

**Table 5.** Error matrix for Malawi, plausibility approach—2 orbits.

Confidence interval of accuracy measures: 95%.

**Table 6.** Error matrix for Austria, blind approach.

Confidence interval of accuracy measures: 95%.

**Table 7.** Error matrix for Austria, plausibility approach.

Confidence interval of accuracy measures: 95%.

Figure 5 shows some mapping examples of windthrow detection for the Austrian test site. The series of S-2 images illustrates the development of forest disturbances occurring at an alpine subset of the test site that was affected by windthrow late in the year 2018 (storm Vaia on 29/30 October 2018). Please note that the depicted sequence does not represent all available imagery, but a selection of cloud-free images of the area of interest. The first image (top left) shows the undisturbed state two weeks before the storm event. A large windthrow area can be identified in the second image (top middle; surrounding the red circle). Forest change detection is complicated by the fact that the timing of the storm exactly coincides with leaf discoloring and leaf-fall for the broadleaf trees at the site and by subsequent bad weather conditions also leading to snow cover. Remaining snow (blueish colored areas) can still be seen in early April 2019 (top right image). Harvesting of damaged forest areas leads to a continuous increase in deforested areas in the subsequent images. Typically, also areas adjacent to completely thrown forest areas are affected to some degree by wind damage, e.g., single tree throws or broken stems. In Austria, all storm damage affected areas are usually harvested directly after the storm event to prevent bark beetles from spreading. The last image in the series shows which pixels are flagged as changed by the algorithm and in which month the changes were detected. Due to frequent cloud and snow cover at the site from mid-November 2018 to April 2019, only few pixels were detected as changed shortly after the storm as only two usable post-storm observations were available until the end of 2018. Many detections are thus delayed until early spring 2019, when the time series of snow-free Sentinel-2 imagery continues and when the damaged forest areas are being cleared (removal of damaged and broken trees). This example shows both the capabilities, but also the limitations of the method when used in near real-time forest change mapping scenarios for Central European/Alpine forests.

**Figure 5.** Examples of near real-time forest change detection in Austria. Storm "Vaia" windthrows on 29/30 October 2018 at an alpine forest site in Styria (~1400 m above sea level).

The software needed to create the presented mapping example and time series plots has been implemented in Python. This Python implementation is capable of processing smaller test sites (up to 10 MP) for development and demonstration purposes in reasonable time. For large-area processing, a performance-optimized C/C++ implementation is also available. However, both versions rely on certain in-house modules and libraries which are subject to licensing restrictions.

#### **4. Discussion**

When comparing our results to those of other studies on near real-time forest change detection, similarly high overall accuracy values can be observed. This is to be expected for a land cover class and test sites that are characterized by a large proportion of unchanged areas. In this case, users' and producers' accuracies of the change class are much better measures to determine the suitability and

practical implementation potential of a monitoring approach. Many recent near real-time forest change monitoring approaches detect and evaluate changes at a minimum mapping unit (MMU) of at least one Landsat 8 pixel (~0.1 ha) as most related studies are still primarily based on Landsat 8 imagery with a lower spatial resolution than Sentinel-2 [29,33–35]. For better comparability, we have also applied a 0.1 ha MMU for the change class, but our change polygons can be of any shape representing at least ten connected Sentinel-2 change pixels.

There are a number of near real-time forest change studies and algorithms that provide accuracy statistics similar to those presented here. The near real-time humid tropical forest monitoring approach of Global Forest Watch was evaluated at national scale in Peru with a reported producer's accuracy of 75.4% and a user's accuracy of 92.2% [33]. However, a direct comparison of results is difficult, since 63% of the detected forest changes in this study were at least one hectare in size, while only 4% were the size of one single Landsat pixel (~0.1 ha). In our test sites, the vast majority of change polygons is smaller than 0.5 ha. Another change detection approach has recently been developed and tested for seasonal tropical forests in Myanmar. It uses a harmonic regression model to account for seasonality and a set of time series disturbance probabilities to detect forest changes [35]. The authors report overall disturbance detection accuracies of 78.3% for Landsat 8 data and 83.6% for Landsat 8 data combined with Sentinel-1 data. The reported users' accuracy for the disturbance class at the single Landsat 8 pixel level (~0.1 ha) is 84.1% and producers' accuracy is 78.6%, but disturbance detections are significantly delayed by 65 days on average. Very high users' and producers' accuracies of 88% and 89% were also reported for a forest change detection approach that combines Landsat 8, Sentinel-1, and ALOS-2 PALSAR-2 data at a dry tropical forest site in Bolivia characterized by distinct dry and wet seasons [34]. Instead of adapting the time series model in near real-time to the observed phenology, the authors apply a-priori spatial normalization to reduce the dry forest seasonality in the time series and then apply the change detection analysis on the normalized data [10]. The method is reported to perform well for extreme events, but we would expect that such a combined normalization and change detection approach could fail in years that behave significantly different from the average, e.g., years with extreme dryness (deviation in magnitude) or a very late start (deviation in time) of the wet season.

Our approach to use the Kalman filter is quite different in that it continuously accounts for inter-annual phenology variations and updates the time series model. The method proved to work properly at both test sites even for cases where the phenology curve significantly differs between years. Figure 6 shows three false-color Sentinel-2 images (band combination B11, B04, B03) of an unchanged forest area at the Malawi test site from three consecutive years, acquired almost on the same day of each year (mid-August). Significant differences in the phenological state can be observed both in the images as well as in the corresponding time series plot below (Figure 6). The time series plot shows the Sentinel-2 red band surface reflectance for a pixel at the center of the red circle in the images. Each year shows time periods with strong deviations from the average IRLS fit over all years. Year 2016 shows an early and stronger than average dry season, while the observations of year 2018 indicate a pronounced prolongation of the spring rainy season and thus late start of the dry season. The observed reflectances for August 2016 and August 2018 deviate by more than 7% of total surface reflectance. A widely used deterministic modelling approach based on harmonic regression (see IRLS fit represented by the green line in Figure 6) results in large and prolonged deviations between observations and the model curve which is disadvantageous for change detection. The Kalman filter is able to track the differences in plant vitality much better. In the given example, only one single observation in December 2018 was flagged as anomalous.

**Figure 6.** The illustration shows 3 false-color images (B11, B04, B03) acquired in three consecutive years in Malawi, almost on the same day (gray dashed lines). The upper time series plot shows the red band surface reflectance for a single pixel located in the center of the red circle together with a harmonic regression fit (IRLS—iteratively reweighted least squares) and the Kalman filtered trajectory. The lower plot shows the trend and seasonal components separately.

This seasonal phenological effect is further illustrated in Figure 7, which compares the residuals of the IRLS fit shown in Figure 6 with the corresponding Kalman filter innovations. The IRLS residuals show a much larger degree of unwanted systematic patterns left in the sequence, which is apparent when comparing the moving mean of the last 15 observations (dashed orange line in Figure 7). These systematic patterns represent the difference in phenology (both in time and magnitude) from the "average" fit. In combination with the CUSUM test, these systematic deviations of the IRLS fit would result in a larger amount of erroneous forest change detections. In case of the Kalman filter innovations, remaining non-random patterns appear when the model receives adjustments to the current phenology (especially the second half of 2018 in Figure 7). However, the drift parameter in Equation (18) avoids that small and short-lived deviations from the zero mean assumption accumulate to a change signal.

**Figure 7.** The residuals of the regression fit (IRLS) and the corresponding Kalman filter innovations for the same plot as in Figure 6. The orange dashed line represents the moving mean based on the last 15 observations.

We could show that the combination of Kalman filter innovations and the CUSUM test is suitable for rapid near real-time detection of abrupt forest changes as shown in Figure 4. With the CUSUM test, changes can be detected with a smaller time lag than with traditional change detection methods based on multiple confirmations [21,35–37]. Many existing algorithms require a fixed number of observed anomalies to signal a change, and some enforce the condition that anomalies also have to be detected consecutively. Certain statistical boundaries or change probabilities specific to the method have to be exceeded multiple times for a change to be recorded. When the spectral footprint of a change event is at the border of detectability, these restrictions likely lead to omission errors or a large time lag between change event and detection. Because of its cumulative nature, the proposed CUSUM test statistic is not restricted in the same away. Any post-change observation can add information to the test statistic, regardless of being flagged by the anomaly test or not. Depending on the spectral footprint of the change event, it may take two, three, or more post-change observations to detect the change. At this point, we cannot give a quantitative evaluation of the typical time lag because appropriate reference data were not available.

The reduction of the time lag between change event and change detection is of major importance to near real-time forest monitoring systems. For this reason, several recent studies have combined optical data with Synthetic Aperture Radar (SAR) data in order to reduce this time lag [34,35]. Regarding optical data, persistent cloud cover during rainy seasons represents a major limitation as the time gap between consecutive valid observations may become very large. Change detection methods that require multiple confirmed change detections will therefore show large time lags between the change event and its detection. Two studies carried out in tropical regions reported mean time lags of 63 days and 70 days for Myanmar and Bolivia respectively, if only optical Landsat 8 data is used [34,35].

At the Malawi site, we also tested if exploiting overlaps of two S-2 orbits delivers higher forest change detection accuracies. Using two orbits has the advantage of providing a much denser time series, but a potential disadvantage stems from the geometric shifts related to the different orbit viewing angles. Especially, the forest border is often slightly misplaced in data from two orbits as the DEM that is used by the European Space Agency (ESA) to orthorectify the imagery does not accurately account for tree height. At the 10 m spatial resolution of Sentinel-2, the effect of orbit-related geometric shifts on the spectral reflectance of single forest border pixels is much stronger than that, for example, in 30 m Landsat 8 pixels. Elevation errors in the DEM may also lead to strong dislocations of up to 20 m at mountain ridges as observed in the alpine test areas in Austria. In Malawi, our test site is characterized by flat to slightly hilly terrain, so errors related to topography are negligible. The results show that for Malawi, the users' accuracy increases from 63.5% to 71.1% when using two orbits instead of only using one orbit. Especially, the commission error is much higher with only one orbit (see mapped areas of 1407 ha versus 756 ha in Tables 3 and 4, respectively). The main reason for this strong improvement may be found by looking at the inter-annual temporal distribution of valid observations. The time series density at the Malawi test site is quite inhomogeneous (see Figure 6). While many clear observations are available during the dry season, the time series is sparse in the rainy season from November to April due to frequent cloud cover. Using two orbits greatly increases the chance of including a few clear observations during the rainy season, which is very important for a proper initialization of the time series model. To our knowledge, such a comparative analysis on using Sentinel-2 data from one or two orbits as input to a change detection approach has not been performed before. Thus, we cannot assess whether our findings are in line with other research results. While our findings suggest that denser time series resulting from orbit overlaps lead to higher forest change mapping accuracies, it is not yet clear if these findings for the dry forests in Malawi are also valid for other forest types and for other regions. Further studies with different data input scenarios and at various forest test sites are still needed to fully answer this question.

The presented approach for detecting forest changes in near real-time using Sentinel-2 imagery has a high potential for operationalization of forest monitoring services, such as improved and automated REDD+ (Reducing Emissions from Deforestation and Forest Degradation) services in the tropics, and windthrow damage assessment or bark beetle monitoring in Central Europe. Future studies to improve the presented forest change detection method and similar near real-time approaches should focus on an integrated and automated separation of different biotic and abiotic forest disturbance agents. Additional developments should include a multivariate analysis of relevant reflectance bands and a spatiotemporal analysis of changes. The near real-time capability should be tested in operational scenarios and the time lag of change detection needs to be assessed with field data or EO data of very high geometric and temporal resolution (e.g., Planet data). For tropical forest monitoring and REDD+ activities, the method should be tested in combination with recent Sentinel-1 SAR forest change detection approaches, such as SAR shadow detection [38], backscatter composite differencing [39], or SAR backscatter thresholding based on the coefficient of variation [40]. With minor adaptations, the presented time series analysis approach might also be directly applicable to SAR data. The approach could also be used for other EO-based applications, such as phenology monitoring in agriculture or grassland monitoring (detection of mowing events).

#### **5. Conclusions**

In this paper, we presented a new algorithm designed for vegetation monitoring and change detection using optical EO data. The approach is largely data-driven and designed especially for near real-time mapping in web- or cloud-based monitoring services. Compared to existing algorithms employing harmonic regression, this study explored methodological improvements, mainly in two aspects:

• Seasonal patterns related to plant phenology can vary strongly between years because of different climate conditions, such as temperature and rainfall. With widely used harmonic regression models, it can be difficult to separate these normal variations from true disturbances due to the deterministic modelling approach. Our algorithm uses structural time series models in state space form which take into account that the trend and seasonal components in a time series can evolve over time. This class of stochastic models is typically used in conjunction with the Kalman filter, which also enables elegant handling of irregular sampling intervals and signal artefacts like un-masked clouds and cloud shadows. We showed that it is possible to track the phenology-related vegetation dynamics more closely without losing the ability to detect disturbances.

• Many existing algorithms require that specific change probabilities or statistical boundaries have to be exceeded multiple times for a change to be confirmed. In our algorithm, the sequence of Kalman filter innovations (i.e., the differences between the one-step-ahead model forecasts and corresponding observations) opens up alternative possibilities to characterize change. If the time series model assumptions are correct and no changes occur, the innovation sequence should be mostly free of temporal autocorrelation and it should have zero mean. We used the CUSUM control chart to monitor these properties and assumed that significant deviations indicate change. Additionally, a separate anomaly test intended to suppress commission errors due to signal artefacts was implemented. Compared to the multiple confirmation approach, our results indicate that the CUSUM approach has the potential to reduce the time lag between change event and detection and a better chance of detecting subtle and gradual changes because of its cumulative nature. However, a better understanding of different change types and especially reference data for validation are required to support this claim.

Two challenging test sites located in Austria and Malawi were selected to test the algorithm in a forest change detection scenario based on Sentinel-2 data. Both sites show pronounced and dynamic seasonal patterns in the Sentinel-2 time series due to plant phenology. The dominant change types in Malawi are deforestation and forest degradation. In Austria, we were mainly interested in changes caused by storm damages or bark beetle infestations, but of course, deforestation also occurs. The validation of the results was performed based on visually interpreted points derived by a stratified random sampling approach. For the forest change class, we reported users' accuracies of 76.8% (Austria) to 88.6% (Malawi), and producers' accuracies of 68.2% (Malawi) to 80.4% (Austria). Due to the low rate of commission errors and large proportion of stable forest, overall accuracies reached over 99%.

In the Malawi site, we further showed that a denser time series with data from two different orbits results in better change detection results compared to using data from only one orbit. The larger number of input images seems to outweigh the possible negative effects of spectral and geometric differences related to the different viewing angles, which especially occur at forest edges. The observed increase in users' accuracy when using two orbits amounted to 7.6%. However, further studies with different data inputs and at various forest test sites are required to confirm these results.

In summary, it can be stated that the combination of structural time series models and Kalman filtering represents an appropriate method for a variety of automated forest monitoring applications. Beside the possibility of detecting abrupt changes, for example caused by storm damages or deforestation, the method also shows a high potential to detect more subtle and continuous changes, such as bark beetle infestations, forest degradation, or drought stress. Regarding the detection of insect infestations, the interesting research question to be dealt with in future is whether an early detection of bark beetle green attack is feasible.

**Author Contributions:** Conceptualization, M.P., J.D., and M.H.; methodology, M.P.; software, A.W.; formal analysis, M.P., J.D., A.W., and U.S.; investigation, M.P. and A.W.; data curation, M.P. and A.W.; writing—original draft preparation, M.P., J.D., and M.H.; writing—review and editing, M.S.; visualization, M.P. and M.H.; funding acquisition, M.H., J.D., and M.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the European Union's Horizon 2020 research and innovation programme, grant numbers 633464 (project DIABOLO) and 685761 (project EOMonDis), and the Austrian Research Promotion Agency (FFG), grant numbers 859764 (project AlpMon) and 878891 (project BEAT IT!).

**Acknowledgments:** We greatly appreciate the free, full, and open data policy adopted by ESA for the COPERNICUS programme. We would also like to thank the reviewers for providing valuable feedback which helped to improve the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).